/[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/paso/src/Solvers/GMRES.c revision 415 by gross, Wed Jan 4 05:37:33 2006 UTC trunk/paso/src/GMRES.c revision 1767 by gross, Mon Sep 8 02:53:50 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /*  /*
17  *  Purpose  *  Purpose
18  *  =======  *  =======
# Line 56  err_t Paso_Solver_GMRES( Line 69  err_t Paso_Solver_GMRES(
69      double * r,      double * r,
70      double * x,      double * x,
71      dim_t *iter,      dim_t *iter,
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) {
74    
75    /* Local variables */     /* Local variables */
76    
77    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];     #ifdef _OPENMP
78    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)];         const int num_threads=omp_get_max_threads();
79    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;     #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;
84      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;    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    
92      
93    /* adapt original routine parameters */    /* adapt original routine parameters */
94      n = Paso_SystemMatrix_getTotalNumRows(A);
95    dim_t n=A->num_cols * A-> col_block_size;    Length_of_mem=MAX(Length_of_recursion,0)+1;
   dim_t Length_of_mem=MAX(Length_of_recursion,0)+1;  
96    
97    /*     Test the input parameters. */    /*     Test the input parameters. */
98    if (restart>0) restart=MAX(Length_of_recursion,restart);    if (restart>0) restart=MAX(Length_of_recursion,restart);
# Line 82  err_t Paso_Solver_GMRES( Line 101  err_t Paso_Solver_GMRES(
101    }    }
102    
103    /*     allocate memory: */    /*     allocate memory: */
104      X_PRES=TMPMEMALLOC(Length_of_mem,double*);
105      R_PRES=TMPMEMALLOC(Length_of_mem,double*);
106      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);
110      R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double);
111      BREAKF=TMPMEMALLOC(Length_of_mem,double);
112      ALPHA=TMPMEMALLOC(Length_of_mem,double);
113    AP=TMPMEMALLOC(n,double);    AP=TMPMEMALLOC(n,double);
114    if (AP==NULL) Status=SOLVER_MEMORY_ERROR;    if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL ||
115    for (i=0;i<Length_of_mem;i++) {           P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL || dots==NULL || loc_dots==NULL) {
116      X_PRES[i]=TMPMEMALLOC(n,double);       Status=SOLVER_MEMORY_ERROR;
117      R_PRES[i]=TMPMEMALLOC(n,double);    } else {
118      P_PRES[i]=TMPMEMALLOC(n,double);       for (i=0;i<Length_of_mem;i++) {
119      if (X_PRES[i]==NULL || R_PRES[i]==NULL ||  P_PRES[i]==NULL) Status=SOLVER_MEMORY_ERROR;         X_PRES[i]=TMPMEMALLOC(n,double);
120           R_PRES[i]=TMPMEMALLOC(n,double);
121           P_PRES[i]=TMPMEMALLOC(n,double);
122           if (X_PRES[i]==NULL || R_PRES[i]==NULL ||  P_PRES[i]==NULL) Status=SOLVER_MEMORY_ERROR;
123         }
124    }    }
125    if ( Status ==SOLVER_NO_ERROR ) {    if ( Status ==SOLVER_NO_ERROR ) {
126    
# Line 97  err_t Paso_Solver_GMRES( Line 128  err_t Paso_Solver_GMRES(
128      maxit=*iter;      maxit=*iter;
129      tol=*tolerance;      tol=*tolerance;
130    
131      #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 */  
132    
133        restartFlag=TRUE;        restartFlag=TRUE;
134        num_iter=0;        num_iter=0;
135        #pragma omp for private(z) schedule(static) nowait  
136        for (z=0; z < n; ++z) AP[z]=0;        #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end)
137        for(i=0;i<Length_of_mem;++i) {        for (th=0;th<num_threads;++th) {
138          #pragma omp for private(z) schedule(static) nowait          local_n=n/num_threads;
139          for (z=0; z < n; ++z) {          rest=n-local_n*num_threads;
140                     P_PRES[i][z]=0;          n_start=local_n*th+MIN(th,rest);
141                     R_PRES[i][z]=0;          n_end=local_n*(th+1)+MIN(th+1,rest);
142                     X_PRES[i][z]=0;          memset(&AP[n_start],0,sizeof(double)*(n_end-n_start));
143            for(i=0;i<Length_of_mem;++i) {
144               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))  {
          #pragma omp barrier  
151           if (restartFlag) {           if (restartFlag) {
              #pragma omp master  
152               BREAKF[0]=ONE;               BREAKF[0]=ONE;
153               #pragma omp for private(z) schedule(static) nowait               #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 145  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              #pragma omp master           #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              R_PRES_dot_P=ZERO;           for (th=0;th<num_threads;++th) {
184              #pragma omp barrier                 local_n=n/num_threads;
185              #pragma omp for private(z) reduction(+:R_PRES_dot_P) schedule(static)                 rest=n-local_n*num_threads;
186              for (z=0;z<n;++z)                 n_start=local_n*th+MIN(th,rest);
187                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                 n_end=local_n*(th+1)+MIN(th+1,rest);
188              #pragma omp master                 if (order==0) {
189              R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                     R_PRES_dot_P=ZERO;
190           } else if (order==1) {                     #pragma ivdep
191              #pragma omp master                     for (z=n_start; z < n_end; ++z) {
192              {                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
193                 R_PRES_dot_P=ZERO;                     }
194                 P_PRES_dot_AP0=ZERO;                     #pragma omp critical
195              }                     {
196              #pragma omp barrier                          loc_dots[0]+=R_PRES_dot_P;
197              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0) schedule(static)                     }
198              for (z=0;z<n;++z) {                 } else if (order==1) {
199                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     R_PRES_dot_P=ZERO;
200                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                     P_PRES_dot_AP0=ZERO;
201              }                     #pragma ivdep
202              #pragma omp master                     for (z=n_start; z < n_end; ++z) {
203              {                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
204                P_PRES_dot_AP[0]=P_PRES_dot_AP0;                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
205                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                     }
206              }                     #pragma omp critical
207           } else if (order==2) {                     {
208              #pragma omp master                          loc_dots[0]+=R_PRES_dot_P;
209              {                          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_AP0=ZERO;
214              #pragma omp barrier                     P_PRES_dot_AP1=ZERO;
215              #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1) schedule(static)                     #pragma ivdep
216              for (z=0;z<n;++z) {                     for (z=n_start; z < n_end; ++z) {
217                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
218                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
219                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
220              }                     }
221              #pragma omp master                     #pragma omp critical
222              {                     {
223                P_PRES_dot_AP[0]=P_PRES_dot_AP0;                          loc_dots[0]+=R_PRES_dot_P;
224                P_PRES_dot_AP[1]=P_PRES_dot_AP1;                          loc_dots[1]+=P_PRES_dot_AP0;
225                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                          loc_dots[2]+=P_PRES_dot_AP1;
226              }                     }
227           } else if (order==3) {                 } else if (order==3) {
228              #pragma omp master                     R_PRES_dot_P=ZERO;
229              {                     P_PRES_dot_AP0=ZERO;
230                 R_PRES_dot_P=ZERO;                     P_PRES_dot_AP1=ZERO;
231                 P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP2=ZERO;
232                 P_PRES_dot_AP1=ZERO;                     #pragma ivdep
233                 P_PRES_dot_AP2=ZERO;                     for (z=n_start; z < n_end; ++z) {
234              }                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
235              #pragma omp barrier                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
236              #pragma omp 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+=P_PRES[1][z]*AP[z];
237              for (z=0;z<n;++z) {                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
238                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                     }
239                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                     #pragma omp critical
240                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                     {
241                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                          loc_dots[0]+=R_PRES_dot_P;
242              }                          loc_dots[1]+=P_PRES_dot_AP0;
243              #pragma omp master                          loc_dots[2]+=P_PRES_dot_AP1;
244              {                          loc_dots[3]+=P_PRES_dot_AP2;
245                P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     }
246                P_PRES_dot_AP[1]=P_PRES_dot_AP1;                 } else if (order==4) {
247                P_PRES_dot_AP[2]=P_PRES_dot_AP2;                     R_PRES_dot_P=ZERO;
248                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                     P_PRES_dot_AP0=ZERO;
249              }                     P_PRES_dot_AP1=ZERO;
250           } else if (order==4) {                     P_PRES_dot_AP2=ZERO;
251              #pragma omp master                     P_PRES_dot_AP3=ZERO;
252              {                     #pragma ivdep
253                R_PRES_dot_P=ZERO;                     for (z=n_start; z < n_end; ++z) {
254                P_PRES_dot_AP0=ZERO;                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
255                P_PRES_dot_AP1=ZERO;                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
256                P_PRES_dot_AP2=ZERO;                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
257                P_PRES_dot_AP3=ZERO;                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
258              }                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
259              #pragma omp barrier                     }
260              #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)                     #pragma omp critical
261              for (z=0;z<n;++z) {                     {
262                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                          loc_dots[0]+=R_PRES_dot_P;
263                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                          loc_dots[1]+=P_PRES_dot_AP0;
264                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                          loc_dots[2]+=P_PRES_dot_AP1;
265                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                          loc_dots[3]+=P_PRES_dot_AP2;
266                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                          loc_dots[4]+=P_PRES_dot_AP3;
267              }                     }
268              #pragma omp master                 } else if (order==5) {
269              {                     R_PRES_dot_P=ZERO;
270                P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     P_PRES_dot_AP0=ZERO;
271                P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     P_PRES_dot_AP1=ZERO;
272                P_PRES_dot_AP[2]=P_PRES_dot_AP2;                     P_PRES_dot_AP2=ZERO;
273                P_PRES_dot_AP[3]=P_PRES_dot_AP3;                     P_PRES_dot_AP3=ZERO;
274                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                     P_PRES_dot_AP4=ZERO;
275              }                     #pragma ivdep
276           } else if (order==5) {                     for (z=n_start; z < n_end; ++z) {
277              #pragma omp master                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
278              {                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
279                R_PRES_dot_P=ZERO;                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
280                P_PRES_dot_AP0=ZERO;                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
281                P_PRES_dot_AP1=ZERO;                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
282                P_PRES_dot_AP2=ZERO;                         P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
283                P_PRES_dot_AP3=ZERO;                     }
284                P_PRES_dot_AP4=ZERO;                     #pragma omp critical
285              }                     {
286              #pragma omp barrier                          loc_dots[0]+=R_PRES_dot_P;
287              #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)                          loc_dots[1]+=P_PRES_dot_AP0;
288              for (z=0;z<n;++z) {                          loc_dots[2]+=P_PRES_dot_AP1;
289                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                          loc_dots[3]+=P_PRES_dot_AP2;
290                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                          loc_dots[4]+=P_PRES_dot_AP3;
291                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                          loc_dots[5]+=P_PRES_dot_AP4;
292                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                     }
293                  P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];                 } else if (order==6) {
294                  P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];                     R_PRES_dot_P=ZERO;
295              }                     P_PRES_dot_AP0=ZERO;
296              #pragma omp master                     P_PRES_dot_AP1=ZERO;
297              {                     P_PRES_dot_AP2=ZERO;
298                P_PRES_dot_AP[0]=P_PRES_dot_AP0;                     P_PRES_dot_AP3=ZERO;
299                P_PRES_dot_AP[1]=P_PRES_dot_AP1;                     P_PRES_dot_AP4=ZERO;
300                P_PRES_dot_AP[2]=P_PRES_dot_AP2;                     P_PRES_dot_AP5=ZERO;
301                P_PRES_dot_AP[3]=P_PRES_dot_AP3;                     #pragma ivdep
302                P_PRES_dot_AP[4]=P_PRES_dot_AP4;                     for (z=n_start; z < n_end; ++z) {
303                R_PRES_dot_P_PRES[0]=R_PRES_dot_P;                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
304              }                         P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];
305           } else if (order==6) {                         P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];
306              #pragma omp master                         P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];
307              {                         P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];
308                R_PRES_dot_P=ZERO;                         P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];
309                P_PRES_dot_AP0=ZERO;                         P_PRES_dot_AP5+=P_PRES[5][z]*AP[z];
310                P_PRES_dot_AP1=ZERO;                     }
311                P_PRES_dot_AP2=ZERO;                     #pragma omp critical
312                P_PRES_dot_AP3=ZERO;                     {
313                P_PRES_dot_AP4=ZERO;                          loc_dots[0]+=R_PRES_dot_P;
314                P_PRES_dot_AP5=ZERO;                          loc_dots[1]+=P_PRES_dot_AP0;
315              }                          loc_dots[2]+=P_PRES_dot_AP1;
316              #pragma omp barrier                          loc_dots[3]+=P_PRES_dot_AP2;
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) schedule(static)                          loc_dots[4]+=P_PRES_dot_AP3;
318              for (z=0;z<n;++z) {                          loc_dots[5]+=P_PRES_dot_AP4;
319                  R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                          loc_dots[6]+=P_PRES_dot_AP5;
320                  P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];                     }
321                  P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];                 } else {
322                  P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];                    R_PRES_dot_P=ZERO;
                 P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];  
                 P_PRES_dot_AP4+=P_PRES[4][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;  
             }  
          } else if (order>6) {  
             #pragma omp master  
             {  
               R_PRES_dot_P=ZERO;  
               P_PRES_dot_AP0=ZERO;  
               P_PRES_dot_AP1=ZERO;  
               P_PRES_dot_AP2=ZERO;  
               P_PRES_dot_AP3=ZERO;  
               P_PRES_dot_AP4=ZERO;  
               P_PRES_dot_AP5=ZERO;  
               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)  
             for (z=0;z<n;++z) {  
                 R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];  
                 P_PRES_dot_AP0+=P_PRES[0][z]*AP[z];  
                 P_PRES_dot_AP1+=P_PRES[1][z]*AP[z];  
                 P_PRES_dot_AP2+=P_PRES[2][z]*AP[z];  
                 P_PRES_dot_AP3+=P_PRES[3][z]*AP[z];  
                 P_PRES_dot_AP4+=P_PRES[4][z]*AP[z];  
                 P_PRES_dot_AP5+=P_PRES[5][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;  
             }  
             for (i=7;i<order;++i) {  
                 #pragma omp barrier  
                 #pragma omp for private(z) reduction(+:P_PRES_dot_AP0) schedule(static)  
                 for (z=0;z<n;++z) P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];  
                 #pragma omp master  
                 {  
                   P_PRES_dot_AP[i]=P_PRES_dot_AP0;  
323                    P_PRES_dot_AP0=ZERO;                    P_PRES_dot_AP0=ZERO;
324                  }                    P_PRES_dot_AP1=ZERO;
325              }                    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           #pragma omp master                  MPI_Allreduce(loc_dots, dots, order+1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
367           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];                  R_PRES_dot_P_PRES[0]=dots[0];
368           #pragma omp barrier                  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    
375           /***   if sum_BREAKF is equal to zero a breakdown occurs.           /***   if sum_BREAKF is equal to zero a breakdown occurs.
376            ***   iteration procedure can be continued but R_PRES is not the            ***   iteration procedure can be continued but R_PRES is not the
377            ***   residual of X_PRES approximation.            ***   residual of X_PRES approximation.
378            ***/            ***/
379           sum_BREAKF=0.;           sum_BREAKF=0.;
380             #pragma ivdep
381           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
382           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));           breakFlag=!((ABS(R_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));
383           if (!breakFlag) {           if (!breakFlag) {
384              breakFlag=FALSE;              breakFlag=FALSE;
385              /***              /***
386              ***** X_PRES and R_PRES are moved to memory:              ***** X_PRES and R_PRES are moved to memory:
387              ***/              ***/
388              #pragma omp master              Factor=0.;
389              {              #pragma ivdep
390                 Factor=0.;              for (i=0;i<order;++i) {
                for (i=0;i<order;++i) {  
391                     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];
392                     Factor+=BREAKF[i]*ALPHA[i];                     Factor+=BREAKF[i]*ALPHA[i];
393                 }              }
394    
395                 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];
396                 save_R_PRES=R_PRES[Length_of_mem-1];               save_R_PRES=R_PRES[Length_of_mem-1];
397                 save_XPRES=X_PRES[Length_of_mem-1];               save_XPRES=X_PRES[Length_of_mem-1];
398                 save_P_PRES=P_PRES[Length_of_mem-1];               save_P_PRES=P_PRES[Length_of_mem-1];
399                 for (i=Length_of_mem-1;i>0;--i) {               #pragma ivdep
400                 for (i=Length_of_mem-1;i>0;--i) {
401                     BREAKF[i]=BREAKF[i-1];                     BREAKF[i]=BREAKF[i-1];
402                     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];
403                     R_PRES[i]=R_PRES[i-1];                     R_PRES[i]=R_PRES[i-1];
404                     X_PRES[i]=X_PRES[i-1];                     X_PRES[i]=X_PRES[i-1];
405                     P_PRES[i]=P_PRES[i-1];                     P_PRES[i]=P_PRES[i-1];
406                 }               }
407                 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES;               R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES;
408                 R_PRES[0]=save_R_PRES;               R_PRES[0]=save_R_PRES;
409                 X_PRES[0]=save_XPRES;               X_PRES[0]=save_XPRES;
410                 P_PRES[0]=save_P_PRES;               P_PRES[0]=save_P_PRES;
411    
412                 if (ABS(Factor)<=ZERO) {               if (ABS(Factor)<=ZERO) {
413                    Factor=1.;                    Factor=1.;
414                    BREAKF[0]=ZERO;                    BREAKF[0]=ZERO;
415                 } else {               } else {
416                    Factor=1./Factor;                    Factor=1./Factor;
417                    BREAKF[0]=ONE;                    BREAKF[0]=ONE;
                }  
                for (i=0;i<order;++i) ALPHA[i]*=Factor;  
418              }              }
419                for (i=0;i<order;++i) ALPHA[i]*=Factor;
420              /*                                                                              /*                                                                
421              ***** update of solution X_PRES and its residual R_PRES:                          ***** update of solution X_PRES and its residual R_PRES:            
422              ***                                                                            ***                                                              
# Line 406  err_t Paso_Solver_GMRES( Line 425  err_t Paso_Solver_GMRES(
425              ***   R_PRES and X_PRES                                                  ***   R_PRES and X_PRES                                    
426              ***                                                                              ***                                                                
427              **/              **/
             #pragma omp barrier  
