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

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;ittt ) { 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;ikind) { 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 }