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

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

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

trunk/esys2/paso/src/Solvers/GMRES.c revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/paso/src/Solvers/GMRES.c revision 415 by gross, Wed Jan 4 05:37:33 2006 UTC
# Line 14  Line 14 
14  *  Arguments  *  Arguments
15  *  =========  *  =========
16  *  *
17  *  r       (input/output) double array, dimension n_row.  *  r       (input/output) double array, dimension n.
18  *          On entry, residual of inital guess X  *          On entry, residual of inital guess X
19  *  *
20  *  x       (input/output) double array, dimension n_col.  *  x       (input/output) double array, dimension n.
21  *          On input, the initial guess.  *          On input, the initial guess.
22  *  *
23  *  iter    (input/output) int  *  iter    (input/output) int
# Line 29  Line 29 
29  *          norm( b - A*x )  *          norm( b - A*x )
30  *          On output, the final value of this measure.  *          On output, the final value of this measure.
31  *  *
32  *  length_of_recursion (input) gives the number of residual to be kept in orthogonalization process  *  Length_of_recursion (input) gives the number of residual to be kept in orthogonalization process
33  *  *
34  *  restart (input) If restart>0, iteration is resterted a after restart steps.  *  restart (input) If restart>0, iteration is resterted a after restart steps.
35  *  *
36  *  INFO    (output) int  *  INFO    (output) int
37  *  *
38  *         =SOLVER_NO_ERROR: Successful exit. num_iterated approximate solution returned.  *         =SOLVER_NO_ERROR: Successful exit. num_iterated approximate solution returned.
39  *         =SOLVER_MAXNUM_ITEr_REACHED  *         =SOLVER_MAXNUM_ITER_REACHED
40  *         =SOLVER_INPUT_ERROR Illegal parameter:  *         =SOLVER_INPUT_ERROR Illegal parameter:
41  *         =SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller  *         =SOLVER_BREAKDOWN: bad luck!
42  *         =SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller  *         =SOLVER_MEMORY_ERROR : no memory available
43  *  *
44  *  ==============================================================  *  ==============================================================
45  */  */
# Line 56  err_t Paso_Solver_GMRES( Line 56  err_t Paso_Solver_GMRES(
56      double * r,      double * r,
57      double * x,      double * x,
58      dim_t *iter,      dim_t *iter,
59      double * tolerance,dim_t length_of_recursion,dim_t restart) {      double * tolerance,dim_t Length_of_recursion,dim_t restart) {
60    
61    /* Local variables */    /* Local variables */
62    
63    double *x_PRES,*r_PRES,*P,*AP,*x_PRES_MEM[MAX(length_of_recursion,0)],*r_PRES_MEM[MAX(length_of_recursion,0)];    double *AP,*X_PRES[MAX(Length_of_recursion,0)+1],*R_PRES[MAX(Length_of_recursion,0)+1],*P_PRES[MAX(Length_of_recursion,0)+1];
64    double r_PRES_MEMdotAP[MAX(length_of_recursion,0)],L2_r_PRES_MEM[MAX(length_of_recursion,0)],BREAKF_MEM[MAX(length_of_recursion,0)];    double P_PRES_dot_AP[MAX(Length_of_recursion,0)],R_PRES_dot_P_PRES[MAX(Length_of_recursion,0)+1],BREAKF[MAX(Length_of_recursion,0)+1],ALPHA[MAX(Length_of_recursion,0)];
65    double r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1,r_PRES_MEMdotAP2,r_PRES_MEMdotAP3,r_PRES_MEMdotAP4,L2_r_PRES,r_PRES_MEMdotAP5;    double P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6,R_PRES_dot_P,abs_RP,breakf0;
66    double *tmp,tol,BREAKF,factor,GAMMA,SC1,SC2,norm_of_residual,diff,L2_R,norm_of_residual_global;    double tol,Factor,sum_BREAKF,gamma,SC1,SC2,norm_of_residual,diff,L2_R,Norm_of_residual_global;
67    dim_t maxit,num_iter_global,num_iter_restart,num_iter;    double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES;
68    dim_t i,z,OLDEST,ORDER;    dim_t maxit,Num_iter_global,num_iter_restart,num_iter;
69      dim_t i,z,order;
70    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;
71    err_t status=SOLVER_NO_ERROR;    err_t Status=SOLVER_NO_ERROR;
72    
73    /* adapt original routine parameters */    /* adapt original routine parameters */
74    
75    dim_t n_col=A->num_cols * A-> col_block_size;    dim_t n=A->num_cols * A-> col_block_size;
76    dim_t n_row=A->num_rows * A-> row_block_size;    dim_t Length_of_mem=MAX(Length_of_recursion,0)+1;
77    
78    /*     Test the input parameters. */    /*     Test the input parameters. */
79    if (restart>0) restart=MAX(length_of_recursion,restart);    if (restart>0) restart=MAX(Length_of_recursion,restart);
80    if (n_col < 0 || n_row<0 || length_of_recursion<=0) {    if (n < 0 || Length_of_recursion<=0) {
81      return SOLVER_INPUT_ERROR;      return SOLVER_INPUT_ERROR;
82    }    }
83    
84    /*     allocate memory: */    /*     allocate memory: */
85    
86    x_PRES=TMPMEMALLOC(n_col,double);    AP=TMPMEMALLOC(n,double);
87    r_PRES=TMPMEMALLOC(n_row,double);    if (AP==NULL) Status=SOLVER_MEMORY_ERROR;
88    P=TMPMEMALLOC(n_col,double);    for (i=0;i<Length_of_mem;i++) {
89    AP=TMPMEMALLOC(n_row,double);      X_PRES[i]=TMPMEMALLOC(n,double);
90    if (x_PRES==NULL || r_PRES==NULL || P==NULL || AP==NULL) status=SOLVER_MEMORY_ERROR;      R_PRES[i]=TMPMEMALLOC(n,double);
91    for (i=0;i<length_of_recursion;i++) {      P_PRES[i]=TMPMEMALLOC(n,double);
92      x_PRES_MEM[i]=TMPMEMALLOC(n_col,double);      if (X_PRES[i]==NULL || R_PRES[i]==NULL ||  P_PRES[i]==NULL) Status=SOLVER_MEMORY_ERROR;
     r_PRES_MEM[i]=TMPMEMALLOC(n_row,double);  
     if (x_PRES_MEM[i]==NULL || r_PRES_MEM[i]==NULL) status=SOLVER_MEMORY_ERROR;  
93    }    }
94    if ( status ==SOLVER_NO_ERROR ) {    if ( Status ==SOLVER_NO_ERROR ) {
95    
96      /* now PRES starts : */      /* now PRES starts : */
97      maxit=*iter;      maxit=*iter;
98      tol=*tolerance;      tol=*tolerance;
     norm_of_residual=tol;  
99    
100      #pragma omp parallel firstprivate(maxit,tol) \      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \
101         private(num_iter,i,num_iter_restart,ORDER,OLDEST,BREAKF,factor,GAMMA,restartFlag,convergeFlag,maxIterFlag,breakFlag)         private(num_iter,i,num_iter_restart,order,sum_BREAKF,gamma,restartFlag,norm_of_residual,abs_RP,breakf0,\
102                   save_XPRES,save_P_PRES,save_R_PRES,save_R_PRES_dot_P_PRES)
103      {      {
104        /* initialization */        /* initialization */
105    
106        restartFlag=TRUE;        restartFlag=TRUE;
107        num_iter=0;        num_iter=0;
108        #pragma omp for private(z) schedule(static) nowait        #pragma omp for private(z) schedule(static) nowait
109        for (z=0; z < n_row; ++z) AP[z]=0;        for (z=0; z < n; ++z) AP[z]=0;
110        #pragma omp for private(z) schedule(static) nowait        for(i=0;i<Length_of_mem;++i) {
       for (z=0; z < n_col; ++z) P[z]=0;  
       for(i=0;i<length_of_recursion;++i) {  
         #pragma omp for private(z) schedule(static) nowait  
         for (z=0; z < n_row; ++z) r_PRES_MEM[i][z]=0;  
111          #pragma omp for private(z) schedule(static) nowait          #pragma omp for private(z) schedule(static) nowait
112          for (z=0; z < n_col; ++z) x_PRES_MEM[i][z]=0;          for (z=0; z < n; ++z) {
113                       P_PRES[i][z]=0;
114                       R_PRES[i][z]=0;
115                       X_PRES[i][z]=0;
116            }
117        }        }
118    
119        next:        while (! (convergeFlag || maxIterFlag || breakFlag))  {
120           #pragma omp barrier           #pragma omp barrier
121           if (restartFlag) {           if (restartFlag) {
122                 #pragma omp master
123                 BREAKF[0]=ONE;
124               #pragma omp for private(z) schedule(static) nowait               #pragma omp for private(z) schedule(static) nowait
125               for (z=0; z < n_row; ++z) r_PRES[z]=r[z];               for (z=0; z < n; ++z) {
126               #pragma omp for private(z) schedule(static) nowait                  R_PRES[0][z]=r[z];
127               for (z=0; z < n_col; ++z) x_PRES[z]=x[z];                  X_PRES[0][z]=x[z];
128                 }
129               num_iter_restart=0;               num_iter_restart=0;
              OLDEST=length_of_recursion-1;  
              BREAKF=ONE;  
130               restartFlag=FALSE;               restartFlag=FALSE;
131           }           }
132           ++num_iter;           ++num_iter;
133           ++num_iter_restart;           ++num_iter_restart;
134           /* ORDER is the dimension of the space on which the residual is minimized: */           /* order is the dimension of the space on which the residual is minimized: */
135           ORDER=MIN(num_iter_restart-1,length_of_recursion);           order=MIN(num_iter_restart,Length_of_recursion);
          /* OLDEST points to the oldest r_PRES and x_PRES in memory :*/  
          OLDEST=(OLDEST==length_of_recursion-1) ? 0 : OLDEST+1;  
   
136           /***                                                                           /***                                                                
137           *** calculate new search direction P from r_PRES           *** calculate new search direction P from R_PRES
138           ***/           ***/
139           Paso_Solver_solvePreconditioner(A,&P[0], &r_PRES[0]);           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);
140           /***                                                                           /***                                                                
141           *** apply A to P to get AP           *** apply A to P to get AP
142           ***/           ***/
143       Paso_SystemMatrix_MatrixVector(ONE, A, &P[0],ZERO, &AP[0]);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);
144           /***                                                                           /***                                                                
145           ***** calculation of the norm of R and the scalar products of                 ***** calculation of the norm of R and the scalar products of      
146           ***   the residuals and A*P:                                                   ***   the residuals and A*P:                                        
147           ***/           ***/
148           if (ORDER==0) {           if (order==0) {
149              #pragma omp master              #pragma omp master
150              {              R_PRES_dot_P=ZERO;
                L2_r_PRES=ZERO;  
                r_PRESdotAP=ZERO;  
             }  
151              #pragma omp barrier              #pragma omp barrier
152              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P) schedule(static)
153              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z)
154                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
155                  r_PRESdotAP+=r_PRES[z]*AP[z];              #pragma omp master
156              }              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
157           } else if (ORDER==1) {           } else if (order==1) {
158              #pragma omp master              #pragma omp master
159              {              {
160                L2_r_PRES=ZERO;                 R_PRES_dot_P=ZERO;
161                r_PRESdotAP=ZERO;                 P_PRES_dot_AP0=ZERO;
               r_PRES_MEMdotAP0=ZERO;  
162              }              }
163              #pragma omp barrier              #pragma omp barrier
164              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0) schedule(static)
165              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
166                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
167                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
                 r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];  
