/[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

revision 1797 by gross, Mon Jul 14 08:55:25 2008 UTC revision 1798 by gross, Wed Sep 17 06:21:12 2008 UTC
# Line 107  err_t Paso_Solver_GMRES2( Line 107  err_t Paso_Solver_GMRES2(
107               Paso_Update(n,1.,v[0],-1./normf0,f0);               Paso_Update(n,1.,v[0],-1./normf0,f0);
108               g[0]=normf0;               g[0]=normf0;
109            }            }
110            while(! (breakFlag || maxIterFlag || convergeFlag || (Status !=SOLVER_NO_ERROR) )) {            while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {
111                 k++;                 k++;
112                 v[k]=TMPMEMALLOC(n,double);                 v[k]=TMPMEMALLOC(n,double);
113                 if (v[k]==NULL) {                 if (v[k]==NULL) {
# Line 116  err_t Paso_Solver_GMRES2( Line 116  err_t Paso_Solver_GMRES2(
116                    /*                    /*
117                    *      call directional derivative function                    *      call directional derivative function
118                    */                    */
119                    Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work);                    Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,TRUE);
120                    normv=Paso_l2(n,v[k],F->mpi_info);                    normv=Paso_l2(n,v[k],F->mpi_info);
121                    /*                    /*
122                     * Modified Gram-Schmidt                     * Modified Gram-Schmidt
# Line 157  err_t Paso_Solver_GMRES2( Line 157  err_t Paso_Solver_GMRES2(
157                     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)]);
158                     if (nu>0) {                     if (nu>0) {
159                         c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;                         c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;
160                         s[k-1]=-h[k,k-1]/nu;                         s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
161                         h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[k,k-1];                         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)];
162                         h[INDEX2(k,k-1,l)]=0;                         h[INDEX2(k,k-1,l)]=0;
163                         ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));                         ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
164                     }                     }
# Line 172  err_t Paso_Solver_GMRES2( Line 172  err_t Paso_Solver_GMRES2(
172        /*        /*
173         * all done and ready for the forward substitution:         * all done and ready for the forward substitution:
174        */        */
175                for (i=k-1;i>=0;--i) {
176        for (i=0;i<k;i++) {             for (j=i+1;j<k;j++) {
177             for (j=0;j<i;j++) {                g[i]-=h[INDEX2(i,j,l)]*g[j];
               g[i]-=h[INDEX2(j,i,l)]*g[j];  
178             }             }
179             g[i]/=h[INDEX2(i,i,l)];             g[i]/=h[INDEX2(i,i,l)];
180             Paso_Update(n,1.,x,g[i],v[i]);             Paso_Update(n,1.,x,g[i],v[i]);

Legend:
Removed from v.1797  
changed lines
  Added in v.1798

  ViewVC Help
Powered by ViewVC 1.1.26