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

Diff of /trunk/paso/src/GMRES2.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 23  err_t Paso_Solver_GMRES2( Line 23  err_t Paso_Solver_GMRES2(
23      Paso_Function * F,      Paso_Function * F,
24      const double* f0,      const double* f0,
25      const double* x0,      const double* x0,
26      double * x,      double * dx,
27      dim_t *iter,      dim_t *iter,
28      double* tolerance,      double* tolerance,
29      Paso_Performance* pp)      Paso_Performance* pp)
# Line 42  err_t Paso_Solver_GMRES2( Line 42  err_t Paso_Solver_GMRES2(
42    if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {    if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
43      return SOLVER_INPUT_ERROR;      return SOLVER_INPUT_ERROR;
44    }    }
45    Paso_zeroes(n,x);    Paso_zeroes(n,dx);
46    /*        /*    
47     *  allocate memory:     *  allocate memory:
48     */     */
# Line 72  err_t Paso_Solver_GMRES2( Line 72  err_t Paso_Solver_GMRES2(
72               Status=SOLVER_MEMORY_ERROR;               Status=SOLVER_MEMORY_ERROR;
73            } else {            } else {
74               Paso_zeroes(n,v[0]);               Paso_zeroes(n,v[0]);
75               Paso_Update(n,1.,v[0],-1./normf0,f0);               Paso_Update(n,1.,v[0],-1./normf0,f0); /* v = -1./normf0*f0 */
76               g[0]=normf0;               g[0]=normf0;
77            }            }
78            while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {            while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {
# Line 91  err_t Paso_Solver_GMRES2( Line 91  err_t Paso_Solver_GMRES2(
91                     */                     */
92                     for (j=0;j<k;j++){                     for (j=0;j<k;j++){
93                        hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);                        hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
94                        Paso_Update(n,1.,v[k],(-hh),v[j]);                        Paso_Update(n,1.,v[k],(-hh),v[j]); /* v[k]-hh*v[j] */
95                        h[INDEX2(j,k-1,l)]=hh;                        h[INDEX2(j,k-1,l)]=hh;
96                     }   /* printf("%d :  %d = %e\n",k,j,hh); */
97               }
98                     normv2=Paso_l2(n,v[k],F->mpi_info);                     normv2=Paso_l2(n,v[k],F->mpi_info);
99                     h[INDEX2(k,k-1,l)]=normv2;                     h[INDEX2(k,k-1,l)]=normv2;
100                     /*                     /*
101                      * reorthogonalize                      * reorthogonalize
102                      */                      */
103                     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {                     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {
104    /* printf("GMRES2: renormalization!"); */
105                          for (j=0;j<k;j++){                          for (j=0;j<k;j++){
106                              hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);                              hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
107    
108                          h[INDEX2(j,k-1,l)]+=hr;                          h[INDEX2(j,k-1,l)]+=hr;
109                              Paso_Update(n,1.,v[k],(-hr),v[j]);                              Paso_Update(n,1.,v[k],(-hr),v[j]);
110                          }                          }
111                          normv2=Paso_l2(n,v[k],F->mpi_info);                          normv2=Paso_l2(n,v[k],F->mpi_info);
112                          h[INDEX2(k,k-1,l)]=normv2;                          h[INDEX2(k,k-1,l)]=normv2;
113                      }                      }
114    {int p,q;
115       for (p=0;p<k+1;p++){
116       for (q=0;q<k+1;q++)printf("%e ",h[INDEX2(p,q,l)]);
117        printf("\n");
118      
119       }
120    }
121                     /*                     /*
122                      *   watch out for happy breakdown                      *   watch out for happy breakdown
123                      */                      */
124                     if(normv2 > 0.) {                     if(normv2 > 0.) {
125                        Paso_Update(n,1./normv2,v[k],0.,v[k]);                        Paso_Update(n,1./normv2,v[k],0.,v[k]); /* normalize v[k] */
126                     }                     }
127                     /*                     /*
128                      *   Form and store the information for the new Givens rotation                      *   Form and store the information for the new Givens rotation
129                      */                      */
130                     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);                     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
131    
132                     /*                     /*
133                      *  Don't divide by zero if solution has  been found                      *  Don't divide by zero if solution has  been found
134                      */                      */
135                     g[k]=0;                     g[k]=0;
136                     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)]);             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)]);
137                     if (nu>0) {                     if (nu>0) {
138                         c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;                         c[k-1]= h[INDEX2(k-1,k-1,l)]/nu;
139                         s[k-1]=-h[INDEX2(k,k-1,l)]/nu;                         s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
140                         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)];                         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)];
141                         h[INDEX2(k,k-1,l)]=0;                         h[INDEX2(k,k-1,l)]=0;
# Line 132  err_t Paso_Solver_GMRES2( Line 143  err_t Paso_Solver_GMRES2(
143                     }                     }
144                     norm_of_residual=fabs(g[k]);                     norm_of_residual=fabs(g[k]);
145                     maxIterFlag = (k>=iter_max);                     maxIterFlag = (k>=iter_max);
146                     convergeFlag = (fabs(g[k]) <= abs_tol);                     convergeFlag = (norm_of_residual <= abs_tol);
147                     printf("GMRES2 step %d: residual %e (abs. tol=%e)\n",k,fabs(g[k]),abs_tol);                     printf("GMRES2 step %d: residual %e (abs. tol=%e)\n",k,fabs(g[k]),abs_tol);
148                }                }
149            }            }
# Line 140  err_t Paso_Solver_GMRES2( Line 151  err_t Paso_Solver_GMRES2(
151        /*        /*
152         * all done and ready for the forward substitution:         * all done and ready for the forward substitution:
153        */        */
154    printf("k =  %d\n", k);
155    
156        for (i=k-1;i>=0;--i) {        for (i=k-1;i>=0;--i) {
157             for (j=i+1;j<k;j++) {             for (j=i+1;j<k;j++) {
158    printf("%d %d : %e\n", i,j,h[INDEX2(i,j,l)]);
159                g[i]-=h[INDEX2(i,j,l)]*g[j];                g[i]-=h[INDEX2(i,j,l)]*g[j];
160             }             }
161             g[i]/=h[INDEX2(i,i,l)];             g[i]/=h[INDEX2(i,i,l)];
162             Paso_Update(n,1.,x,g[i],v[i]);             Paso_Update(n,1.,dx,g[i],v[i]); /* dx = dx+g[i]*v[i] */
163        }        }
164   }   }
165   /*       /*    

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

  ViewVC Help
Powered by ViewVC 1.1.26