168              }              }
169              #pragma omp master              #pragma omp master
170              {              {
171                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
172                  R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
173              }              }
174               } else if (order==2) {
          } else if (ORDER==2) {  
175              #pragma omp master              #pragma omp master
176              {              {
177                 L2_r_PRES=ZERO;                R_PRES_dot_P=ZERO;
178                 r_PRESdotAP=ZERO;                P_PRES_dot_AP0=ZERO;
179                 r_PRES_MEMdotAP0=ZERO;                P_PRES_dot_AP1=ZERO;
                r_PRES_MEMdotAP1=ZERO;  
180              }              }
181              #pragma omp barrier              #pragma omp barrier
182              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1) schedule(static)
183              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
184                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
185                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
186                  r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
                 r_PRES_MEMdotAP1+=r_PRES_MEM[1][z]*AP[z];  
187              }              }
188              #pragma omp master              #pragma omp master
189              {              {
190                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
191                r_PRES_MEMdotAP[1]=r_PRES_MEMdotAP1;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
192                  R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
193              }              }
194             } else if (order==3) {
          } else if (ORDER==3) {  
195              #pragma omp master              #pragma omp master
196              {              {
197                L2_r_PRES=ZERO;                 R_PRES_dot_P=ZERO;
198                r_PRESdotAP=ZERO;                 P_PRES_dot_AP0=ZERO;
199                r_PRES_MEMdotAP0=ZERO;                 P_PRES_dot_AP1=ZERO;
200                r_PRES_MEMdotAP1=ZERO;                 P_PRES_dot_AP2=ZERO;
               r_PRES_MEMdotAP2=ZERO;  
201              }              }
202              #pragma omp barrier              #pragma omp barrier
203              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1,r_PRES_MEMdotAP2) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2) schedule(static)
204              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
205                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
206                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
207                  r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
208                  r_PRES_MEMdotAP1+=r_PRES_MEM[1][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
                 r_PRES_MEMdotAP2+=r_PRES_MEM[2][z]*AP[z];  
209              }              }
210              #pragma omp master              #pragma omp master
211              {              {
212                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
213                r_PRES_MEMdotAP[1]=r_PRES_MEMdotAP1;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
214                r_PRES_MEMdotAP[2]=r_PRES_MEMdotAP2;                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
215                  R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
216              }              }
217           } else if (ORDER==4) {           } else if (order==4) {
218              #pragma omp master              #pragma omp master
219              {              {
220                L2_r_PRES=ZERO;                R_PRES_dot_P=ZERO;
221                r_PRESdotAP=ZERO;                P_PRES_dot_AP0=ZERO;
222                r_PRES_MEMdotAP0=ZERO;                P_PRES_dot_AP1=ZERO;
223                r_PRES_MEMdotAP1=ZERO;                P_PRES_dot_AP2=ZERO;
224                r_PRES_MEMdotAP2=ZERO;                P_PRES_dot_AP3=ZERO;
               r_PRES_MEMdotAP3=ZERO;  
225              }              }
226              #pragma omp barrier              #pragma omp barrier
227              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1,r_PRES_MEMdotAP2,r_PRES_MEMdotAP3) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3) schedule(static)
228              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
229                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
230                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
231                  r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
232                  r_PRES_MEMdotAP1+=r_PRES_MEM[1][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
233                  r_PRES_MEMdotAP2+=r_PRES_MEM[2][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
                 r_PRES_MEMdotAP3+=r_PRES_MEM[3][z]*AP[z];  
234              }              }
235              #pragma omp master              #pragma omp master
236              {              {
237                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
238                r_PRES_MEMdotAP[1]=r_PRES_MEMdotAP1;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
239                r_PRES_MEMdotAP[2]=r_PRES_MEMdotAP2;                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
240                r_PRES_MEMdotAP[3]=r_PRES_MEMdotAP3;                P_PRES_dot_AP[3]=P_PRES_dot_AP3;
241              }                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
242           } else if (ORDER==5) {              }
243              #pragma omp master           } else if (order==5) {
244              {              #pragma omp master
245                L2_r_PRES=ZERO;              {
246                r_PRESdotAP=ZERO;                R_PRES_dot_P=ZERO;
247                r_PRES_MEMdotAP0=ZERO;                P_PRES_dot_AP0=ZERO;
248                r_PRES_MEMdotAP1=ZERO;                P_PRES_dot_AP1=ZERO;
249                r_PRES_MEMdotAP2=ZERO;                P_PRES_dot_AP2=ZERO;
250                r_PRES_MEMdotAP3=ZERO;                P_PRES_dot_AP3=ZERO;
251                r_PRES_MEMdotAP4=ZERO;                P_PRES_dot_AP4=ZERO;
252              }              }
253              #pragma omp barrier              #pragma omp barrier
254              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1,r_PRES_MEMdotAP2,r_PRES_MEMdotAP3,r_PRES_MEMdotAP4) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4) schedule(static)
255              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
256                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
257                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
258                  r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
259                  r_PRES_MEMdotAP1+=r_PRES_MEM[1][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
260                  r_PRES_MEMdotAP2+=r_PRES_MEM[2][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
261                  r_PRES_MEMdotAP3+=r_PRES_MEM[3][z]*AP[z];                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
262                  r_PRES_MEMdotAP4+=r_PRES_MEM[4][z]*AP[z];              }
              }  
