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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1643 - (show annotations)
Tue Jul 15 05:24:14 2008 UTC (11 years, 6 months ago) by phornby
File MIME type: text/plain
File size: 2157 byte(s)
Missing include file. A

ssume that a NULL F passed to Paso_FunctionCall() is an error.
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 #include "Solver.h"
18 /*
19 * numerical calculation of the directional derivative J0w if F at x0 in the direction w. f0 is the value of F at x0.
20 * setoff is workspace
21 */
22
23 err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff)
24 {
25 err_t err=0;
26 dim_t n=F->n;
27 double norm_w,epsnew,norm_x0;
28 epsnew=10.*sqrt(EPSILON);
29 norm_w=Paso_l2(n,w,F->mpi_info);
30
31 if (norm_w>0) {
32 Paso_zeroes(n,J0w);
33 } else {
34 epsnew = epsnew/norm_w;
35 norm_x0=Paso_l2(n,x0,F->mpi_info);
36 if (norm_x0>0) epsnew*=norm_x0;
37 Paso_LinearCombination(n,setoff,1.,x0,epsnew,w);
38 err=Paso_FunctionCall(F,J0w,setoff);
39 if (err==NO_ERROR) {
40 Paso_Update(n,1./epsnew,J0w,-1./epsnew,f0); /* J0w = (J0w - f0)/epsnew; */
41 }
42 }
43 return err;
44 }
45
46 /*
47 * sets value=F(arg)
48 *
49 */
50 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg)
51 {
52 if (F!=NULL) {
53 switch(F->kind) {
54 case LINEAR_SYSTEM:
55 return Paso_Function_LinearSystem_call(F, value, arg);
56 case FCT:
57 return Paso_Function_FCT_call(F, value, arg);
58 default:
59 return SYSTEM_ERROR;
60 }
61 }
62 /* Added by PGH, assume a null pointe is an error */
63 return SYSTEM_ERROR;
64 }
65 /*
66 * clear Paso_Function
67 */
68 void Paso_Function_free(Paso_Function * F) {
69 if (F!=NULL) {
70
71 switch(F->kind) {
72 case LINEAR_SYSTEM:
73 Paso_Function_LinearSystem_free(F);
74 break;
75 case FCT:
76 Paso_Function_FCT_free(F);
77 break;
78 default:
79 MEMFREE(F);
80 }
81 }
82 }

  ViewVC Help
Powered by ViewVC 1.1.26