428              breakf0=BREAKF[0];              breakf0=BREAKF[0];
429              if (order==0) {              #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end)
430                #pragma omp for private(z) schedule(static)              for (th=0;th<num_threads;++th) {
431                for (z=0;z<n;++z) {                 local_n=n/num_threads;
432                    R_PRES[0][z]= Factor*       AP[z];                 rest=n-local_n*num_threads;
433                    X_PRES[0][z]=-Factor*P_PRES[1][z];                 n_start=local_n*th+MIN(th,rest);
434                }                 n_end=local_n*(th+1)+MIN(th+1,rest);
435              } else if (order==1) {                 if (order==0) {
436                #pragma omp for private(z) schedule(static)                     #pragma ivdep
437                for (z=0;z<n;++z) {                     for (z=n_start; z < n_end; ++z) {
438                    R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z];                        R_PRES[0][z]= Factor*       AP[z];
439                    X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];                        X_PRES[0][z]=-Factor*P_PRES[1][z];
440                }                     }
441              } else if (order==2) {                 } else if (order==1) {
442                #pragma omp for private(z) schedule(static)                     #pragma ivdep
443                for (z=0;z<n;++z) {                     for (z=n_start; z < n_end; ++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];                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z];
446                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     }
447                                                    +ALPHA[1]*X_PRES[2][z];                 } else if (order==2) {
448                }                     #pragma ivdep
449              } else if (order==3) {                     for (z=n_start; z < n_end; ++z) {
450                #pragma omp for private(z) schedule(static)                         R_PRES[0][z]= Factor*       AP[z]+ALPHA[0]*R_PRES[1][z]
451                for (z=0;z<n;++z) {                                                          +ALPHA[1]*R_PRES[2][z];
452                   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]
453                                             +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[1]*X_PRES[2][z];
454                                             +ALPHA[2]*R_PRES[3][z];                     }
455                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                 } else if (order==3) {
456                                                    +ALPHA[1]*X_PRES[2][z]                     #pragma ivdep
457                                                    +ALPHA[2]*X_PRES[3][z];                     for (z=n_start; z < n_end; ++z) {
458                }                         R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
459              } else if (order==4) {                                                   +ALPHA[1]*R_PRES[2][z]
460                #pragma omp for private(z) schedule(static)                                                   +ALPHA[2]*R_PRES[3][z];
461                for (z=0;z<n;++z) {                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
462                   R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[1]*X_PRES[2][z]
463                                             +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[2]*X_PRES[3][z];
464                                             +ALPHA[2]*R_PRES[3][z]                     }
465                                             +ALPHA[3]*R_PRES[4][z];                 } else if (order==4) {
466                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     #pragma ivdep
467                                                    +ALPHA[1]*X_PRES[2][z]                     for (z=n_start; z < n_end; ++z) {
468                                                    +ALPHA[2]*X_PRES[3][z]                         R_PRES[0][z]= Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
469                                                    +ALPHA[3]*X_PRES[4][z];                                                   +ALPHA[1]*R_PRES[2][z]
470                }                                                   +ALPHA[2]*R_PRES[3][z]
471              } else if (order==5) {                                                   +ALPHA[3]*R_PRES[4][z];
472                #pragma omp for private(z) schedule(static)                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
473                for (z=0;z<n;++z) {                                                          +ALPHA[1]*X_PRES[2][z]
474                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[2]*X_PRES[3][z]
475                                            +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[3]*X_PRES[4][z];
476                                            +ALPHA[2]*R_PRES[3][z]                     }
477                                            +ALPHA[3]*R_PRES[4][z]                 } else if (order==5) {
478                                            +ALPHA[4]*R_PRES[5][z];                     #pragma ivdep
479                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                     for (z=n_start; z < n_end; ++z) {
480                                                    +ALPHA[1]*X_PRES[2][z]                         R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
481                                                    +ALPHA[2]*X_PRES[3][z]                                                  +ALPHA[1]*R_PRES[2][z]
482                                                    +ALPHA[3]*X_PRES[4][z]                                                  +ALPHA[2]*R_PRES[3][z]
483                                                    +ALPHA[4]*X_PRES[5][z];                                                  +ALPHA[3]*R_PRES[4][z]
484                }                                                  +ALPHA[4]*R_PRES[5][z];
485              } else if (order==6) {                         X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
486                #pragma omp for private(z) schedule(static)                                                          +ALPHA[1]*X_PRES[2][z]
487                for (z=0;z<n;++z) {                                                          +ALPHA[2]*X_PRES[3][z]
488                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                          +ALPHA[3]*X_PRES[4][z]
489                                            +ALPHA[1]*R_PRES[2][z]                                                          +ALPHA[4]*X_PRES[5][z];
490                                            +ALPHA[2]*R_PRES[3][z]                      }
491                                            +ALPHA[3]*R_PRES[4][z]                 } else if (order==6) {
492                                            +ALPHA[4]*R_PRES[5][z]                     #pragma ivdep
493                                            +ALPHA[5]*R_PRES[6][z];                     for (z=n_start; z < n_end; ++z) {
494                   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]
495                                                    +ALPHA[1]*X_PRES[2][z]                                                 +ALPHA[1]*R_PRES[2][z]
496                                                    +ALPHA[2]*X_PRES[3][z]                                                 +ALPHA[2]*R_PRES[3][z]
497                                                    +ALPHA[3]*X_PRES[4][z]                                                 +ALPHA[3]*R_PRES[4][z]
498                                                    +ALPHA[4]*X_PRES[5][z]                                                 +ALPHA[4]*R_PRES[5][z]
499                                                    +ALPHA[5]*X_PRES[6][z];                                                 +ALPHA[5]*R_PRES[6][z];
500                }                        X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
501              } else if (order>6) {                                                         +ALPHA[1]*X_PRES[2][z]
502                #pragma omp for private(z) schedule(static)                                                         +ALPHA[2]*X_PRES[3][z]
503                for (z=0;z<n;++z) {                                                         +ALPHA[3]*X_PRES[4][z]
504                   R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]                                                         +ALPHA[4]*X_PRES[5][z]
505                                            +ALPHA[1]*R_PRES[2][z]                                                         +ALPHA[5]*X_PRES[6][z];
506                                            +ALPHA[2]*R_PRES[3][z]                     }
507                                            +ALPHA[3]*R_PRES[4][z]                 } else {
508                                            +ALPHA[4]*R_PRES[5][z]                     #pragma ivdep
509                                            +ALPHA[5]*R_PRES[6][z]                     for (z=n_start; z < n_end; ++z) {
510                                            +ALPHA[6]*R_PRES[7][z];                          R_PRES[0][z]=Factor*AP[z]+ALPHA[0]*R_PRES[1][z]
511                   X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]                                                   +ALPHA[1]*R_PRES[2][z]
512                                                    +ALPHA[1]*X_PRES[2][z]                                                   +ALPHA[2]*R_PRES[3][z]
513                                                    +ALPHA[2]*X_PRES[3][z]                                                   +ALPHA[3]*R_PRES[4][z]
514                                                    +ALPHA[3]*X_PRES[4][z]                                                   +ALPHA[4]*R_PRES[5][z]
515                                                    +ALPHA[4]*X_PRES[5][z]                                                   +ALPHA[5]*R_PRES[6][z]
516                                                    +ALPHA[5]*X_PRES[6][z]                                                   +ALPHA[6]*R_PRES[7][z];
517                                                    +ALPHA[6]*X_PRES[7][z];                          X_PRES[0][z]=-Factor*P_PRES[1][z]+ALPHA[0]*X_PRES[1][z]
518                }                                                           +ALPHA[1]*X_PRES[2][z]
519                for (i=7;i<order;++i) {                                                           +ALPHA[2]*X_PRES[3][z]
520                  #pragma omp for private(z) schedule(static)                                                           +ALPHA[3]*X_PRES[4][z]
521                  for (z=0;z<n;++z) {                                                           +ALPHA[4]*X_PRES[5][z]
522                      R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];                                                           +ALPHA[5]*X_PRES[6][z]
523                      X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];                                                           +ALPHA[6]*X_PRES[7][z];
524                  }                     }
525                }                     for (i=7;i<order;++i) {
526                          #pragma ivdep
527                          for (z=n_start; z < n_end; ++z) {
528                             R_PRES[0][z]+=ALPHA[i]*R_PRES[i+1][z];
529                             X_PRES[0][z]+=ALPHA[i]*X_PRES[i+1][z];
530                         }
531                       }
532                   }
533              }              }
534              if (breakf0>0.) {              if (breakf0>0.) {
535                /***                /***
536                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:                ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}:
537                ***/                ***/
538                #pragma omp master                memset(loc_dots,0,sizeof(double)*3);
539                {                #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end,diff,SC1,SC2)
540                   SC1=ZERO;                for (th=0;th<num_threads;++th) {
541                   SC2=ZERO;                    local_n=n/num_threads;
542                   L2_R=ZERO;                    rest=n-local_n*num_threads;
543                }                    n_start=local_n*th+MIN(th,rest);
544                #pragma omp barrier                    n_end=local_n*(th+1)+MIN(th+1,rest);
545                #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static)                    SC1=ZERO;
546                for (z=0;z<n;++z) {                    SC2=ZERO;
547                  diff=R_PRES[0][z]-r[z];                    #pragma ivdep
548                  SC1+=diff*diff;                    for (z=n_start; z < n_end; ++z) {
549                  SC2+=diff*r[z];                         diff=R_PRES[0][z]-r[z];
550                }                         SC1+=diff*diff;
551                           SC2+=diff*r[z];
552                      }
553                      #pragma omp critical
554                      {
555                           loc_dots[0]+=SC1;
556                           loc_dots[1]+=SC2;
557                      }
558                  }
559                  #ifdef PASO_MPI
560                      MPI_Allreduce(loc_dots, dots, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
561                      SC1=dots[0];
562                      SC2=dots[1];
563                  #else
564                      SC1=loc_dots[0];
565                      SC2=loc_dots[1];
566                  #endif
567                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;
568                #pragma omp 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)
569                for (z=0;z<n;++z) {                for (th=0;th<num_threads;++th) {
570                   x[z]+=gamma*(X_PRES[0][z]-x[z]);                    local_n=n/num_threads;
571                   r[z]+=gamma*(R_PRES[0][z]-r[z]);                    rest=n-local_n*num_threads;
572                   L2_R+=r[z]*r[z];                    n_start=local_n*th+MIN(th,rest);
573                }                    n_end=local_n*(th+1)+MIN(th+1,rest);
574                      L2_R=ZERO;
575                      #pragma ivdep
576                      for (z=n_start; z < n_end; ++z) {
577                          x[z]+=gamma*(X_PRES[0][z]-x[z]);
578                          r[z]+=gamma*(R_PRES[0][z]-r[z]);
579                          L2_R+=r[z]*r[z];
580                      }
581                      #pragma omp critical
582                      {
583                           loc_dots[2]+=L2_R;
584                      }
585                  }
586                  #ifdef PASO_MPI
587                      MPI_Allreduce(&loc_dots[2], &dots[2], 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
588                      L2_R=dots[2];
589                  #else
590                      L2_R=loc_dots[2];
591                  #endif
592                norm_of_residual=sqrt(L2_R);                norm_of_residual=sqrt(L2_R);
593                convergeFlag = (norm_of_residual <= tol);                convergeFlag = (norm_of_residual <= tol);
594                if (restart>0) restartFlag=(num_iter_restart >= restart);                if (restart>0) restartFlag=(num_iter_restart >= restart);
# Line 541  err_t Paso_Solver_GMRES( Line 600  err_t Paso_Solver_GMRES(
600           }           }
601          }          }
602          /* end of iteration */          /* end of iteration */
603          #pragma omp master          Norm_of_residual_global=norm_of_residual;
604          {      Num_iter_global=num_iter;
605             Norm_of_residual_global=norm_of_residual;      if (maxIterFlag) {
        Num_iter_global=num_iter;  
        if (maxIterFlag) {  
606             Status = SOLVER_MAXITER_REACHED;             Status = SOLVER_MAXITER_REACHED;
607         } else if (breakFlag) {         } else if (breakFlag) {
608             Status = SOLVER_BREAKDOWN;             Status = SOLVER_BREAKDOWN;
609        }      }
610          }      }
     }  /* end of parallel region */  
     TMPMEMFREE(AP);  
611      for (i=0;i<Length_of_recursion;i++) {      for (i=0;i<Length_of_recursion;i++) {
612         TMPMEMFREE(X_PRES[i]);         TMPMEMFREE(X_PRES[i]);
613         TMPMEMFREE(R_PRES[i]);         TMPMEMFREE(R_PRES[i]);
614         TMPMEMFREE(P_PRES[i]);         TMPMEMFREE(P_PRES[i]);
615      }      }
616        TMPMEMFREE(AP);
617        TMPMEMFREE(X_PRES);
618        TMPMEMFREE(R_PRES);
619        TMPMEMFREE(P_PRES);
620        TMPMEMFREE(P_PRES_dot_AP);
621        TMPMEMFREE(R_PRES_dot_P_PRES);
622        TMPMEMFREE(BREAKF);
623        TMPMEMFREE(ALPHA);
624        TMPMEMFREE(dots);
625        TMPMEMFREE(loc_dots);
626      *iter=Num_iter_global;      *iter=Num_iter_global;
627      *tolerance=Norm_of_residual_global;      *tolerance=Norm_of_residual_global;
   }  
628    return Status;    return Status;
629  }  }
   

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

  ViewVC Help
Powered by ViewVC 1.1.26