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

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

  ViewVC Help
Powered by ViewVC 1.1.26