263              #pragma omp master              #pragma omp master
264              {              {
265                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
266                r_PRES_MEMdotAP[1]=r_PRES_MEMdotAP1;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
267                r_PRES_MEMdotAP[2]=r_PRES_MEMdotAP2;                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
268                r_PRES_MEMdotAP[3]=r_PRES_MEMdotAP3;                P_PRES_dot_AP[3]=P_PRES_dot_AP3;
269                r_PRES_MEMdotAP[4]=r_PRES_MEMdotAP4;                P_PRES_dot_AP[4]=P_PRES_dot_AP4;
270              }                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
271           } else if (ORDER==6) {              }
272              #pragma omp master           } else if (order==6) {
273              {              #pragma omp master
274                L2_r_PRES=ZERO;              {
275                r_PRESdotAP=ZERO;                R_PRES_dot_P=ZERO;
276                r_PRES_MEMdotAP0=ZERO;                P_PRES_dot_AP0=ZERO;
277                r_PRES_MEMdotAP1=ZERO;                P_PRES_dot_AP1=ZERO;
278                r_PRES_MEMdotAP2=ZERO;                P_PRES_dot_AP2=ZERO;
279                r_PRES_MEMdotAP3=ZERO;                P_PRES_dot_AP3=ZERO;
280                r_PRES_MEMdotAP4=ZERO;                P_PRES_dot_AP4=ZERO;
281                r_PRES_MEMdotAP5=ZERO;                P_PRES_dot_AP5=ZERO;
282              }              }
283              #pragma omp barrier              #pragma omp barrier
284              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP,r_PRES_MEMdotAP0,r_PRES_MEMdotAP1,r_PRES_MEMdotAP2,r_PRES_MEMdotAP3,r_PRES_MEMdotAP4,r_PRES_MEMdotAP5) schedule(static)              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5) schedule(static)
285              for (z=0;z<n_row;++z) {              for (z=0;z<n;++z) {
286                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
287                  r_PRESdotAP+=r_PRES[z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
288                  r_PRES_MEMdotAP0+=r_PRES_MEM[0][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
289                  r_PRES_MEMdotAP1+=r_PRES_MEM[1][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
290                  r_PRES_MEMdotAP2+=r_PRES_MEM[2][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
291                  r_PRES_MEMdotAP3+=r_PRES_MEM[3][z]*AP[z];                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
292                  r_PRES_MEMdotAP4+=r_PRES_MEM[4][z]*AP[z];                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
                 r_PRES_MEMdotAP5+=r_PRES_MEM[5][z]*AP[z];  
293               }               }
294              #pragma omp master              #pragma omp master
295              {              {
296                r_PRES_MEMdotAP[0]=r_PRES_MEMdotAP0;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
297                r_PRES_MEMdotAP[1]=r_PRES_MEMdotAP1;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
298                r_PRES_MEMdotAP[2]=r_PRES_MEMdotAP2;                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
299                r_PRES_MEMdotAP[3]=r_PRES_MEMdotAP3;                P_PRES_dot_AP[3]=P_PRES_dot_AP3;
300                r_PRES_MEMdotAP[4]=r_PRES_MEMdotAP4;                P_PRES_dot_AP[4]=P_PRES_dot_AP4;
301                r_PRES_MEMdotAP[5]=r_PRES_MEMdotAP5;                P_PRES_dot_AP[5]=P_PRES_dot_AP5;
302              }                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
303           } else {              }
304             } else if (order>6) {
305                #pragma omp master
306                {
307                  R_PRES_dot_P=ZERO;
308                  P_PRES_dot_AP0=ZERO;
309                  P_PRES_dot_AP1=ZERO;
310                  P_PRES_dot_AP2=ZERO;
311                  P_PRES_dot_AP3=ZERO;
312                  P_PRES_dot_AP4=ZERO;
313                  P_PRES_dot_AP5=ZERO;
314                  P_PRES_dot_AP6=ZERO;
315                }
316                #pragma omp barrier
317                #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6) schedule(static)
318                for (z=0;z<n;++z) {
319                    R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
320                    P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
321                    P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
322                    P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
323                    P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
324                    P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
325                    P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
326                    P_PRES_dot_AP6+=P_PRES[6][z]*AP[z];
327                 }
328              #pragma omp master              #pragma omp master
329              {              {
330                L2_r_PRES=ZERO;                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
331                r_PRESdotAP=ZERO;                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
332                r_PRES_MEMdotAP0=ZERO;                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
333                  P_PRES_dot_AP[3]=P_PRES_dot_AP3;
334                  P_PRES_dot_AP[4]=P_PRES_dot_AP4;
335                  P_PRES_dot_AP[5]=P_PRES_dot_AP5;
336                  P_PRES_dot_AP[6]=P_PRES_dot_AP6;
337                  R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
338    
339                  P_PRES_dot_AP0=ZERO;
340              }              }
341              #pragma omp barrier              for (i=7;i<order;++i) {
342              #pragma omp for private(z) reduction(+:L2_r_PRES,r_PRESdotAP) schedule(static)                  #pragma omp barrier
343              for (z=0;z<n_row;++z) {                  #pragma omp for private(z) reduction(+:P_PRES_dot_AP0) schedule(static)
344                  L2_r_PRES+=r_PRES[z]*r_PRES[z];                  for (z=0;z<n;++z) P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];
                 r_PRESdotAP+=r_PRES[z]*AP[z];  
             }  
             for (i=0;i<ORDER;++i) {  
                 #pragma omp for private(z) reduction(+:r_PRES_MEMdotAP0) schedule(static)  
                 for (z=0;z<n_row;++z) r_PRES_MEMdotAP0+=r_PRES_MEM[i][z]*AP[z];  
345                  #pragma omp master                  #pragma omp master
346                  {                  {
347                    r_PRES_MEMdotAP[i]=r_PRES_MEMdotAP0;                    P_PRES_dot_AP[i]=P_PRES_dot_AP0;
348                    r_PRES_MEMdotAP0=ZERO;                    P_PRES_dot_AP0=ZERO;
349                  }                  }
                 #pragma omp barrier  
350              }              }
351           }           }
352           /* if L2_r_PRES=0 there is a fatal breakdown */           /* this fixes a problem with the intel compiler */
353           if (L2_r_PRES<=ZERO) {           #pragma omp master
354              breakFlag=TRUE;           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];
355           } else {           #pragma omp barrier
356              /***                                                                           /***   if sum_BREAKF is equal to zero a breakdown occurs.
357              ***** calculation of the weights for the update of X:                          ***   iteration procedure can be continued but R_PRES is not the
358              */                                                                          ***   residual of X_PRES approximation.
359              #pragma omp master            ***/
360             sum_BREAKF=0.;
361             for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
362             breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));
363             if (!breakFlag) {
364                breakFlag=FALSE;
365                /***
366                ***** X_PRES and R_PRES are moved to memory:
367                ***/
368                #pragma omp master
369              {              {
370                 r_PRESdotAP*= -ONE/L2_r_PRES;                 Factor=0.;
371                 for (i=0;i<ORDER;++i) r_PRES_MEMdotAP[i]*=-ONE/L2_r_PRES_MEM[i];                 for (i=0;i<order;++i) {
372                       ALPHA[i]=-P_PRES_dot_AP[i]/R_PRES_dot_P_PRES[i];
373                       Factor+=BREAKF[i]*ALPHA[i];
374                   }
375    
376                   save_R_PRES_dot_P_PRES=R_PRES_dot_P_PRES[Length_of_mem-1];
377                   save_R_PRES=R_PRES[Length_of_mem-1];
378                   save_XPRES=X_PRES[Length_of_mem-1];
379                   save_P_PRES=P_PRES[Length_of_mem-1];
380                   for (i=Length_of_mem-1;i>0;--i) {
381                       BREAKF[i]=BREAKF[i-1];
382                       R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1];
383                       R_PRES[i]=R_PRES[i-1];
384                       X_PRES[i]=X_PRES[i-1];
385                       P_PRES[i]=P_PRES[i-1];
386                   }
387                   R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES;
388                   R_PRES[0]=save_R_PRES;
389                   X_PRES[0]=save_XPRES;
390                   P_PRES[0]=save_P_PRES;
391    
392                   if (ABS(Factor)<=ZERO) {
393                      Factor=1.;
394                      BREAKF[0]=ZERO;
395                   } else {
396                      Factor=1./Factor;
397                      BREAKF[0]=ONE;
398                   }
399                   for (i=0;i<order;++i) ALPHA[i]*=Factor;
400              }              }
401              /*                                                                              /*                                                                
402              ***** update of solution x_PRES and its residual r_PRES:                          ***** update of solution X_PRES and its residual R_PRES:            
403              ***                                                                            ***                                                              
404              ***   P is used to accumulate X and AP to accumulate R. X and R                  ***   P is used to accumulate X and AP to accumulate R. X and R    
405              ***   are still needed until they are put into the X and R memory                ***   are still needed until they are put into the X and R memory  
406              ***   r_PRES_MEM and x_PRES_MEM                                                  ***   R_PRES and X_PRES                                    
407              ***                                                                              ***                                                                
408              **/              **/
409              #pragma omp barrier              #pragma omp barrier
410              if (ORDER==0) {              breakf0=BREAKF[0];
411                #pragma omp for private(z) schedule(static)              if (order==0) {
               for (z=0;z<n_row;++z)  
                   AP[z]+=r_PRESdotAP*r_PRES[z];  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_col;++z)  
                   P[z]=-P[z]+r_PRESdotAP*x_PRES[z];  
             } else if (ORDER==1) {  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_row;++z)  
                  AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z];  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_col;++z)  
                 P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z];  
             } else if (ORDER==2) {  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_row;++z)  
                  AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z]  
                                              +r_PRES_MEMdotAP[1]*r_PRES_MEM[1][z];  
