/[escript]/trunk/paso/src/Functions.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 }

  ViewVC Help
Powered by ViewVC 1.1.26