/[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 1758 by phornby, Fri Jul 11 13:12:46 2008 UTC revision 1759 by gross, Mon Sep 8 02:31:22 2008 UTC
# Line 72  err_t Paso_Solver_GMRES( Line 72  err_t Paso_Solver_GMRES(
72      double * tolerance,dim_t Length_of_recursion,dim_t restart,      double * tolerance,dim_t Length_of_recursion,dim_t restart,
73      Paso_Performance* pp) {      Paso_Performance* pp) {
74    
75    /* Local variables */     /* Local variables */
76    
77    double *AP,**X_PRES,**R_PRES,**P_PRES;     #ifdef _OPENMP
78           const int num_threads=omp_get_max_threads();
79       #else
80           const int num_threads=1;
81       #endif
82      double *AP,**X_PRES,**R_PRES,**P_PRES, *dots, *loc_dots;
83    double *P_PRES_dot_AP,*R_PRES_dot_P_PRES,*BREAKF,*ALPHA;    double *P_PRES_dot_AP,*R_PRES_dot_P_PRES,*BREAKF,*ALPHA;
84    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,breakf0;    double R_PRES_dot_AP0,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,breakf0;
85    double tol,Factor,sum_BREAKF,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;
86    double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES;    double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES;
87    dim_t maxit,Num_iter_global,num_iter_restart,num_iter;    dim_t maxit,Num_iter_global,num_iter_restart,num_iter;
88    dim_t i,z,order,n, Length_of_mem;    dim_t i,z,order,n, Length_of_mem, th, local_n , rest, n_start ,n_end;
89    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;
90    err_t Status=SOLVER_NO_ERROR;    err_t Status=SOLVER_NO_ERROR;
91    
# Line 99  err_t Paso_Solver_GMRES( Line 104  err_t Paso_Solver_GMRES(
104    X_PRES=TMPMEMALLOC(Length_of_mem,double*);    X_PRES=TMPMEMALLOC(Length_of_mem,double*);
105    R_PRES=TMPMEMALLOC(Length_of_mem,double*);    R_PRES=TMPMEMALLOC(Length_of_mem,double*);
106    P_PRES=TMPMEMALLOC(Length_of_mem,double*);    P_PRES=TMPMEMALLOC(Length_of_mem,double*);
107      loc_dots=TMPMEMALLOC(MAX(Length_of_mem+1,3),double);
108      dots=TMPMEMALLOC(MAX(Length_of_mem+1,3),double);
109    P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double);    P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double);
110    R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double);    R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double);
111    BREAKF=TMPMEMALLOC(Length_of_mem,double);    BREAKF=TMPMEMALLOC(Length_of_mem,double);
112    ALPHA=TMPMEMALLOC(Length_of_mem,double);    ALPHA=TMPMEMALLOC(Length_of_mem,double);
113    AP=TMPMEMALLOC(n,double);    AP=TMPMEMALLOC(n,double);
114    if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL ||    if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL ||
115           P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL) {           P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL || dots==NULL || loc_dots==NULL) {
116       Status=SOLVER_MEMORY_ERROR;       Status=SOLVER_MEMORY_ERROR;
117    } else {    } else {
118       for (i=0;i<Length_of_mem;i++) {       for (i=0;i<Length_of_mem;i++) {
# Line 125  err_t Paso_Solver_GMRES( Line 132  err_t Paso_Solver_GMRES(
132    
133        restartFlag=TRUE;        restartFlag=TRUE;
134        num_iter=0;        num_iter=0;
135        #pragma omp parallel private(i)  
136        {        #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end)
137           #pragma omp for private(z) schedule(static) nowait        for (th=0;th<num_threads;++th) {
138           for (z=0; z < n; ++z) AP[z]=0;          local_n=n/num_threads;
139           for(i=0;i<Length_of_mem;++i) {          rest=n-local_n*num_threads;
140             #pragma omp for private(z) schedule(static) nowait          n_start=local_n*th+MIN(th,rest);
141             for (z=0; z < n; ++z) {          n_end=local_n*(th+1)+MIN(th+1,rest);
142                        P_PRES[i][z]=0;          memset(&AP[n_start],0,sizeof(double)*(n_end-n_start));
143                        R_PRES[i][z]=0;          for(i=0;i<Length_of_mem;++i) {
144                        X_PRES[i][z]=0;             memset(&P_PRES[i][n_start],0,sizeof(double)*(n_end-n_start));
145             }             memset(&R_PRES[i][n_start],0,sizeof(double)*(n_end-n_start));
146           }             memset(&X_PRES[i][n_start],0,sizeof(double)*(n_end-n_start));
147            }
148        }        }
149    
150        while (! (convergeFlag || maxIterFlag || breakFlag))  {        while (! (convergeFlag || maxIterFlag || breakFlag))  {
151           if (restartFlag) {           if (restartFlag) {
152               BREAKF[0]=ONE;               BREAKF[0]=ONE;
153               #pragma omp parallel for private(z) schedule(static)               #pragma omp parallel for private(th,z, local_n, rest, n_start, n_end)
154               for (z=0; z < n; ++z) {               for (th=0;th<num_threads;++th) {
155                  R_PRES[0][z]=r[z];                 local_n=n/num_threads;
156                  X_PRES[0][z]=x[z];                 rest=n-local_n*num_threads;
157                   n_start=local_n*th+MIN(th,rest);
158                   n_end=local_n*(th+1)+MIN(th+1,rest);
159                   memcpy(&R_PRES[0][n_start],&r[n_start],sizeof(double)*(n_end-n_start));
160                   memcpy(&X_PRES[0][n_start],&x[n_start],sizeof(double)*(n_end-n_start));
161               }               }
162               num_iter_restart=0;               num_iter_restart=0;
163               restartFlag=FALSE;               restartFlag=FALSE;
# Line 166  err_t Paso_Solver_GMRES( Line 178  err_t Paso_Solver_GMRES(
178           ***** calculation of the norm of R and the scalar products of                 ***** calculation of the norm of R and the scalar products of      
179           ***   the residuals and A*P:                                                   ***   the residuals and A*P:                                        
180           ***/           ***/
181           if (order==0) {           memset(loc_dots,0,sizeof(double)*(order+1));
182              R_PRES_dot_P=ZERO;           #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end, 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)
183              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P) schedule(static)           for (th=0;th<num_threads;++th) {
184              for (z=0;z<n;++z) R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                 local_n=n/num_threads;
185              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                 rest=n-local_n*num_threads;
186           } else if (order==1) {                 n_start=local_n*th+MIN(th,rest);
187              R_PRES_dot_P=ZERO;                 n_end=local_n*(th+1)+MIN(th+1,rest);
188              P_PRES_dot_AP0=ZERO;                 if (order==0) {
189              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0) schedule(static)                     R_PRES_dot_P=ZERO;
190              for (z=0;z<n;++z) {                     #pragma ivdep
191                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
192                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
193              }                     }
194              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
195              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                     {
196           } else if (order==2) {                          loc_dots[0]+=R_PRES_dot_P;
197              R_PRES_dot_P=ZERO;                     }
198              P_PRES_dot_AP0=ZERO;                 } else if (order==1) {
199              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
200              #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1) schedule(static)                     P_PRES_dot_AP0=ZERO;
201              for (z=0;z<n;++z) {                     #pragma ivdep
202                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
203                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
204                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
205              }                     }
206              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
207              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
208              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[0]+=R_PRES_dot_P;
209           } else if (order==3) {                          loc_dots[1]+=P_PRES_dot_AP0;
210              R_PRES_dot_P=ZERO;                     }
211              P_PRES_dot_AP0=ZERO;                 } else if (order==2) {
212              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
213              P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP0=ZERO;
214              #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_AP1=ZERO;
215              for (z=0;z<n;++z) {                     #pragma ivdep
216                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
217                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
218                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
219                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
220              }                     }
221              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
222              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
223              P_PRES_dot_AP[2]=P_PRES_dot_AP2;                          loc_dots[0]+=R_PRES_dot_P;
224              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[1]+=P_PRES_dot_AP0;
225           } else if (order==4) {                          loc_dots[2]+=P_PRES_dot_AP1;
226              R_PRES_dot_P=ZERO;                     }
227              P_PRES_dot_AP0=ZERO;                 } else if (order==3) {
228              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
229              P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP0=ZERO;
230              P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP1=ZERO;
231              #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_AP2=ZERO;
232              for (z=0;z<n;++z) {                     #pragma ivdep
233                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
234                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
235                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
236                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
237                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
238              }                     }
239              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
240              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
241              P_PRES_dot_AP[2]=P_PRES_dot_AP2;                          loc_dots[0]+=R_PRES_dot_P;
242              P_PRES_dot_AP[3]=P_PRES_dot_AP3;                          loc_dots[1]+=P_PRES_dot_AP0;
243              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[2]+=P_PRES_dot_AP1;
244           } else if (order==5) {                          loc_dots[3]+=P_PRES_dot_AP2;
245              R_PRES_dot_P=ZERO;                     }
246              P_PRES_dot_AP0=ZERO;                 } else if (order==4) {
247              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
248              P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP0=ZERO;
249              P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP1=ZERO;
250              P_PRES_dot_AP4=ZERO;                     P_PRES_dot_AP2=ZERO;
251              #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_AP3=ZERO;
252              for (z=0;z<n;++z) {                     #pragma ivdep
253                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
254                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
255                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
256                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
257                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
258                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
259              }                     }
260              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
261              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
262              P_PRES_dot_AP[2]=P_PRES_dot_AP2;                          loc_dots[0]+=R_PRES_dot_P;
263              P_PRES_dot_AP[3]=P_PRES_dot_AP3;                          loc_dots[1]+=P_PRES_dot_AP0;
264              P_PRES_dot_AP[4]=P_PRES_dot_AP4;                          loc_dots[2]+=P_PRES_dot_AP1;
265              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[3]+=P_PRES_dot_AP2;
266           } else if (order==6) {                          loc_dots[4]+=P_PRES_dot_AP3;
267              R_PRES_dot_P=ZERO;                     }
268              P_PRES_dot_AP0=ZERO;                 } else if (order==5) {
269              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
270              P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP0=ZERO;
271              P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP1=ZERO;
272              P_PRES_dot_AP4=ZERO;                     P_PRES_dot_AP2=ZERO;
273              P_PRES_dot_AP5=ZERO;                     P_PRES_dot_AP3=ZERO;
274              #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_AP4=ZERO;
275              for (z=0;z<n;++z) {                     #pragma ivdep
276                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
277                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
278                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
279                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
280                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
281                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
282                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];                         P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
283              }                     }
284              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
285              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
286              P_PRES_dot_AP[2]=P_PRES_dot_AP2;                          loc_dots[0]+=R_PRES_dot_P;
287              P_PRES_dot_AP[3]=P_PRES_dot_AP3;                          loc_dots[1]+=P_PRES_dot_AP0;
288              P_PRES_dot_AP[4]=P_PRES_dot_AP4;                          loc_dots[2]+=P_PRES_dot_AP1;
289              P_PRES_dot_AP[5]=P_PRES_dot_AP5;                          loc_dots[3]+=P_PRES_dot_AP2;
290              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[4]+=P_PRES_dot_AP3;
291           } else if (order>6) {                          loc_dots[5]+=P_PRES_dot_AP4;
292              R_PRES_dot_P=ZERO;                     }
293              P_PRES_dot_AP0=ZERO;                 } else if (order==6) {
294              P_PRES_dot_AP1=ZERO;                     R_PRES_dot_P=ZERO;
295              P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP0=ZERO;
296              P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP1=ZERO;
297              P_PRES_dot_AP4=ZERO;                     P_PRES_dot_AP2=ZERO;
298              P_PRES_dot_AP5=ZERO;                     P_PRES_dot_AP3=ZERO;
299              P_PRES_dot_AP6=ZERO;                     P_PRES_dot_AP4=ZERO;
300              #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_AP5=ZERO;
301              for (z=0;z<n;++z) {                     #pragma ivdep
302                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     for (z=n_start; z < n_end; ++z) {
303                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
304                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
305                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
306                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
307                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
308                  P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];                         P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
309                  P_PRES_dot_AP6+=P_PRES[6][z]*AP[z];                         P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
310              }                     }
311              P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     #pragma omp critical
312              P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     {
313              P_PRES_dot_AP[2]=P_PRES_dot_AP2;                          loc_dots[0]+=R_PRES_dot_P;
314              P_PRES_dot_AP[3]=P_PRES_dot_AP3;                          loc_dots[1]+=P_PRES_dot_AP0;
315              P_PRES_dot_AP[4]=P_PRES_dot_AP4;                          loc_dots[2]+=P_PRES_dot_AP1;
316              P_PRES_dot_AP[5]=P_PRES_dot_AP5;                          loc_dots[3]+=P_PRES_dot_AP2;
317              P_PRES_dot_AP[6]=P_PRES_dot_AP6;                          loc_dots[4]+=P_PRES_dot_AP3;
318              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[5]+=P_PRES_dot_AP4;
319                            loc_dots[6]+=P_PRES_dot_AP5;
320              P_PRES_dot_AP0=ZERO;                     }
321              for (i=7;i<order;++i) {                 } else {
322                  #pragma omp parallel for private(z) reduction(+:P_PRES_dot_AP0) schedule(static)                    R_PRES_dot_P=ZERO;
323                  for (z=0;z<n;++z) P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];                    P_PRES_dot_AP0=ZERO;
324                  P_PRES_dot_AP[i]=P_PRES_dot_AP0;                    P_PRES_dot_AP1=ZERO;
325                  P_PRES_dot_AP0=ZERO;                    P_PRES_dot_AP2=ZERO;
326              }                    P_PRES_dot_AP3=ZERO;
327                      P_PRES_dot_AP4=ZERO;
328                      P_PRES_dot_AP5=ZERO;
329                      P_PRES_dot_AP6=ZERO;
330                      #pragma ivdep
331                      for (z=n_start; z < n_end; ++z) {
332                          R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
333                          P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
334                          P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
335                          P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
336                          P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
337                          P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
338                          P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
339                          P_PRES_dot_AP6+=P_PRES[6][z]*AP[z];
340                      }
341                      #pragma omp critical
342                      {
343                           loc_dots[0]+=R_PRES_dot_P;
344                           loc_dots[1]+=P_PRES_dot_AP0;
345                           loc_dots[2]+=P_PRES_dot_AP1;
346                           loc_dots[3]+=P_PRES_dot_AP2;
347                           loc_dots[4]+=P_PRES_dot_AP3;
348                           loc_dots[5]+=P_PRES_dot_AP4;
349                           loc_dots[6]+=P_PRES_dot_AP5;
350                           loc_dots[7]+=P_PRES_dot_AP6;
351                      }
352                      for (i=7;i<order;++i) {
353                           P_PRES_dot_AP0=ZERO;
354                           #pragma ivdep
355                           for (z=n_start; z < n_end; ++z) {
356                                 P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];
357                           }
358                           #pragma omp critical
359                           {
360                               loc_dots[i+1]+=P_PRES_dot_AP0;
361                           }
362                      }
363                   }
364           }           }
365           /* this fixes a problem with the intel compiler */           #ifdef PASO_MPI
366           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];                  MPI_Allreduce(loc_dots, dots, order+1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
367                    R_PRES_dot_P_PRES[0]=dots[0];
368                    memcpy(P_PRES_dot_AP,&dots[1],sizeof(double)*order);
369             #else
370                    R_PRES_dot_P_PRES[0]=loc_dots[0];
371                    memcpy(P_PRES_dot_AP,&loc_dots[1],sizeof(double)*order);
372             #endif
373             R_PRES_dot_AP0=R_PRES_dot_P_PRES[0];
374           /***   if sum_BREAKF is equal to zero a breakdown occurs.           /***   if sum_BREAKF is equal to zero a breakdown occurs.
375            ***   iteration procedure can be continued but R_PRES is not the            ***   iteration procedure can be continued but R_PRES is not the
376            ***   residual of X_PRES approximation.            ***   residual of X_PRES approximation.
377            ***/            ***/
378           sum_BREAKF=0.;           sum_BREAKF=0.;
379             #pragma ivdep
380           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
381           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));           breakFlag=!((ABS(R_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));
382           if (!breakFlag) {           if (!breakFlag) {
383              breakFlag=FALSE;              breakFlag=FALSE;
384              /***              /***
385              ***** X_PRES and R_PRES are moved to memory:              ***** X_PRES and R_PRES are moved to memory:
386              ***/              ***/
387              Factor=0.;              Factor=0.;
388                #pragma ivdep
389              for (i=0;i<order;++i) {              for (i=0;i<order;++i) {
390                     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];
391                     Factor+=BREAKF[i]*ALPHA[i];                     Factor+=BREAKF[i]*ALPHA[i];
# Line 337  err_t Paso_Solver_GMRES( Line 395  err_t Paso_Solver_GMRES(
395               save_R_PRES=R_PRES[Length_of_mem-1];               save_R_PRES=R_PRES[Length_of_mem-1];
396               save_XPRES=X_PRES[Length_of_mem-1];               save_XPRES=X_PRES[Length_of_mem-1];
397               save_P_PRES=P_PRES[Length_of_mem-1];               save_P_PRES=P_PRES[Length_of_mem-1];
398                 #pragma ivdep
399               for (i=Length_of_mem-1;i>0;--i) {               for (i=Length_of_mem-1;i>0;--i) {
400                     BREAKF[i]=BREAKF[i-1];                     BREAKF[i]=BREAKF[i-1];
401                     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];
# Line 366  err_t Paso_Solver_GMRES( Line 425  err_t Paso_Solver_GMRES(
425              ***                                                                              ***                                                                
426              **/              **/
427              breakf0=BREAKF[0];              breakf0=BREAKF[0];
428              if (order==0) {              #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end)
429                #pragma omp parallel for private(z) schedule(static)              for (th=0;th<num_threads;++th) {
430                for (z=0;z<n;++z) {                 local_n=n/num_threads;
431                    R_PRES[0][z]= Factor*       AP[z];                 rest=n-local_n*num_threads;
432                    X_PRES[0][z]=-Factor*P_PRES[1][z];                 n_start=local_n*th+MIN(th,rest);
433                }                 n_end=local_n*(th+1)+MIN(th+1,rest);
434              } else if (order==1) {                 if (order==0) {
435                #pragma omp parallel for private(z) schedule(static)                     #pragma ivdep
436                for (z=0;z<n;++z) {                     for (z=n_start; z < n_end; ++z) {
437                    R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z];                        R_PRES[0][z]= Factor*       AP[z];
438                    X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];                        X_PRES[0][z]=-Factor*P_PRES[1][z];
439                }                     }
440              } else if (order==2) {                 } else if (order==1) {
441                #pragma omp parallel for private(z) schedule(static)                     #pragma ivdep
442                for (z=0;z<n;++z) {                     for (z=n_start; z < n_end; ++z) {
443                   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];
444                                                    +ALPHA[1]*R_PRES[2][z];                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];
445                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     }
446                                                    +ALPHA[1]*X_PRES[2][z];                 } else if (order==2) {
447                }                     #pragma ivdep
448              } else if (order==3) {                     for (z=n_start; z < n_end; ++z) {
449                #pragma omp parallel for private(z) schedule(static)                         R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z]
450                for (z=0;z<n;++z) {                                                          +ALPHA[1]*R_PRES[2][z];
451                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
452                                             +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[1]*X_PRES[2][z];
453                                             +ALPHA[2]*R_PRES[3][z];                     }
454                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                 } else if (order==3) {
455                                                    +ALPHA[1]*X_PRES[2][z]                     #pragma ivdep
456                                                    +ALPHA[2]*X_PRES[3][z];                     for (z=n_start; z < n_end; ++z) {
457                }                         R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
458              } else if (order==4) {                                                   +ALPHA[1]*R_PRES[2][z]
459                #pragma omp parallel for private(z) schedule(static)                                                   +ALPHA[2]*R_PRES[3][z];
460                for (z=0;z<n;++z) {                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
461                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[1]*X_PRES[2][z]
462                                             +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[2]*X_PRES[3][z];
463                                             +ALPHA[2]*R_PRES[3][z]                     }
464                                             +ALPHA[3]*R_PRES[4][z];                 } else if (order==4) {
465                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     #pragma ivdep
466                                                    +ALPHA[1]*X_PRES[2][z]                     for (z=n_start; z < n_end; ++z) {
467                                                    +ALPHA[2]*X_PRES[3][z]                         R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
468                                                    +ALPHA[3]*X_PRES[4][z];                                                   +ALPHA[1]*R_PRES[2][z]
469                }                                                   +ALPHA[2]*R_PRES[3][z]
470              } else if (order==5) {                                                   +ALPHA[3]*R_PRES[4][z];
471                #pragma omp parallel for private(z) schedule(static)                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
472                for (z=0;z<n;++z) {                                                          +ALPHA[1]*X_PRES[2][z]
473                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[2]*X_PRES[3][z]
474                                            +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[3]*X_PRES[4][z];
475                                            +ALPHA[2]*R_PRES[3][z]                     }
476                                            +ALPHA[3]*R_PRES[4][z]                 } else if (order==5) {
477                                            +ALPHA[4]*R_PRES[5][z];                     #pragma ivdep
478                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     for (z=n_start; z < n_end; ++z) {
479                                                    +ALPHA[1]*X_PRES[2][z]                         R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
480                                                    +ALPHA[2]*X_PRES[3][z]                                                  +ALPHA[1]*R_PRES[2][z]
481                                                    +ALPHA[3]*X_PRES[4][z]                                                  +ALPHA[2]*R_PRES[3][z]
482                                                    +ALPHA[4]*X_PRES[5][z];                                                  +ALPHA[3]*R_PRES[4][z]
483                }                                                  +ALPHA[4]*R_PRES[5][z];
484              } else if (order==6) {                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
485                #pragma omp parallel for private(z) schedule(static)                                                          +ALPHA[1]*X_PRES[2][z]
486                for (z=0;z<n;++z) {                                                          +ALPHA[2]*X_PRES[3][z]
487                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[3]*X_PRES[4][z]
488                                            +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[4]*X_PRES[5][z];
489                                            +ALPHA[2]*R_PRES[3][z]                      }
490                                            +ALPHA[3]*R_PRES[4][z]                 } else if (order==6) {
491                                            +ALPHA[4]*R_PRES[5][z]                     #pragma ivdep
492                                            +ALPHA[5]*R_PRES[6][z];                     for (z=n_start; z < n_end; ++z) {
493                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                        R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
494                                                    +ALPHA[1]*X_PRES[2][z]                                                 +ALPHA[1]*R_PRES[2][z]
495                                                    +ALPHA[2]*X_PRES[3][z]                                                 +ALPHA[2]*R_PRES[3][z]
496                                                    +ALPHA[3]*X_PRES[4][z]                                                 +ALPHA[3]*R_PRES[4][z]
497                                                    +ALPHA[4]*X_PRES[5][z]                                                 +ALPHA[4]*R_PRES[5][z]
498                                                    +ALPHA[5]*X_PRES[6][z];                                                 +ALPHA[5]*R_PRES[6][z];
499                }                        X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
500              } else if (order>6) {                                                         +ALPHA[1]*X_PRES[2][z]
501                #pragma omp parallel for private(z) schedule(static)                                                         +ALPHA[2]*X_PRES[3][z]
502                for (z=0;z<n;++z) {                                                         +ALPHA[3]*X_PRES[4][z]
503                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                         +ALPHA[4]*X_PRES[5][z]
504                                            +ALPHA[1]*R_PRES[2][z]                                                         +ALPHA[5]*X_PRES[6][z];
505                                            +ALPHA[2]*R_PRES[3][z]                     }
506                                            +ALPHA[3]*R_PRES[4][z]                 } else {
507                                            +ALPHA[4]*R_PRES[5][z]                     #pragma ivdep
508                                            +ALPHA[5]*R_PRES[6][z]                     for (z=n_start; z < n_end; ++z) {
509                                            +ALPHA[6]*R_PRES[7][z];                          R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
510                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                                                   +ALPHA[1]*R_PRES[2][z]
511                                                    +ALPHA[1]*X_PRES[2][z]                                                   +ALPHA[2]*R_PRES[3][z]
512                                                    +ALPHA[2]*X_PRES[3][z]                                                   +ALPHA[3]*R_PRES[4][z]
513                                                    +ALPHA[3]*X_PRES[4][z]                                                   +ALPHA[4]*R_PRES[5][z]
514                                                    +ALPHA[4]*X_PRES[5][z]                                                   +ALPHA[5]*R_PRES[6][z]
515                                                    +ALPHA[5]*X_PRES[6][z]                                                   +ALPHA[6]*R_PRES[7][z];
516                                                    +ALPHA[6]*X_PRES[7][z];                          X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
517                }                                                           +ALPHA[1]*X_PRES[2][z]
518                for (i=7;i<order;++i) {                                                           +ALPHA[2]*X_PRES[3][z]
519                  #pragma omp parallel for private(z) schedule(static)                                                           +ALPHA[3]*X_PRES[4][z]
520                  for (z=0;z<n;++z) {                                                           +ALPHA[4]*X_PRES[5][z]
521                      R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];                                                           +ALPHA[5]*X_PRES[6][z]
522                      X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];                                                           +ALPHA[6]*X_PRES[7][z];
523                  }                     }
524                }                     for (i=7;i<order;++i) {
525                          #pragma ivdep
526                          for (z=n_start; z < n_end; ++z) {
527                             R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];
528                             X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];
529                         }
530                       }
531                   }
532              }              }
533              if (breakf0>0.) {              if (breakf0>0.) {
534                /***                /***
535                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:
536                ***/                ***/
537                SC1=ZERO;                memset(loc_dots,0,sizeof(double)*3);
538                SC2=ZERO;                #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end,diff,SC1,SC2)
539                L2_R=ZERO;                for (th=0;th<num_threads;++th) {
540                #pragma omp parallel for private(z,diff) reduction(+:SC1,SC2) schedule(static)                    local_n=n/num_threads;
541                for (z=0;z<n;++z) {                    rest=n-local_n*num_threads;
542                  diff=R_PRES[0][z]-r[z];                    n_start=local_n*th+MIN(th,rest);
543                  SC1+=diff*diff;                    n_end=local_n*(th+1)+MIN(th+1,rest);
544                  SC2+=diff*r[z];                    SC1=ZERO;
545                }                    SC2=ZERO;
546                      #pragma ivdep
547                      for (z=n_start; z < n_end; ++z) {
548                           diff=R_PRES[0][z]-r[z];
549                           SC1+=diff*diff;
550                           SC2+=diff*r[z];
551                      }
552                      #pragma omp critical
553                      {
554                           loc_dots[0]+=SC1;
555                           loc_dots[1]+=SC2;
556                      }
557                  }
558                  #ifdef PASO_MPI
559                      MPI_Allreduce(loc_dots, dots, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
560                      SC1=dots[0];
561                      SC2=dots[1];
562                  #else
563                      SC1=loc_dots[0];
564                      SC2=loc_dots[1];
565                  #endif
566                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;
567                #pragma omp parallel for private(z) reduction(+:L2_R) schedule(static)                #pragma omp parallel for private(th,z,local_n, rest, n_start, n_end,diff,L2_R)
568                for (z=0;z<n;++z) {                for (th=0;th<num_threads;++th) {
569                   x[z]+=gamma*(X_PRES[0][z]-x[z]);                    local_n=n/num_threads;
570                   r[z]+=gamma*(R_PRES[0][z]-r[z]);                    rest=n-local_n*num_threads;
571                   L2_R+=r[z]*r[z];                    n_start=local_n*th+MIN(th,rest);
572                }                    n_end=local_n*(th+1)+MIN(th+1,rest);
573                      L2_R=ZERO;
574                      #pragma ivdep
575                      for (z=n_start; z < n_end; ++z) {
576                          x[z]+=gamma*(X_PRES[0][z]-x[z]);
577                          r[z]+=gamma*(R_PRES[0][z]-r[z]);
578                          L2_R+=r[z]*r[z];
579                      }
580                      #pragma omp critical
581                      {
582                           loc_dots[2]+=L2_R;
583                      }
584                  }
585                  #ifdef PASO_MPI
586                      MPI_Allreduce(&loc_dots[2], &dots[2], 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
587                      L2_R=dots[2];
588                  #else
589                      L2_R=loc_dots[2];
590                  #endif
591                norm_of_residual=sqrt(L2_R);                norm_of_residual=sqrt(L2_R);
592                convergeFlag = (norm_of_residual <= tol);                convergeFlag = (norm_of_residual <= tol);
593                if (restart>0) restartFlag=(num_iter_restart >= restart);                if (restart>0) restartFlag=(num_iter_restart >= restart);
# Line 516  err_t Paso_Solver_GMRES( Line 620  err_t Paso_Solver_GMRES(
620      TMPMEMFREE(R_PRES_dot_P_PRES);      TMPMEMFREE(R_PRES_dot_P_PRES);
621      TMPMEMFREE(BREAKF);      TMPMEMFREE(BREAKF);
622      TMPMEMFREE(ALPHA);      TMPMEMFREE(ALPHA);
623        TMPMEMFREE(dots);
624        TMPMEMFREE(loc_dots);
625      *iter=Num_iter_global;      *iter=Num_iter_global;
626      *tolerance=Norm_of_residual_global;      *tolerance=Norm_of_residual_global;
627    return Status;    return Status;
Legend:
Removed from v.1758  
changed lines
  Added in v.1759
  ViewVC Help
Powered by ViewVC 1.1.26