412                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
413                for (z=0;z<n_col;++z)                for (z=0;z<n;++z) {
414                   P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z]                    R_PRES[0][z]= Factor*       AP[z];
415                                              +r_PRES_MEMdotAP[1]*x_PRES_MEM[1][z];                    X_PRES[0][z]=-Factor*P_PRES[1][z];
416              } else if (ORDER==3) {                }
417                #pragma omp for private(z) schedule(static)              } else if (order==1) {
               for (z=0;z<n_row;++z)  
                  AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z]  
                                              +r_PRES_MEMdotAP[1]*r_PRES_MEM[1][z]  
                                              +r_PRES_MEMdotAP[2]*r_PRES_MEM[2][z];  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_col;++z)  
                  P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z]  
                                             +r_PRES_MEMdotAP[1]*x_PRES_MEM[1][z]  
                                             +r_PRES_MEMdotAP[2]*x_PRES_MEM[2][z];  
             } else if (ORDER==4) {  
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_row;++z)  
                  AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z]  
                                              +r_PRES_MEMdotAP[1]*r_PRES_MEM[1][z]  
                                              +r_PRES_MEMdotAP[2]*r_PRES_MEM[2][z]  
                                              +r_PRES_MEMdotAP[3]*r_PRES_MEM[3][z];  
