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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 8 months ago) by gross
File MIME type: text/plain
File size: 2864 byte(s)
new implementation of FCT solver with some modifications to the python interface
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 /*
20 * numerical calculation of the directional derivative J0w if F at x0 in the direction w. f0 is the value of F at x0.
21 * setoff is workspace
22 */
23
24 err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff, Paso_Performance* pp)
25 {
26 err_t err=SOLVER_NO_ERROR;
27 const dim_t n=F->n;
28 dim_t i;
29 register double aw;
30 const double epsnew=sqrt(EPSILON);
31 double /*norm_x0,*/ ttt, s=epsnew, local_s, norm_w=0.;
32
33 /*norm_x0=Paso_lsup(n,x0,F->mpi_info);*/
34 norm_w=Paso_lsup(n,w,F->mpi_info);
35 ttt=sqrt(EPSILON)*norm_w;
36 #pragma omp parallel private(local_s)
37 {
38 local_s=s;
39 #pragma omp for private(i, aw)
40 for (i=0;i<n;++i) {
41 aw=fabs(w[i]);
42 if ( aw>ttt ) {
43 local_s=MAX(local_s,fabs(x0[i])/aw);
44 }
45 }
46 #pragma omp critical
47 {
48 s=MAX(s,local_s);
49 }
50
51 }
52 #ifdef ESYS_MPI
53 {
54 double local_v[2], v[2];
55 local_v[0]=s;
56 local_v[1]=norm_w;
57 MPI_Allreduce(local_v,v, 2, MPI_DOUBLE, MPI_MAX, F->mpi_info->comm);
58 s=v[0];
59 norm_w=v[1];
60 }
61 #endif
62 /* printf("s :: = %e, %e \n",s, norm_x0/norm_w); */
63 if (norm_w>0) {
64 s=s*epsnew;
65 /* printf("s = %e\n",s); */
66 Paso_LinearCombination(n,setoff,1.,x0,s,w);
67 err=Paso_FunctionCall(F,J0w,setoff,pp);
68 if (err==SOLVER_NO_ERROR) {
69 Paso_Update(n,1./s,J0w,-1./s,f0); /* J0w = (J0w - f0)/epsnew; */
70 /* {
71 int i;
72 for (i=0;i<n; i++) printf("df[%d]=%e %e\n",i,J0w[i],w[i]);
73 }
74 */
75 }
76 } else {
77 Paso_zeroes(n,J0w);
78 }
79 return err;
80 }
81
82 /*
83 * sets value=F(arg)
84 *
85 */
86 err_t Paso_FunctionCall(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp)
87 {
88 if (F!=NULL) {
89 switch(F->kind) {
90 case LINEAR_SYSTEM:
91 return Paso_Function_LinearSystem_call(F, value, arg,pp);
92 break;
93 default:
94 return SYSTEM_ERROR;
95 }
96 }
97 /* Added by PGH, assume a null pointer is an error */
98 return SYSTEM_ERROR;
99 }
100
101 /*
102 * clear Paso_Function
103 */
104 void Paso_Function_free(Paso_Function * F) {
105 if (F!=NULL) {
106
107 switch(F->kind) {
108 case LINEAR_SYSTEM:
109 Paso_Function_LinearSystem_free(F);
110 break;
111 default:
112 MEMFREE(F);
113 }
114 }
115 }
116

  ViewVC Help
Powered by ViewVC 1.1.26