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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3005 - (show annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 5 months ago) by gross
File MIME type: text/plain
File size: 3048 byte(s)
early call of setPreconditioner in FCT solver removed.
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 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 norm_w=Paso_lsup(n,w,F->mpi_info);
36 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 if (norm_w>0) {
66 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 }
78 } else {
79 Paso_zeroes(n,J0w);
80 }
81 return err;
82 }
83
84 /*
85 * sets value=F(arg)
86 *
87 */
88 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp)
89 {
90 if (F!=NULL) {
91 switch(F->kind) {
92 case LINEAR_SYSTEM:
93 return Paso_Function_LinearSystem_call(F, value, arg,pp);
94 break;
95 case FCT:
96 return Paso_FCTSolver_Function_call(F, value, arg, pp);
97 break;
98 default:
99 return SYSTEM_ERROR;
100 }
101 }
102 /* Added by PGH, assume a null pointe is an error */
103 return SYSTEM_ERROR;
104 }
105 /*
106 * clear Paso_Function
107 */
108 void Paso_Function_free(Paso_Function * F) {
109 if (F!=NULL) {
110
111 switch(F->kind) {
112 case LINEAR_SYSTEM:
113 Paso_Function_LinearSystem_free(F);
114 break;
115 case FCT:
116 Paso_FCTSolver_Function_free(F);
117 break;
118 default:
119 MEMFREE(F);
120 }
121 }
122 }

  ViewVC Help
Powered by ViewVC 1.1.26