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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3004 by gross, Tue Mar 16 01:32:43 2010 UTC revision 3005 by gross, Thu Apr 22 05:59:31 2010 UTC
# Line 25  Line 25 
25  err_t Paso_FunctionDerivative(double* J0w, const double* w, Paso_Function* F, const double *f0, const double *x0, double* setoff, Paso_Performance* pp)  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;     err_t err=SOLVER_NO_ERROR;
28     dim_t n=F->n;     const dim_t n=F->n;
29     double norm_w,epsnew,norm_x0;     dim_t i;
30     epsnew=10.*sqrt(EPSILON);     register double aw;
31     norm_w=Paso_lsup(n,w,F->mpi_info);     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);     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) {     if (norm_w>0) {
66         epsnew = epsnew/norm_w;          s=s*epsnew;
67             if (norm_x0>0) epsnew*=norm_x0;  printf("s = %e\n",s);
68             Paso_LinearCombination(n,setoff,1.,x0,epsnew,w);          Paso_LinearCombination(n,setoff,1.,x0,s,w);
69             err=Paso_FunctionCall(F,J0w,setoff,pp);          err=Paso_FunctionCall(F,J0w,setoff,pp);
70             if (err==SOLVER_NO_ERROR) {          if (err==SOLVER_NO_ERROR) {
71                Paso_Update(n,1./epsnew,J0w,-1./epsnew,f0); /* J0w = (J0w - f0)/epsnew; */                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 {     } else {
79         Paso_zeroes(n,J0w);         Paso_zeroes(n,J0w);

Legend:
Removed from v.3004  
changed lines
  Added in v.3005

  ViewVC Help
Powered by ViewVC 1.1.26