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

revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1556 by gross, Mon May 12 00:54:58 2008 UTC
# Line 121  err_t Paso_Solver_GMRES( Line 121  err_t Paso_Solver_GMRES(
121      maxit=*iter;      maxit=*iter;
122      tol=*tolerance;      tol=*tolerance;
123    
124      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \       /* initialization */
        private(num_iter,i,num_iter_restart,order,sum_BREAKF,gamma,restartFlag,norm_of_residual,abs_RP,breakf0,\  
                save_XPRES,save_P_PRES,save_R_PRES,save_R_PRES_dot_P_PRES)  
     {  
       /* initialization */  
125    
126        restartFlag=TRUE;        restartFlag=TRUE;
127        num_iter=0;        num_iter=0;
128        #pragma omp for private(z) schedule(static) nowait        #pragma omp parallel private(i)
129        for (z=0; z < n; ++z) AP[z]=0;        {
130        for(i=0;i<Length_of_mem;++i) {           #pragma omp for private(z) schedule(static) nowait
131          #pragma omp for private(z) schedule(static) nowait           for (z=0; z < n; ++z) AP[z]=0;
132          for (z=0; z < n; ++z) {           for(i=0;i<Length_of_mem;++i) {
133                     P_PRES[i][z]=0;             #pragma omp for private(z) schedule(static) nowait
134                     R_PRES[i][z]=0;             for (z=0; z < n; ++z) {
135                     X_PRES[i][z]=0;                        P_PRES[i][z]=0;
136          }                        R_PRES[i][z]=0;
137                          X_PRES[i][z]=0;
138               }
139             }
140        }        }
141    
142        while (! (convergeFlag || maxIterFlag || breakFlag))  {        while (! (convergeFlag || maxIterFlag || breakFlag))  {
          #pragma omp barrier  
143           if (restartFlag) {           if (restartFlag) {
              #pragma omp master  
144               BREAKF[0]=ONE;               BREAKF[0]=ONE;
145               #pragma omp for private(z) schedule(static) nowait               #pragma omp parallel for private(z) schedule(static)
146               for (z=0; z < n; ++z) {               for (z=0; z < n; ++z) {
147                  R_PRES[0][z]=r[z];                  R_PRES[0][z]=r[z];
148                  X_PRES[0][z]=x[z];                  X_PRES[0][z]=x[z];
# Line 160  err_t Paso_Solver_GMRES( Line 157  err_t Paso_Solver_GMRES(
157           /***                                                                           /***                                                                
158           *** calculate new search direction P from R_PRES           *** calculate new search direction P from R_PRES
159           ***/           ***/
          #pragma omp barrier  
160           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);
161           /***                                                                           /***                                                                
162           *** apply A to P to get AP           *** apply A to P to get AP
163           ***/           ***/
          #pragma omp barrier  
164       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);
165           /***                                                                           /***                                                                
166           ***** calculation of the norm of R and the scalar products of                 ***** calculation of the norm of R and the scalar products of      
167           ***   the residuals and A*P:                                                   ***   the residuals and A*P:                                        
168           ***/           ***/
169           if (order==0) {           if (order==0) {
             #pragma omp master  
170              R_PRES_dot_P=ZERO;              R_PRES_dot_P=ZERO;
171              #pragma omp barrier              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P) schedule(static)
172              #pragma omp for private(z) reduction(+:R_PRES_dot_P) schedule(static)              for (z=0;z<n;++z) R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
             for (z=0;z<n;++z)  
                 R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];  
             #pragma omp master  
173              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
174           } else if (order==1) {           } else if (order==1) {
175              #pragma omp master              R_PRES_dot_P=ZERO;
176              {              P_PRES_dot_AP0=ZERO;
177                 R_PRES_dot_P=ZERO;              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0) schedule(static)
                P_PRES_dot_AP0=ZERO;  
             }  
             #pragma omp barrier  
             #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0) schedule(static)  
178              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
179                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
180                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
181              }              }
182              #pragma omp master              P_PRES_dot_AP[0]=P_PRES_dot_AP0;
183              {              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
               P_PRES_dot_AP[0]=P_PRES_dot_AP0;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
             }  
184           } else if (order==2) {           } else if (order==2) {
185              #pragma omp master              R_PRES_dot_P=ZERO;
186              {              P_PRES_dot_AP0=ZERO;
187                R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
188                P_PRES_dot_AP0=ZERO;              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1) schedule(static)
               P_PRES_dot_AP1=ZERO;  
             }  
             #pragma omp barrier  
             #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1) schedule(static)  
189              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
190                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
191                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
192                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
193              }              }
194              #pragma omp master              P_PRES_dot_AP[0]=P_PRES_dot_AP0;
195              {              P_PRES_dot_AP[1]=P_PRES_dot_AP1;
196                P_PRES_dot_AP[0]=P_PRES_dot_AP0;              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
               P_PRES_dot_AP[1]=P_PRES_dot_AP1;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
             }  
