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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3490 - (show annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 3033 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

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,*/ 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 ttt=sqrt(EPSILON)*norm_w;
37 #pragma omp parallel private(local_s)
38 {
39 local_s=s;
40 #pragma omp for private(i, aw)
41 for (i=0;i<n;++i) {
42 aw=fabs(w[i]);
43 if ( aw>ttt ) {
44 local_s=MAX(local_s,fabs(x0[i])/aw);
45 }
46 }
47 #pragma omp critical
48 {
49 s=MAX(s,local_s);
50 }
51
52 }
53 #ifdef ESYS_MPI
54 {
55 double local_v[2], v[2];
56 local_v[0]=s;
57 local_v[1]=norm_w;
58 MPI_Allreduce(local_v,v, 2, MPI_DOUBLE, MPI_MAX, F->mpi_info->comm);
59 s=v[0];
60 norm_w=v[1];
61 }
62 #endif
63 /* printf("s :: = %e, %e \n",s, norm_x0/norm_w); */
64 if (norm_w>0) {
65 s=s*epsnew;
66 /* printf("s = %e\n",s); */
67 Paso_LinearCombination(n,setoff,1.,x0,s,w);
68 err=Paso_FunctionCall(F,J0w,setoff,pp);
69 if (err==SOLVER_NO_ERROR) {
70 Paso_Update(n,1./s,J0w,-1./s,f0); /* J0w = (J0w - f0)/epsnew; */
71 /* {
72 int i;
73 for (i=0;i<n; i++) printf("df[%d]=%e %e\n",i,J0w[i],w[i]);
74 }
75 */
76 }
77 } else {
78 Paso_zeroes(n,J0w);
79 }
80 return err;
81 }
82
83 /*
84 * sets value=F(arg)
85 *
86 */
87 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp)
88 {
89 if (F!=NULL) {
90 switch(F->kind) {
91 case LINEAR_SYSTEM:
92 return Paso_Function_LinearSystem_call(F, value, arg,pp);
93 break;
94 case FCT:
95 return Paso_FCTSolver_Function_call(F, value, arg, pp);
96 break;
97 default:
98 return SYSTEM_ERROR;
99 }
100 }
101 /* Added by PGH, assume a null pointe is an error */
102 return SYSTEM_ERROR;
103 }
104 /*
105 * clear Paso_Function
106 */
107 void Paso_Function_free(Paso_Function * F) {
108 if (F!=NULL) {
109
110 switch(F->kind) {
111 case LINEAR_SYSTEM:
112 Paso_Function_LinearSystem_free(F);
113 break;
114 case FCT:
115 Paso_FCTSolver_Function_free(F);
116 break;
117 default:
118 MEMFREE(F);
119 }
120 }
121 }

  ViewVC Help
Powered by ViewVC 1.1.26