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

Revision 2987 - (show annotations)
Tue Mar 16 01:32:43 2010 UTC (9 years, 7 months ago) by gross
File MIME type: text/plain
File size: 2282 byte(s)
```FCT solver rewritten
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 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 15 #include "Common.h" 16 #include "Functions.h" 17 #include "PasoUtil.h" 18 #include "Solver.h" 19 #include "FCTSolver.h" 20 /* 21 * numerical calculation of the directional derivative J0w if F at x0 in the direction w. f0 is the value of F at x0. 22 * setoff is workspace 23 */ 24 25 err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff, Paso_Performance* pp) 26 { 27 err_t err=SOLVER_NO_ERROR; 28 dim_t n=F->n; 29 double norm_w,epsnew,norm_x0; 30 epsnew=10.*sqrt(EPSILON); 31 norm_w=Paso_lsup(n,w,F->mpi_info); 32 norm_x0=Paso_lsup(n,x0,F->mpi_info); 33 if (norm_w>0) { 34 epsnew = epsnew/norm_w; 35 if (norm_x0>0) epsnew*=norm_x0; 36 Paso_LinearCombination(n,setoff,1.,x0,epsnew,w); 37 err=Paso_FunctionCall(F,J0w,setoff,pp); 38 if (err==SOLVER_NO_ERROR) { 39 Paso_Update(n,1./epsnew,J0w,-1./epsnew,f0); /* J0w = (J0w - f0)/epsnew; */ 40 } 41 } else { 42 Paso_zeroes(n,J0w); 43 } 44 return err; 45 } 46 47 /* 48 * sets value=F(arg) 49 * 50 */ 51 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp) 52 { 53 if (F!=NULL) { 54 switch(F->kind) { 55 case LINEAR_SYSTEM: 56 return Paso_Function_LinearSystem_call(F, value, arg,pp); 57 break; 58 case FCT: 59 return Paso_FCTSolver_Function_call(F, value, arg, pp); 60 break; 61 default: 62 return SYSTEM_ERROR; 63 } 64 } 65 /* Added by PGH, assume a null pointe is an error */ 66 return SYSTEM_ERROR; 67 } 68 /* 69 * clear Paso_Function 70 */ 71 void Paso_Function_free(Paso_Function * F) { 72 if (F!=NULL) { 73 74 switch(F->kind) { 75 case LINEAR_SYSTEM: 76 Paso_Function_LinearSystem_free(F); 77 break; 78 case FCT: 79 Paso_FCTSolver_Function_free(F); 80 break; 81 default: 82 MEMFREE(F); 83 } 84 } 85 }