197           } else if (order==3) {           } else if (order==3) {
198              #pragma omp master              R_PRES_dot_P=ZERO;
199              {              P_PRES_dot_AP0=ZERO;
200                 R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
201                 P_PRES_dot_AP0=ZERO;              P_PRES_dot_AP2=ZERO;
202                 P_PRES_dot_AP1=ZERO;              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2) schedule(static)
                P_PRES_dot_AP2=ZERO;  
             }  
             #pragma omp barrier  
             #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2) schedule(static)  
203              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
204                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
205                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
206                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
207                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
208              }              }
209              #pragma omp master              P_PRES_dot_AP[0]=P_PRES_dot_AP0;
210              {              P_PRES_dot_AP[1]=P_PRES_dot_AP1;
211                P_PRES_dot_AP[0]=P_PRES_dot_AP0;              P_PRES_dot_AP[2]=P_PRES_dot_AP2;
212                P_PRES_dot_AP[1]=P_PRES_dot_AP1;              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
               P_PRES_dot_AP[2]=P_PRES_dot_AP2;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
             }  
213           } else if (order==4) {           } else if (order==4) {
214              #pragma omp master              R_PRES_dot_P=ZERO;
215              {              P_PRES_dot_AP0=ZERO;
216                R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
217                P_PRES_dot_AP0=ZERO;              P_PRES_dot_AP2=ZERO;
218                P_PRES_dot_AP1=ZERO;              P_PRES_dot_AP3=ZERO;
219                P_PRES_dot_AP2=ZERO;              #pragma omp parallel 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)
               P_PRES_dot_AP3=ZERO;  
             }  
             #pragma omp barrier  
             #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)  
220              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
221                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
222                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
# Line 258  err_t Paso_Solver_GMRES( Line 224  err_t Paso_Solver_GMRES(
224                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
225                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
226              }              }
227              #pragma omp master              P_PRES_dot_AP[0]=P_PRES_dot_AP0;
228              {              P_PRES_dot_AP[1]=P_PRES_dot_AP1;
229                P_PRES_dot_AP[0]=P_PRES_dot_AP0;              P_PRES_dot_AP[2]=P_PRES_dot_AP2;
230                P_PRES_dot_AP[1]=P_PRES_dot_AP1;              P_PRES_dot_AP[3]=P_PRES_dot_AP3;
231                P_PRES_dot_AP[2]=P_PRES_dot_AP2;              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
               P_PRES_dot_AP[3]=P_PRES_dot_AP3;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
             }  
232           } else if (order==5) {           } else if (order==5) {
233              #pragma omp master              R_PRES_dot_P=ZERO;
234              {              P_PRES_dot_AP0=ZERO;
235                R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
236                P_PRES_dot_AP0=ZERO;              P_PRES_dot_AP2=ZERO;
237                P_PRES_dot_AP1=ZERO;              P_PRES_dot_AP3=ZERO;
238                P_PRES_dot_AP2=ZERO;              P_PRES_dot_AP4=ZERO;
239                P_PRES_dot_AP3=ZERO;              #pragma omp parallel 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)
               P_PRES_dot_AP4=ZERO;  
             }  
             #pragma omp barrier  
             #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)  
240              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
241                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
242                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
# Line 286  err_t Paso_Solver_GMRES( Line 245  err_t Paso_Solver_GMRES(
245                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
246                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
247              }              }
248              #pragma omp master              P_PRES_dot_AP[0]=P_PRES_dot_AP0;
249              {              P_PRES_dot_AP[1]=P_PRES_dot_AP1;
250                P_PRES_dot_AP[0]=P_PRES_dot_AP0;              P_PRES_dot_AP[2]=P_PRES_dot_AP2;
251                P_PRES_dot_AP[1]=P_PRES_dot_AP1;              P_PRES_dot_AP[3]=P_PRES_dot_AP3;
252                P_PRES_dot_AP[2]=P_PRES_dot_AP2;              P_PRES_dot_AP[4]=P_PRES_dot_AP4;
253                P_PRES_dot_AP[3]=P_PRES_dot_AP3;              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
               P_PRES_dot_AP[4]=P_PRES_dot_AP4;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
             }  