418                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
419                for (z=0;z<n_col;++z)                for (z=0;z<n;++z) {
420                   P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z]                    R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z];
421                                              +r_PRES_MEMdotAP[1]*x_PRES_MEM[1][z]                    X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];
422                                              +r_PRES_MEMdotAP[2]*x_PRES_MEM[2][z]                }
423                                              +r_PRES_MEMdotAP[3]*x_PRES_MEM[3][z];              } else if (order==2) {
             } else if (ORDER==5) {  
424                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
425                for (z=0;z<n_row;++z)                for (z=0;z<n;++z) {
426                   AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z]                   R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z]
427                                               +r_PRES_MEMdotAP[1]*r_PRES_MEM[1][z]                                                    +ALPHA[1]*R_PRES[2][z];
428                                               +r_PRES_MEMdotAP[2]*r_PRES_MEM[2][z]                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
429                                               +r_PRES_MEMdotAP[3]*r_PRES_MEM[3][z]                                                    +ALPHA[1]*X_PRES[2][z];
430                                               +r_PRES_MEMdotAP[4]*r_PRES_MEM[4][z];                }
431                } else if (order==3) {
432                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
433                for (z=0;z<n_col;++z)                for (z=0;z<n;++z) {
434                   P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z]                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
435                                              +r_PRES_MEMdotAP[1]*x_PRES_MEM[1][z]                                             +ALPHA[1]*R_PRES[2][z]
436                                              +r_PRES_MEMdotAP[2]*x_PRES_MEM[2][z]                                             +ALPHA[2]*R_PRES[3][z];
437                                              +r_PRES_MEMdotAP[3]*x_PRES_MEM[3][z]                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
438                                              +r_PRES_MEMdotAP[4]*x_PRES_MEM[4][z];                                                    +ALPHA[1]*X_PRES[2][z]
439              } else if (ORDER==6) {                                                    +ALPHA[2]*X_PRES[3][z];
440                  }
441                } else if (order==4) {
442                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
443                for (z=0;z<n_row;++z)                for (z=0;z<n;++z) {
444                   AP[z]+=r_PRESdotAP*r_PRES[z]+r_PRES_MEMdotAP[0]*r_PRES_MEM[0][z]                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
445                                               +r_PRES_MEMdotAP[1]*r_PRES_MEM[1][z]                                             +ALPHA[1]*R_PRES[2][z]
446                                               +r_PRES_MEMdotAP[2]*r_PRES_MEM[2][z]                                             +ALPHA[2]*R_PRES[3][z]
447                                               +r_PRES_MEMdotAP[3]*r_PRES_MEM[3][z]                                             +ALPHA[3]*R_PRES[4][z];
448                                               +r_PRES_MEMdotAP[4]*r_PRES_MEM[4][z]                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
449                                               +r_PRES_MEMdotAP[5]*r_PRES_MEM[5][z];                                                    +ALPHA[1]*X_PRES[2][z]
450                                                      +ALPHA[2]*X_PRES[3][z]
451                                                      +ALPHA[3]*X_PRES[4][z];
452                  }
453                } else if (order==5) {
454                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
455                for (z=0;z<n_col;++z)                for (z=0;z<n;++z) {
456                   P[z]=-P[z]+r_PRESdotAP*x_PRES[z]+r_PRES_MEMdotAP[0]*x_PRES_MEM[0][z]                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
457                                              +r_PRES_MEMdotAP[1]*x_PRES_MEM[1][z]                                            +ALPHA[1]*R_PRES[2][z]
458                                              +r_PRES_MEMdotAP[2]*x_PRES_MEM[2][z]                                            +ALPHA[2]*R_PRES[3][z]
459                                              +r_PRES_MEMdotAP[3]*x_PRES_MEM[3][z]                                            +ALPHA[3]*R_PRES[4][z]
460                                              +r_PRES_MEMdotAP[4]*x_PRES_MEM[4][z]                                            +ALPHA[4]*R_PRES[5][z];
461                                              +r_PRES_MEMdotAP[5]*x_PRES_MEM[5][z];                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
462              } else {                                                    +ALPHA[1]*X_PRES[2][z]
463                                                      +ALPHA[2]*X_PRES[3][z]
464                                                      +ALPHA[3]*X_PRES[4][z]
465                                                      +ALPHA[4]*X_PRES[5][z];
466                  }
467                } else if (order==6) {
468                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
469                for (z=0;z<n_row;++z)                for (z=0;z<n;++z) {
470                  AP[z]+=r_PRESdotAP*r_PRES[z];                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
471                                              +ALPHA[1]*R_PRES[2][z]
472                                              +ALPHA[2]*R_PRES[3][z]
473                                              +ALPHA[3]*R_PRES[4][z]
474                                              +ALPHA[4]*R_PRES[5][z]
475                                              +ALPHA[5]*R_PRES[6][z];
476                     X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
477                                                      +ALPHA[1]*X_PRES[2][z]
478                                                      +ALPHA[2]*X_PRES[3][z]
479                                                      +ALPHA[3]*X_PRES[4][z]
480                                                      +ALPHA[4]*X_PRES[5][z]
481                                                      +ALPHA[5]*X_PRES[6][z];
482                  }
483                } else if (order>6) {
484                #pragma omp for private(z) schedule(static)                #pragma omp for private(z) schedule(static)
485                for (z=0;z<n_col;++z)                for (z=0;z<n;++z) {
486                  P[z]=-P[z]+r_PRESdotAP*x_PRES[z];                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
487                                              +ALPHA[1]*R_PRES[2][z]
488                for (i=0;i<ORDER;++i) {                                            +ALPHA[2]*R_PRES[3][z]
489                  #pragma omp for private(z) schedule(static)                                            +ALPHA[3]*R_PRES[4][z]
490                  for (z=0;z<n_row;++z)                                            +ALPHA[4]*R_PRES[5][z]
491                      AP[z]+=r_PRES_MEMdotAP[i]*r_PRES_MEM[i][z];                                            +ALPHA[5]*R_PRES[6][z]
492                                              +ALPHA[6]*R_PRES[7][z];
493                     X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
494                                                      +ALPHA[1]*X_PRES[2][z]
495                                                      +ALPHA[2]*X_PRES[3][z]
496                                                      +ALPHA[3]*X_PRES[4][z]
497                                                      +ALPHA[4]*X_PRES[5][z]
498                                                      +ALPHA[5]*X_PRES[6][z]
499                                                      +ALPHA[6]*X_PRES[7][z];
500                  }
501                  for (i=7;i<order;++i) {
502                  #pragma omp for private(z) schedule(static)                  #pragma omp for private(z) schedule(static)
503                  for (z=0;z<n_col;++z)                  for (z=0;z<n;++z) {
504                      P[z]+=r_PRES_MEMdotAP[i]*x_PRES_MEM[i][z];                      R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];
505                        X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];
506                    }
507                }                }
508              }              }
509              /***  factor scales AP and P to make it a residual and              if (breakf0>0.) {
             ***   as solution approximation  
             ***  
             ***   if factor is equal to zero a breakdown occurs. the  
             ***   iteration procedure can be continued but r_PRES is not the  
             ***   residual of x_PRES approximation.  
             ***/  
             factor=BREAKF*r_PRESdotAP;  
             for (i=0;i<ORDER;++i) factor +=BREAKF_MEM[i]*r_PRES_MEMdotAP[i];  
             /***  
             ***** x_PRES and r_PRES are moved to memory:  
             ***/  
             #pragma omp master  
             {  
                L2_r_PRES_MEM[OLDEST]=L2_r_PRES;  
                BREAKF_MEM[OLDEST]=BREAKF;  
                tmp=x_PRES; x_PRES=x_PRES_MEM[OLDEST];x_PRES_MEM[OLDEST]=tmp;  
                tmp=r_PRES; r_PRES=r_PRES_MEM[OLDEST];r_PRES_MEM[OLDEST]=tmp;  
             }  
             if (ABS(factor)<=ZERO) {  
                 /* in case of a break down: */  
                 BREAKF=ZERO;  
                 #pragma omp for private(z) schedule(static)  
                 for (z=0;z<n_row;++z) r_PRES[z]=AP[z];  
                 #pragma omp for private(z) schedule(static)  
                 for (z=0;z<n_col;++z) x_PRES[z]=P[z];  
                 /* is there any progress */  
                 breakFlag=TRUE;  
                 #pragma omp barrier  
                 for (i=0;i<ORDER;++i) if (BREAKF_MEM[i]>ZERO) breakFlag=FALSE;  
                 convergeFlag = FALSE;  
                 maxIterFlag = FALSE;  
             } else {  
               BREAKF=ONE;  
               breakFlag=FALSE;  
510                /***                /***
511                *** rescale P an AP and apply smooting:                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:
               ***  
               ***** calculate GAMMA from min_(GAMMA){|R+GAMMA*(r_PRES-R)|_2}:  
512                ***/                ***/
513                #pragma omp master                #pragma omp master
514                {                {
# Line 502  err_t Paso_Solver_GMRES( Line 516  err_t Paso_Solver_GMRES(
516                   SC2=ZERO;                   SC2=ZERO;
517                   L2_R=ZERO;                   L2_R=ZERO;
518                }                }
               factor=ONE/factor;  
519                #pragma omp barrier                #pragma omp barrier
520                #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static)                #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static)
521                for (z=0;z<n_row;++z) {                for (z=0;z<n;++z) {
522                  r_PRES[z]=factor*AP[z];                  diff=R_PRES[0][z]-r[z];
                 diff=r_PRES[z]-r[z];  
523                  SC1+=diff*diff;                  SC1+=diff*diff;
524                  SC2+=diff*r[z];                  SC2+=diff*r[z];
525                }                }
526                GAMMA=(SC1<=ZERO) ? ZERO : -SC2/SC1;                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;
               #pragma omp for private(z) schedule(static)  
               for (z=0;z<n_col;++z) {  
                  x_PRES[z]=factor* P[z];  
                  x[z]+=GAMMA*(x_PRES[z]-x[z]);  
               }  
527                #pragma omp for private(z) reduction(+:L2_R) schedule(static)                #pragma omp for private(z) reduction(+:L2_R) schedule(static)
528                for (z=0;z<n_row;++z) {                for (z=0;z<n;++z) {
529                   r[z]+=GAMMA*(r_PRES[z]-r[z]);                   x[z]+=gamma*(X_PRES[0][z]-x[z]);
530                     r[z]+=gamma*(R_PRES[0][z]-r[z]);
531                   L2_R+=r[z]*r[z];                   L2_R+=r[z]*r[z];
532                }                }
533                norm_of_residual=sqrt(L2_R);                norm_of_residual=sqrt(L2_R);
534                convergeFlag = (norm_of_residual <= tol);                convergeFlag = (norm_of_residual <= tol);
535                maxIterFlag = (num_iter >= maxit);                if (restart>0) restartFlag=(num_iter_restart >= restart);
536                } else {
537                  convergeFlag=FALSE;
538                  restartFlag=FALSE;
539              }              }
540           }              maxIterFlag = (num_iter >= maxit);
541           if (restart>0) restartFlag=(num_iter_restart >= restart);           }
542           if (! (convergeFlag || maxIterFlag || breakFlag))  goto next;          }
543           /* end of iteration */          /* end of iteration */
544           #pragma omp master          #pragma omp master
545           {          {
546             norm_of_residual_global=norm_of_residual;             Norm_of_residual_global=norm_of_residual;
547         num_iter_global=num_iter;         Num_iter_global=num_iter;
548         if (maxIterFlag) {         if (maxIterFlag) {
549             status = SOLVER_MAXITER_REACHED;             Status = SOLVER_MAXITER_REACHED;
550         } else if (breakFlag) {         } else if (breakFlag) {
551             status = SOLVER_BREAKDOWN;             Status = SOLVER_BREAKDOWN;
552         }        }
553          }          }
554      }  /* end of parallel region */      }  /* end of parallel region */
     TMPMEMFREE(x_PRES);  
     TMPMEMFREE(r_PRES);  
     TMPMEMFREE(P);  
555      TMPMEMFREE(AP);      TMPMEMFREE(AP);
556      for (i=0;i<length_of_recursion;i++) {      for (i=0;i<Length_of_recursion;i++) {
557         TMPMEMFREE(x_PRES_MEM[i]);         TMPMEMFREE(X_PRES[i]);
558         TMPMEMFREE(r_PRES_MEM[i]);         TMPMEMFREE(R_PRES[i]);
559           TMPMEMFREE(P_PRES[i]);
560      }      }
561      *iter=num_iter_global;      *iter=Num_iter_global;
562      *tolerance=norm_of_residual_global;      *tolerance=Norm_of_residual_global;
563    }    }
564    return status;    return Status;
565  }  }
566    
   
 /*  
  * $Log$  
  * Revision 1.2  2005/09/15 03:44:40  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.1.2.1  2005/09/05 06:29:49  gross  
  * These files have been extracted from finley to define a stand alone libray for iterative  
  * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  
  * has not been tested yet.  
  *  
  *  
  */  
   
     

Legend:
Removed from v.150  
changed lines
  Added in v.415

  ViewVC Help
Powered by ViewVC 1.1.26