# Diff of /trunk/paso/src/GMRES2.c

trunk/paso/src/FGMRES.c revision 1624 by phornby, Wed Apr 23 05:31:54 2008 UTC trunk/paso/src/GMRES2.c revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 14  Line 14
14  *  Purpose  *  Purpose
15  *  =======  *  =======
16  *  *
17  *  FGMRES solves the non-linear system f0+J_0*d=0  *  GMRES2 solves the non-linear system f0+J_0*d=0
18  *  where f0=F(x0), J_0 is the jacobian of F at x0.  *  where f0=F(x0), J_0 is the jacobian of F at x0.
19  *  *
20  *  Convergence test: norm(f0+J_0*d)<=tolerance*norm(f0)  *  Convergence test: norm(f0+J_0*d)<=tolerance*norm(f0)
# Line 49  Line 49
49
50  #include "Common.h"  #include "Common.h"
51  #include "Solver.h"  #include "Solver.h"
52  #include "Util.h"  #include "PasoUtil.h"
53  #ifdef _OPENMP  #ifdef _OPENMP
54  #include <omp.h>  #include <omp.h>
55  #endif  #endif
56
57  err_t Paso_Solver_NLGMRES(  err_t Paso_Solver_GMRES2(
58      Paso_Function * F,      Paso_Function * F,
59      const double* f0,      const double* f0,
60      const double* x0,      const double* x0,
# Line 64  err_t Paso_Solver_NLGMRES( Line 64  err_t Paso_Solver_NLGMRES(
64      Paso_Performance* pp)      Paso_Performance* pp)
65  {  {
66    double static RENORMALIZATION_CONST=0.001;    double static RENORMALIZATION_CONST=0.001;
67    dim_t l=(*iter)+1, iter_max=*iter,k=0,n=F->local_n,i,j;    dim_t l=(*iter)+1, iter_max=*iter,k=0,i,j;
68      const dim_t n=F->n;
69    double rel_tol=*tolerance, abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;    double rel_tol=*tolerance, abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
70    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
71    double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;    double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;
# Line 87  err_t Paso_Solver_NLGMRES( Line 88  err_t Paso_Solver_NLGMRES(
88
89    if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {    if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {
90       Status=SOLVER_MEMORY_ERROR;       Status=SOLVER_MEMORY_ERROR;
} else {
for (i=0;i<iter_max;i++) v[i]=NULL;
for (i=0;i<iter_max;i++) {
v[i]=TMPMEMALLOC(n,double);
if (v[i]==NULL) {
Status=SOLVER_MEMORY_ERROR;
break;
}
}
91    }    }
92    if ( Status ==SOLVER_NO_ERROR ) {    if ( Status ==SOLVER_NO_ERROR ) {
93         for (i=0;i<iter_max;i++) v[i]=NULL;
94       /*           /*
95        *  the show begins:        *  the show begins:
96        */        */
# Line 106  err_t Paso_Solver_NLGMRES( Line 99  err_t Paso_Solver_NLGMRES(
99        convergeFlag=(ABS(normf0)<=0);        convergeFlag=(ABS(normf0)<=0);
100        if (! convergeFlag) {        if (! convergeFlag) {
101            abs_tol=rel_tol*normf0;            abs_tol=rel_tol*normf0;
102            Paso_zeroes(n,v[0]);            v[0]=TMPMEMALLOC(n,double);
103            Paso_Update(n,1.,v[0],-1./normf0,f0);            if (v[0]==NULL) {
104            g[0]=normf0;               Status=SOLVER_MEMORY_ERROR;
105            while(! (breakFlag || maxIterFlag || convergeFlag)) {            } else {
106                 Paso_zeroes(n,v[0]);
107                 Paso_Update(n,1.,v[0],-1./normf0,f0);
108                 g[0]=normf0;
109              }
110              while(! (breakFlag || maxIterFlag || convergeFlag || (Status !=SOLVER_NO_ERROR) )) {
111                 k++;                 k++;
112                 /*                 v[k]=TMPMEMALLOC(n,double);
113                 *      call directional derivative function                 if (v[k]==NULL) {
114                 */                    Status=SOLVER_MEMORY_ERROR;
115                 Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work);                 } else {
116                 normv=Paso_l2(n,v[k],F->mpi_info);                    /*
117                 /*                    *      call directional derivative function
118                  * Modified Gram-Schmidt                    */
119                  */                    Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work);
120                  for (j=0;j<k;j++){                    normv=Paso_l2(n,v[k],F->mpi_info);
121                     hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);                    /*
122                     Paso_Update(n,1.,v[k],(-hh),v[j]);                     * Modified Gram-Schmidt
123                     h[INDEX2(j,k-1,l)]=hh;                     */
124                  }                     for (j=0;j<k;j++){
125                  normv2=Paso_l2(n,v[k],F->mpi_info);                        hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
126                  h[INDEX2(k,k-1,l)]=normv2;                        Paso_Update(n,1.,v[k],(-hh),v[j]);
127                  /*                        h[INDEX2(j,k-1,l)]=hh;
128                   * reorthogonalize                     }
129                   */                     normv2=Paso_l2(n,v[k],F->mpi_info);
130                  if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {                     h[INDEX2(k,k-1,l)]=normv2;
131                       for (j=0;j<k;j++){                     /*
132                           hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);                      * reorthogonalize
133                       h[INDEX2(j,k-1,l)]+=hr;                      */
134                           Paso_Update(n,1.,v[k],(-hr),v[j]);                     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {
135                       }                          for (j=0;j<k;j++){
136                       normv2=Paso_l2(n,v[k],F->mpi_info);                              hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
137                       h[INDEX2(k,k-1,l)]=normv2;                          h[INDEX2(j,k-1,l)]+=hr;
138                   }                              Paso_Update(n,1.,v[k],(-hr),v[j]);
139                  /*                          }
140                   *   watch out for happy breakdown                          normv2=Paso_l2(n,v[k],F->mpi_info);
141                   */                          h[INDEX2(k,k-1,l)]=normv2;
142                  if(normv2 > 0.) {                      }
143                     Paso_Update(n,1./normv2,v[k],0.,v[k]);                     /*
144                  }                      *   watch out for happy breakdown
145                  /*                      */
146                   *   Form and store the information for the new Givens rotation                     if(normv2 > 0.) {
147                   */                        Paso_Update(n,1./normv2,v[k],0.,v[k]);
148                  ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);                     }
149                  /*                     /*
150                   *  Don't divide by zero if solution has  been found                      *   Form and store the information for the new Givens rotation
151                   */                      */
152                  g[k]=0;                     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
153                  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)]);                     /*
154                  if (nu>0) {                      *  Don't divide by zero if solution has  been found
155                      c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;                      */
156                      s[k-1]=-h[k,k-1]/nu;                     g[k]=0;
157                      h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[k,k-1];                     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)]);
158                      h[INDEX2(k,k-1,l)]=0;                     if (nu>0) {
159                      ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));                         c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;
160                  }                         s[k-1]=-h[k,k-1]/nu;
161                  norm_of_residual=fabs(g[k]);                         h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[k,k-1];
162                  maxIterFlag = (k>=iter_max);                         h[INDEX2(k,k-1,l)]=0;
163                  convergeFlag = (fabs(g[k]) <= abs_tol);                         ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
164                  printf("FGMRES step %d: error %e (tol=%e)\n",k,fabs(g[k]),abs_tol);                     }
165                       norm_of_residual=fabs(g[k]);
166                       maxIterFlag = (k>=iter_max);
167                       convergeFlag = (fabs(g[k]) <= abs_tol);
168                       printf("GMRES2 step %d: error %e (tol=%e)\n",k,fabs(g[k]),abs_tol);
169                  }
170            }            }
171        }        }
172        /*        /*
# Line 182  err_t Paso_Solver_NLGMRES( Line 185  err_t Paso_Solver_NLGMRES(
185    *  clean up:    *  clean up:
186    */    */
187    if ( v !=NULL) {    if ( v !=NULL) {
188      for (i=0;i<iter_max;i++) TMPMEMFREE(v);      for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);
189    }    }
190    TMPMEMFREE(h);    TMPMEMFREE(h);
191    TMPMEMFREE(v);    TMPMEMFREE(v);

Legend:
 Removed from v.1624 changed lines Added in v.1639