254           } else if (order==6) {           } else if (order==6) {
255              #pragma omp master              R_PRES_dot_P=ZERO;
256              {              P_PRES_dot_AP0=ZERO;
257                R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
258                P_PRES_dot_AP0=ZERO;              P_PRES_dot_AP2=ZERO;
259                P_PRES_dot_AP1=ZERO;              P_PRES_dot_AP3=ZERO;
260                P_PRES_dot_AP2=ZERO;              P_PRES_dot_AP4=ZERO;
261                P_PRES_dot_AP3=ZERO;              P_PRES_dot_AP5=ZERO;
262                P_PRES_dot_AP4=ZERO;              #pragma omp parallel 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)
               P_PRES_dot_AP5=ZERO;  
             }  
             #pragma omp barrier  
             #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)  
263              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
264                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
265                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
# Line 316  err_t Paso_Solver_GMRES( Line 268  err_t Paso_Solver_GMRES(
268                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
269                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
270                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
              }  
             #pragma omp master  
             {  
               P_PRES_dot_AP[0]=P_PRES_dot_AP0;  
               P_PRES_dot_AP[1]=P_PRES_dot_AP1;  
               P_PRES_dot_AP[2]=P_PRES_dot_AP2;  
               P_PRES_dot_AP[3]=P_PRES_dot_AP3;  
               P_PRES_dot_AP[4]=P_PRES_dot_AP4;  
               P_PRES_dot_AP[5]=P_PRES_dot_AP5;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
271              }              }
272                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
273                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
274                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
275                P_PRES_dot_AP[3]=P_PRES_dot_AP3;
276                P_PRES_dot_AP[4]=P_PRES_dot_AP4;
277                P_PRES_dot_AP[5]=P_PRES_dot_AP5;
278                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
279           } else if (order>6) {           } else if (order>6) {
280              #pragma omp master              R_PRES_dot_P=ZERO;
281              {              P_PRES_dot_AP0=ZERO;
282                R_PRES_dot_P=ZERO;              P_PRES_dot_AP1=ZERO;
283                P_PRES_dot_AP0=ZERO;              P_PRES_dot_AP2=ZERO;
284                P_PRES_dot_AP1=ZERO;              P_PRES_dot_AP3=ZERO;
285                P_PRES_dot_AP2=ZERO;              P_PRES_dot_AP4=ZERO;
286                P_PRES_dot_AP3=ZERO;              P_PRES_dot_AP5=ZERO;
287                P_PRES_dot_AP4=ZERO;              P_PRES_dot_AP6=ZERO;
288                P_PRES_dot_AP5=ZERO;              #pragma omp parallel 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)
               P_PRES_dot_AP6=ZERO;  
             }  
             #pragma omp barrier  
             #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)  
289              for (z=0;z<n;++z) {              for (z=0;z<n;++z) {
290                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
291                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
# Line 350  err_t Paso_Solver_GMRES( Line 295  err_t Paso_Solver_GMRES(
295                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
296                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
297                  P_PRES_dot_AP6+=P_PRES[6][z]*AP[z];                  P_PRES_dot_AP6+=P_PRES[6][z]*AP[z];
              }  
             #pragma omp master  
             {  
               P_PRES_dot_AP[0]=P_PRES_dot_AP0;  
               P_PRES_dot_AP[1]=P_PRES_dot_AP1;  
               P_PRES_dot_AP[2]=P_PRES_dot_AP2;  
               P_PRES_dot_AP[3]=P_PRES_dot_AP3;  
               P_PRES_dot_AP[4]=P_PRES_dot_AP4;  
               P_PRES_dot_AP[5]=P_PRES_dot_AP5;  
               P_PRES_dot_AP[6]=P_PRES_dot_AP6;  
               R_PRES_dot_P_PRES[0]=R_PRES_dot_P;  
   
               P_PRES_dot_AP0=ZERO;  
298              }              }
299                P_PRES_dot_AP[0]=P_PRES_dot_AP0;
300                P_PRES_dot_AP[1]=P_PRES_dot_AP1;
301                P_PRES_dot_AP[2]=P_PRES_dot_AP2;
302                P_PRES_dot_AP[3]=P_PRES_dot_AP3;
303                P_PRES_dot_AP[4]=P_PRES_dot_AP4;
304                P_PRES_dot_AP[5]=P_PRES_dot_AP5;
305                P_PRES_dot_AP[6]=P_PRES_dot_AP6;
306                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;
307    
308                P_PRES_dot_AP0=ZERO;
309              for (i=7;i<order;++i) {              for (i=7;i<order;++i) {
310                  #pragma omp barrier                  #pragma omp parallel for private(z) reduction(+:P_PRES_dot_AP0) schedule(static)
                 #pragma omp for private(z) reduction(+:P_PRES_dot_AP0) schedule(static)  
311                  for (z=0;z<n;++z) P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];                  for (z=0;z<n;++z) P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];
312                  #pragma omp master                  P_PRES_dot_AP[i]=P_PRES_dot_AP0;
313                  {                  P_PRES_dot_AP0=ZERO;
                   P_PRES_dot_AP[i]=P_PRES_dot_AP0;  
                   P_PRES_dot_AP0=ZERO;  
                 }  
314              }              }
315           }           }
316           /* this fixes a problem with the intel compiler */           /* this fixes a problem with the intel compiler */
          #pragma omp master  
