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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3005 - (hide annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/plain
File size: 3048 byte(s)
early call of setPreconditioner in FCT solver removed.
1 gross 1476
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 gross 1476
14 ksteube 1811
15 gross 1476 #include "Common.h"
16     #include "Functions.h"
17 gross 1639 #include "PasoUtil.h"
18 phornby 1643 #include "Solver.h"
19 gross 2987 #include "FCTSolver.h"
20 gross 1476 /*
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 gross 2987 err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff, Paso_Performance* pp)
26 gross 1476 {
27 gross 2987 err_t err=SOLVER_NO_ERROR;
28 gross 3005 const dim_t n=F->n;
29     dim_t i;
30     register double aw;
31     const double epsnew=sqrt(EPSILON);
32     double norm_x0, tt, ttt, s=epsnew, local_s, norm_w=0.;
33    
34     norm_x0=Paso_lsup(n,x0,F->mpi_info);
35 gross 2987 norm_w=Paso_lsup(n,w,F->mpi_info);
36 gross 3005 tt=EPSILON*norm_x0;
37     ttt=sqrt(EPSILON)*norm_w;
38     #pragma omp parallel private(local_s, local_norm_w)
39     {
40     local_s=s;
41     #pragma omp for private(i, aw)
42     for (i=0;i<n;++i) {
43     aw=fabs(w[i]);
44     if ( aw>ttt ) {
45     local_s=MAX(local_s,fabs(x0[i])/aw);
46     }
47     }
48     #pragma omp critical
49     {
50     s=MAX(s,local_s);
51     }
52    
53     }
54     #ifdef PASO_MPI
55     {
56     double local_v[2], v[2];
57     local_v[0]=s;
58     local_v[1]=norm_w;
59     MPI_Allreduce(local_v,v, 2, MPI_DOUBLE, MPI_MAX, F->mpi_info);
60     s=v[0];
61     norm_w=v[1];
62     }
63     #endif
64     printf("s :: = %e, %e \n",s, norm_x0/norm_w);
65 gross 1476 if (norm_w>0) {
66 gross 3005 s=s*epsnew;
67     printf("s = %e\n",s);
68     Paso_LinearCombination(n,setoff,1.,x0,s,w);
69     err=Paso_FunctionCall(F,J0w,setoff,pp);
70     if (err==SOLVER_NO_ERROR) {
71     Paso_Update(n,1./s,J0w,-1./s,f0); /* J0w = (J0w - f0)/epsnew; */
72     /* {
73     int i;
74     for (i=0;i<n; i++) printf("df[%d]=%e %e\n",i,J0w[i],w[i]);
75     }
76     */
77 gross 2987 }
78 gross 1798 } else {
79     Paso_zeroes(n,J0w);
80 gross 1476 }
81     return err;
82     }
83    
84 gross 1639 /*
85     * sets value=F(arg)
86     *
87     */
88 gross 1804 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp)
89 gross 1476 {
90 gross 1639 if (F!=NULL) {
91     switch(F->kind) {
92     case LINEAR_SYSTEM:
93 gross 1804 return Paso_Function_LinearSystem_call(F, value, arg,pp);
94 gross 2987 break;
95     case FCT:
96     return Paso_FCTSolver_Function_call(F, value, arg, pp);
97     break;
98 gross 1639 default:
99     return SYSTEM_ERROR;
100     }
101     }
102 phornby 1643 /* Added by PGH, assume a null pointe is an error */
103     return SYSTEM_ERROR;
104 gross 1639 }
105     /*
106     * clear Paso_Function
107     */
108     void Paso_Function_free(Paso_Function * F) {
109     if (F!=NULL) {
110 gross 1476
111 gross 1639 switch(F->kind) {
112     case LINEAR_SYSTEM:
113     Paso_Function_LinearSystem_free(F);
114     break;
115 gross 2987 case FCT:
116     Paso_FCTSolver_Function_free(F);
117     break;
118 gross 1639 default:
119     MEMFREE(F);
120     }
121     }
122 gross 1476 }

  ViewVC Help
Powered by ViewVC 1.1.26