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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1798 - (hide annotations)
Wed Sep 17 06:21:12 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/plain
File size: 2255 byte(s)
Fixes for the JacobeanFreeNewton scheme. Still needs to be tested under OPENMP but runs under MPI.


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

  ViewVC Help
Powered by ViewVC 1.1.26