/[escript]/trunk/paso/src/GMRES2.cpp
ViewVC logotype

Diff of /trunk/paso/src/GMRES2.cpp

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

branches/doubleplusgood/paso/src/GMRES2.c revision 4257 by jfenwick, Wed Feb 27 03:42:40 2013 UTC trunk/paso/src/GMRES2.cpp revision 4856 by caltinay, Wed Apr 9 06:46:55 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
# Line 17  Line 18 
18  #include "Common.h"  #include "Common.h"
19  #include "Solver.h"  #include "Solver.h"
20  #include "PasoUtil.h"  #include "PasoUtil.h"
21  #ifdef _OPENMP  
22  #include <omp.h>  namespace paso {
23  #endif  
24    err_t Solver_GMRES2(Function* F, const double* f0, const double* x0,
25  err_t Paso_Solver_GMRES2(                      double* dx, dim_t* iter, double* tolerance,
26      Paso_Function * F,                      Paso_Performance* pp)
     const double* f0,  
     const double* x0,  
     double * dx,  
     dim_t *iter,  
     double* tolerance,  
     Paso_Performance* pp)  
27  {  {
28    static double RENORMALIZATION_CONST=0.001;      static double RENORMALIZATION_CONST=0.001;
29    const dim_t l=(*iter)+1, iter_max=*iter;      const dim_t l=(*iter)+1, iter_max=*iter;
30    dim_t k=0,i,j;      dim_t k=0,i,j;
31    const dim_t n=F->n;      const dim_t n=F->getLen();
32    const double rel_tol=*tolerance;      const double rel_tol=*tolerance;
33    double abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;      double abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
34    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;      bool breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
35    double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;  
36    err_t Status=SOLVER_NO_ERROR;      if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
37            return SOLVER_INPUT_ERROR;
38    /*     Test the input parameters. */      }
39    if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {  
40      return SOLVER_INPUT_ERROR;      err_t Status=SOLVER_NO_ERROR;
41    }  
42    Paso_zeroes(n,dx);      double* h = new double[l*l];
43    /*          double** v = new double*[l];
44     *  allocate memory:      double* c = new double[l];
45     */      double* s = new double[l];
46    h=TMPMEMALLOC(l*l,double);      double* g = new double[l];
47    v=TMPMEMALLOC(l,double*);      double* work = new double[n];
48    c=TMPMEMALLOC(l,double);  
49    s=TMPMEMALLOC(l,double);      for (i=0; i<iter_max; i++)
50    g=TMPMEMALLOC(l,double);          v[i]=NULL;
51    work=TMPMEMALLOC(n,double);  
52        util::zeroes(n,dx);
53    if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {  
54       Status=SOLVER_MEMORY_ERROR;      /*    
55    }       *  the show begins:
56    if ( Status ==SOLVER_NO_ERROR ) {       */
57       for (i=0;i<iter_max;i++) v[i]=NULL;      normf0 = util::l2(n,f0,F->mpi_info);
58       /*          k = 0;
59        *  the show begins:      convergeFlag = (ABS(normf0)<=0);
60        */      if (!convergeFlag) {
61        normf0=Paso_l2(n,f0,F->mpi_info);          abs_tol = rel_tol*normf0;
62        k=0;          printf("GMRES2 initial residual norm %e (rel. tol=%e)\n",normf0,rel_tol);
63        convergeFlag=(ABS(normf0)<=0);          v[0] = new double[n];
64        if (! convergeFlag) {          util::zeroes(n, v[0]);
65            abs_tol=rel_tol*normf0;          util::update(n, 1., v[0], -1./normf0, f0); // v = -1./normf0*f0
66            printf("GMRES2 initial residual norm %e (rel. tol=%e)\n",normf0,rel_tol);          g[0] = normf0;
67            v[0]=TMPMEMALLOC(n,double);          while (!breakFlag && !maxIterFlag && !convergeFlag && Status==SOLVER_NO_ERROR) {
68            if (v[0]==NULL) {              k++;
69               Status=SOLVER_MEMORY_ERROR;              v[k]=new double[n];
70            } else {              /*
71               Paso_zeroes(n,v[0]);               * call directional derivative function
72               Paso_Update(n,1.,v[0],-1./normf0,f0); /* v = -1./normf0*f0 */               */
73               g[0]=normf0;              F->derivative(v[k], v[k-1], f0, x0, work, pp);
74            }              normv=util::l2(n,v[k],F->mpi_info);
75            while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {              /*
76                 k++;               * Modified Gram-Schmidt
77                 v[k]=TMPMEMALLOC(n,double);               */
78                 if (v[k]==NULL) {              for (j=0; j<k; j++){
79                    Status=SOLVER_MEMORY_ERROR;                  hh = util::innerProduct(n,v[j],v[k],F->mpi_info);
80                 } else {                  util::update(n,1.,v[k],(-hh),v[j]); // v[k]-hh*v[j]
81                    /*                  h[INDEX2(j,k-1,l)]=hh;
82                    *      call directional derivative function                  //printf("%d :  %d = %e\n",k,j,hh);
83                    */              }
84                    Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,pp);              normv2=util::l2(n,v[k],F->mpi_info);
85                    normv=Paso_l2(n,v[k],F->mpi_info);              h[INDEX2(k,k-1,l)]=normv2;
86                    /*              /*
87                     * Modified Gram-Schmidt               * reorthogonalize
88                     */               */
89                     for (j=0;j<k;j++){              if (!(normv + RENORMALIZATION_CONST*normv2 > normv)) {
90                        hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);                  // printf("GMRES2: renormalization!");
91                        Paso_Update(n,1.,v[k],(-hh),v[j]); /* v[k]-hh*v[j] */                  for (j=0; j<k; j++) {
92                        h[INDEX2(j,k-1,l)]=hh;                      hr = util::innerProduct(n,v[j],v[k],F->mpi_info);
93   /* printf("%d :  %d = %e\n",k,j,hh); */  
94             }                      h[INDEX2(j,k-1,l)]+=hr;
95                     normv2=Paso_l2(n,v[k],F->mpi_info);                      util::update(n,1.,v[k],(-hr),v[j]);
96                     h[INDEX2(k,k-1,l)]=normv2;                  }
97                     /*                  normv2=util::l2(n,v[k],F->mpi_info);
98                      * reorthogonalize                  h[INDEX2(k,k-1,l)]=normv2;
99                      */              }
100                     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {              /*
101  /* printf("GMRES2: renormalization!"); */               * watch out for happy breakdown
102                          for (j=0;j<k;j++){               */
103                              hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);              if (normv2 > 0.) {
104                    util::update(n,1./normv2,v[k],0.,v[k]); /* normalize v[k] */
105                          h[INDEX2(j,k-1,l)]+=hr;              }
106                              Paso_Update(n,1.,v[k],(-hr),v[j]);              /*
107                          }               * Form and store the information for the new Givens rotation
108                          normv2=Paso_l2(n,v[k],F->mpi_info);               */
109                          h[INDEX2(k,k-1,l)]=normv2;              util::applyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
110                      }  
111  /*              /*
112  {int p,q;               * Don't divide by zero if solution has been found
113     for (p=0;p<k+1;p++){               */
114     for (q=0;q<k+1;q++)printf("%e ",h[INDEX2(p,q,l)]);              g[k] = 0;
115      printf("\n");              nu = sqrt(h[INDEX2(k-1,k-1,l)]*h[INDEX2(k-1,k-1,l)]+h[INDEX2(k,k-1,l)]*h[INDEX2(k,k-1,l)]);
116                  if (nu > 0) {
117     }                  c[k-1]= h[INDEX2(k-1,k-1,l)]/nu;
118  }                  s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
119  */                  h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[INDEX2(k,k-1,l)];
120                     /*                  h[INDEX2(k,k-1,l)]=0;
121                      *   watch out for happy breakdown                  util::applyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
122                      */              }
123                     if(normv2 > 0.) {              norm_of_residual = fabs(g[k]);
124                        Paso_Update(n,1./normv2,v[k],0.,v[k]); /* normalize v[k] */              maxIterFlag = (k >= iter_max);
125                     }              convergeFlag = (norm_of_residual <= abs_tol);
126                     /*              printf("GMRES2 step %d: residual %e (abs. tol=%e)\n", k, fabs(g[k]), abs_tol);
127                      *   Form and store the information for the new Givens rotation          }
128                      */      }
129                     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);  
130        // all done and ready for the forward substitution:
131                     /*      for (i=k-1;i>=0;--i) {
132                      *  Don't divide by zero if solution has  been found          for (j=i+1;j<k;j++) {
133                      */              g[i]-=h[INDEX2(i,j,l)]*g[j];
134                     g[k]=0;          }
135             nu=sqrt(h[INDEX2(k-1,k-1,l)]*h[INDEX2(k-1,k-1,l)]+h[INDEX2(k,k-1,l)]*h[INDEX2(k,k-1,l)]);          g[i] /= h[INDEX2(i,i,l)];
136                     if (nu>0) {          util::update(n, 1., dx, g[i], v[i]); // dx = dx+g[i]*v[i]
137                         c[k-1]= h[INDEX2(k-1,k-1,l)]/nu;      }
138                         s[k-1]=-h[INDEX2(k,k-1,l)]/nu;      if ( v != NULL) {
139                         h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[INDEX2(k,k-1,l)];          for (i=0; i<iter_max; i++)
140                         h[INDEX2(k,k-1,l)]=0;              delete[] v[i];
141                         ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));      }
142                     }      delete[] h;
143                     norm_of_residual=fabs(g[k]);      delete[] v;
144                     maxIterFlag = (k>=iter_max);      delete[] c;
145                     convergeFlag = (norm_of_residual <= abs_tol);      delete[] s;
146                     printf("GMRES2 step %d: residual %e (abs. tol=%e)\n",k,fabs(g[k]),abs_tol);      delete[] g;
147                }      delete[] work;
148            }      *iter=k;
149        }      *tolerance=norm_of_residual;
150        /*      return Status;
        * all done and ready for the forward substitution:  
       */  
   
       for (i=k-1;i>=0;--i) {  
            for (j=i+1;j<k;j++) {  
               g[i]-=h[INDEX2(i,j,l)]*g[j];  
            }  
            g[i]/=h[INDEX2(i,i,l)];  
            Paso_Update(n,1.,dx,g[i],v[i]); /* dx = dx+g[i]*v[i] */  
       }  
  }  
  /*      
   *  clean up:  
   */  
   if ( v !=NULL) {  
     for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);  
   }  
   TMPMEMFREE(h);  
   TMPMEMFREE(v);  
   TMPMEMFREE(c);  
   TMPMEMFREE(s);  
   TMPMEMFREE(g);  
   TMPMEMFREE(work);  
   *iter=k;  
   *tolerance=norm_of_residual;  
   return Status;  
151  }  }
152    
153    } // namespace paso
154    

Legend:
Removed from v.4257  
changed lines
  Added in v.4856

  ViewVC Help
Powered by ViewVC 1.1.26