317           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];
318           /***   if sum_BREAKF is equal to zero a breakdown occurs.           /***   if sum_BREAKF is equal to zero a breakdown occurs.
319            ***   iteration procedure can be continued but R_PRES is not the            ***   iteration procedure can be continued but R_PRES is not the
320            ***   residual of X_PRES approximation.            ***   residual of X_PRES approximation.
321            ***/            ***/
          #pragma omp barrier  
322           sum_BREAKF=0.;           sum_BREAKF=0.;
323           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
324           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));
# Line 391  err_t Paso_Solver_GMRES( Line 327  err_t Paso_Solver_GMRES(
327              /***              /***
328              ***** X_PRES and R_PRES are moved to memory:              ***** X_PRES and R_PRES are moved to memory:
329              ***/              ***/
330              #pragma omp master              Factor=0.;
331              {              for (i=0;i<order;++i) {
                Factor=0.;  
                for (i=0;i<order;++i) {  
332                     ALPHA[i]=-P_PRES_dot_AP[i]/R_PRES_dot_P_PRES[i];                     ALPHA[i]=-P_PRES_dot_AP[i]/R_PRES_dot_P_PRES[i];
333                     Factor+=BREAKF[i]*ALPHA[i];                     Factor+=BREAKF[i]*ALPHA[i];
334                 }              }
335    
336                 save_R_PRES_dot_P_PRES=R_PRES_dot_P_PRES[Length_of_mem-1];               save_R_PRES_dot_P_PRES=R_PRES_dot_P_PRES[Length_of_mem-1];
337                 save_R_PRES=R_PRES[Length_of_mem-1];               save_R_PRES=R_PRES[Length_of_mem-1];
338                 save_XPRES=X_PRES[Length_of_mem-1];               save_XPRES=X_PRES[Length_of_mem-1];
339                 save_P_PRES=P_PRES[Length_of_mem-1];               save_P_PRES=P_PRES[Length_of_mem-1];
340                 for (i=Length_of_mem-1;i>0;--i) {               for (i=Length_of_mem-1;i>0;--i) {
341                     BREAKF[i]=BREAKF[i-1];                     BREAKF[i]=BREAKF[i-1];
342                     R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1];                     R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1];
343                     R_PRES[i]=R_PRES[i-1];                     R_PRES[i]=R_PRES[i-1];
344                     X_PRES[i]=X_PRES[i-1];                     X_PRES[i]=X_PRES[i-1];
345                     P_PRES[i]=P_PRES[i-1];                     P_PRES[i]=P_PRES[i-1];
346                 }               }
347                 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES;               R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES;
348                 R_PRES[0]=save_R_PRES;               R_PRES[0]=save_R_PRES;
349                 X_PRES[0]=save_XPRES;               X_PRES[0]=save_XPRES;
350                 P_PRES[0]=save_P_PRES;               P_PRES[0]=save_P_PRES;
351    
352                 if (ABS(Factor)<=ZERO) {               if (ABS(Factor)<=ZERO) {
353                    Factor=1.;                    Factor=1.;
354                    BREAKF[0]=ZERO;                    BREAKF[0]=ZERO;
355                 } else {               } else {
356                    Factor=1./Factor;                    Factor=1./Factor;
357                    BREAKF[0]=ONE;                    BREAKF[0]=ONE;
                }  
                for (i=0;i<order;++i) ALPHA[i]*=Factor;  
358              }              }
359                for (i=0;i<order;++i) ALPHA[i]*=Factor;
360              /*                                                                              /*                                                                
361              ***** update of solution X_PRES and its residual R_PRES:                          ***** update of solution X_PRES and its residual R_PRES:            
362              ***                                                                            ***                                                              
# Line 432  err_t Paso_Solver_GMRES( Line 365  err_t Paso_Solver_GMRES(
365              ***   R_PRES and X_PRES                                                  ***   R_PRES and X_PRES                                    
366              ***                                                                              ***                                                                
367              **/              **/
             #pragma omp barrier  
