# Contents of /trunk/paso/src/Functions.c

Revision 1639 - (show annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/plain
File size: 2057 byte(s)

 1 /* \$Id:\$ */ 2 3 /******************************************************* 4 * 5 * Copyright 2008 by University of Queensland 6 * 7 * http://esscc.uq.edu.au 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 * http://www.opensource.org/licenses/osl-3.0.php 11 * 12 *******************************************************/ 13 14 #include "Common.h" 15 #include "Functions.h" 16 #include "PasoUtil.h" 17 /* 18 * numerical calculation of the directional derivative J0w if F at x0 in the direction w. f0 is the value of F at x0. 19 * setoff is workspace 20 */ 21 22 err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff) 23 { 24 err_t err=0; 25 dim_t n=F->n; 26 double norm_w,epsnew,norm_x0; 27 epsnew=10.*sqrt(EPSILON); 28 norm_w=Paso_l2(n,w,F->mpi_info); 29 30 if (norm_w>0) { 31 Paso_zeroes(n,J0w); 32 } else { 33 epsnew = epsnew/norm_w; 34 norm_x0=Paso_l2(n,x0,F->mpi_info); 35 if (norm_x0>0) epsnew*=norm_x0; 36 Paso_LinearCombination(n,setoff,1.,x0,epsnew,w); 37 err=Paso_FunctionCall(F,J0w,setoff); 38 if (err==NO_ERROR) { 39 Paso_Update(n,1./epsnew,J0w,-1./epsnew,f0); /* J0w = (J0w - f0)/epsnew; */ 40 } 41 } 42 return err; 43 } 44 45 /* 46 * sets value=F(arg) 47 * 48 */ 49 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg) 50 { 51 if (F!=NULL) { 52 switch(F->kind) { 53 case LINEAR_SYSTEM: 54 return Paso_Function_LinearSystem_call(F, value, arg); 55 case FCT: 56 return Paso_Function_FCT_call(F, value, arg); 57 default: 58 return SYSTEM_ERROR; 59 } 60 } 61 } 62 /* 63 * clear Paso_Function 64 */ 65 void Paso_Function_free(Paso_Function * F) { 66 if (F!=NULL) { 67 68 switch(F->kind) { 69 case LINEAR_SYSTEM: 70 Paso_Function_LinearSystem_free(F); 71 break; 72 case FCT: 73 Paso_Function_FCT_free(F); 74 break; 75 default: 76 MEMFREE(F); 77 } 78 } 79 }