368              breakf0=BREAKF[0];              breakf0=BREAKF[0];
369              if (order==0) {              if (order==0) {
370                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
371                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
372                    R_PRES[0][z]= Factor*       AP[z];                    R_PRES[0][z]= Factor*       AP[z];
373                    X_PRES[0][z]=-Factor*P_PRES[1][z];                    X_PRES[0][z]=-Factor*P_PRES[1][z];
374                }                }
375              } else if (order==1) {              } else if (order==1) {
376                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
377                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
378                    R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z];                    R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z];
379                    X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];                    X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];
380                }                }
381              } else if (order==2) {              } else if (order==2) {
382                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
383                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
384                   R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z]
385                                                    +ALPHA[1]*R_PRES[2][z];                                                    +ALPHA[1]*R_PRES[2][z];
# Line 455  err_t Paso_Solver_GMRES( Line 387  err_t Paso_Solver_GMRES(
387                                                    +ALPHA[1]*X_PRES[2][z];                                                    +ALPHA[1]*X_PRES[2][z];
388                }                }
389              } else if (order==3) {              } else if (order==3) {
390                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
391                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
392                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
393                                             +ALPHA[1]*R_PRES[2][z]                                             +ALPHA[1]*R_PRES[2][z]
# Line 465  err_t Paso_Solver_GMRES( Line 397  err_t Paso_Solver_GMRES(
397                                                    +ALPHA[2]*X_PRES[3][z];                                                    +ALPHA[2]*X_PRES[3][z];
398                }                }
399              } else if (order==4) {              } else if (order==4) {
400                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
401                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
402                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
403                                             +ALPHA[1]*R_PRES[2][z]                                             +ALPHA[1]*R_PRES[2][z]
# Line 477  err_t Paso_Solver_GMRES( Line 409  err_t Paso_Solver_GMRES(
409                                                    +ALPHA[3]*X_PRES[4][z];                                                    +ALPHA[3]*X_PRES[4][z];
410                }                }
411              } else if (order==5) {              } else if (order==5) {
412                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
413                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
414                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
415                                            +ALPHA[1]*R_PRES[2][z]                                            +ALPHA[1]*R_PRES[2][z]
# Line 491  err_t Paso_Solver_GMRES( Line 423  err_t Paso_Solver_GMRES(
423                                                    +ALPHA[4]*X_PRES[5][z];                                                    +ALPHA[4]*X_PRES[5][z];
424                }                }
425              } else if (order==6) {              } else if (order==6) {
426                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
427                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
428                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
429                                            +ALPHA[1]*R_PRES[2][z]                                            +ALPHA[1]*R_PRES[2][z]
# Line 507  err_t Paso_Solver_GMRES( Line 439  err_t Paso_Solver_GMRES(
439                                                    +ALPHA[5]*X_PRES[6][z];                                                    +ALPHA[5]*X_PRES[6][z];
440                }                }
441              } else if (order>6) {              } else if (order>6) {
442                #pragma omp for private(z) schedule(static)                #pragma omp parallel for private(z) schedule(static)
443                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
444                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
445                                            +ALPHA[1]*R_PRES[2][z]                                            +ALPHA[1]*R_PRES[2][z]
# Line 525  err_t Paso_Solver_GMRES( Line 457  err_t Paso_Solver_GMRES(
457                                                    +ALPHA[6]*X_PRES[7][z];                                                    +ALPHA[6]*X_PRES[7][z];
458                }                }
459                for (i=7;i<order;++i) {                for (i=7;i<order;++i) {
460                  #pragma omp for private(z) schedule(static)                  #pragma omp parallel for private(z) schedule(static)
461                  for (z=0;z<n;++z) {                  for (z=0;z<n;++z) {
462                      R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];                      R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];
463                      X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];                      X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];
# Line 536  err_t Paso_Solver_GMRES( Line 468  err_t Paso_Solver_GMRES(
468                /***                /***
469                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:
470                ***/                ***/
471                #pragma omp master                SC1=ZERO;
472                {                SC2=ZERO;
473                   SC1=ZERO;                L2_R=ZERO;
474                   SC2=ZERO;                #pragma omp parallel for private(z,diff) reduction(+:SC1,SC2) schedule(static)
                  L2_R=ZERO;  
               }  
               #pragma omp barrier  
               #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static)  
475                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
476                  diff=R_PRES[0][z]-r[z];                  diff=R_PRES[0][z]-r[z];
477                  SC1+=diff*diff;                  SC1+=diff*diff;
478                  SC2+=diff*r[z];                  SC2+=diff*r[z];
479                }                }
480                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;
481                #pragma omp for private(z) reduction(+:L2_R) schedule(static)                #pragma omp parallel for private(z) reduction(+:L2_R) schedule(static)
482                for (z=0;z<n;++z) {                for (z=0;z<n;++z) {
483                   x[z]+=gamma*(X_PRES[0][z]-x[z]);                   x[z]+=gamma*(X_PRES[0][z]-x[z]);
484                   r[z]+=gamma*(R_PRES[0][z]-r[z]);                   r[z]+=gamma*(R_PRES[0][z]-r[z]);
# Line 567  err_t Paso_Solver_GMRES( Line 495  err_t Paso_Solver_GMRES(
495           }           }
496          }          }
497          /* end of iteration */          /* end of iteration */
498          #pragma omp master          Norm_of_residual_global=norm_of_residual;
499          {      Num_iter_global=num_iter;
500             Norm_of_residual_global=norm_of_residual;      if (maxIterFlag) {
        Num_iter_global=num_iter;  
        if (maxIterFlag) {  
501             Status = SOLVER_MAXITER_REACHED;             Status = SOLVER_MAXITER_REACHED;
502         } else if (breakFlag) {         } else if (breakFlag) {
503             Status = SOLVER_BREAKDOWN;             Status = SOLVER_BREAKDOWN;
504        }      }
505          }      }
     }  /* end of parallel region */  
506      for (i=0;i<Length_of_recursion;i++) {      for (i=0;i<Length_of_recursion;i++) {
507         TMPMEMFREE(X_PRES[i]);         TMPMEMFREE(X_PRES[i]);
508         TMPMEMFREE(R_PRES[i]);         TMPMEMFREE(R_PRES[i]);
# Line 593  err_t Paso_Solver_GMRES( Line 518  err_t Paso_Solver_GMRES(
518      TMPMEMFREE(ALPHA);      TMPMEMFREE(ALPHA);
519      *iter=Num_iter_global;      *iter=Num_iter_global;
520      *tolerance=Norm_of_residual_global;      *tolerance=Norm_of_residual_global;
   }  
521    return Status;    return Status;
522  }  }

Legend:
Removed from v.1388  
changed lines
  Added in v.1556

  ViewVC Help
Powered by ViewVC 1.1.26