/[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 1973 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 1974 by jfenwick, Thu Nov 6 02:40:10 2008 UTC
# Line 148  err_t Paso_Solver_GMRES( Line 148  err_t Paso_Solver_GMRES(
148    
149        while (! (convergeFlag || maxIterFlag || breakFlag))  {        while (! (convergeFlag || maxIterFlag || breakFlag))  {
150           if (restartFlag) {           if (restartFlag) {
151               BREAKF[0]=ONE;               BREAKF[0]=PASO_ONE;
152               #pragma omp parallel for private(th,z, local_n, rest, n_start, n_end)               #pragma omp parallel for private(th,z, local_n, rest, n_start, n_end)
153               for (th=0;th<num_threads;++th) {               for (th=0;th<num_threads;++th) {
154                 local_n=n/num_threads;                 local_n=n/num_threads;
# Line 172  err_t Paso_Solver_GMRES( Line 172  err_t Paso_Solver_GMRES(
172           /***                                                                           /***                                                                
173           *** apply A to P to get AP           *** apply A to P to get AP
174           ***/           ***/
175       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, &P_PRES[0][0],PASO_ZERO, &AP[0]);
176           /***                                                                           /***                                                                
177           ***** calculation of the norm of R and the scalar products of                 ***** calculation of the norm of R and the scalar products of      
178           ***   the residuals and A*P:                                                   ***   the residuals and A*P:                                        
# Line 185  err_t Paso_Solver_GMRES( Line 185  err_t Paso_Solver_GMRES(
185                 n_start=local_n*th+MIN(th,rest);                 n_start=local_n*th+MIN(th,rest);
186                 n_end=local_n*(th+1)+MIN(th+1,rest);                 n_end=local_n*(th+1)+MIN(th+1,rest);
187                 if (order==0) {                 if (order==0) {
188                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
189                     #pragma ivdep                     #pragma ivdep
190                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
191                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 195  err_t Paso_Solver_GMRES( Line 195  err_t Paso_Solver_GMRES(
195                          loc_dots[0]+=R_PRES_dot_P;                          loc_dots[0]+=R_PRES_dot_P;
196                     }                     }
197                 } else if (order==1) {                 } else if (order==1) {
198                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
199                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
200                     #pragma ivdep                     #pragma ivdep
201                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
202                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 208  err_t Paso_Solver_GMRES( Line 208  err_t Paso_Solver_GMRES(
208                          loc_dots[1]+=P_PRES_dot_AP0;                          loc_dots[1]+=P_PRES_dot_AP0;
209                     }                     }
210                 } else if (order==2) {                 } else if (order==2) {
211                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
212                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
213                     P_PRES_dot_AP1=ZERO;                     P_PRES_dot_AP1=PASO_ZERO;
214                     #pragma ivdep                     #pragma ivdep
215                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
216                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 224  err_t Paso_Solver_GMRES( Line 224  err_t Paso_Solver_GMRES(
224                          loc_dots[2]+=P_PRES_dot_AP1;                          loc_dots[2]+=P_PRES_dot_AP1;
225                     }                     }
226                 } else if (order==3) {                 } else if (order==3) {
227                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
228                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
229                     P_PRES_dot_AP1=ZERO;                     P_PRES_dot_AP1=PASO_ZERO;
230                     P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP2=PASO_ZERO;
231                     #pragma ivdep                     #pragma ivdep
232                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
233                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 243  err_t Paso_Solver_GMRES( Line 243  err_t Paso_Solver_GMRES(
243                          loc_dots[3]+=P_PRES_dot_AP2;                          loc_dots[3]+=P_PRES_dot_AP2;
244                     }                     }
245                 } else if (order==4) {                 } else if (order==4) {
246                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
247                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
248                     P_PRES_dot_AP1=ZERO;                     P_PRES_dot_AP1=PASO_ZERO;
249                     P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP2=PASO_ZERO;
250                     P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP3=PASO_ZERO;
251                     #pragma ivdep                     #pragma ivdep
252                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
253                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 265  err_t Paso_Solver_GMRES( Line 265  err_t Paso_Solver_GMRES(
265                          loc_dots[4]+=P_PRES_dot_AP3;                          loc_dots[4]+=P_PRES_dot_AP3;
266                     }                     }
267                 } else if (order==5) {                 } else if (order==5) {
268                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
269                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
270                     P_PRES_dot_AP1=ZERO;                     P_PRES_dot_AP1=PASO_ZERO;
271                     P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP2=PASO_ZERO;
272                     P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP3=PASO_ZERO;
273                     P_PRES_dot_AP4=ZERO;                     P_PRES_dot_AP4=PASO_ZERO;
274                     #pragma ivdep                     #pragma ivdep
275                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
276                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 290  err_t Paso_Solver_GMRES( Line 290  err_t Paso_Solver_GMRES(
290                          loc_dots[5]+=P_PRES_dot_AP4;                          loc_dots[5]+=P_PRES_dot_AP4;
291                     }                     }
292                 } else if (order==6) {                 } else if (order==6) {
293                     R_PRES_dot_P=ZERO;                     R_PRES_dot_P=PASO_ZERO;
294                     P_PRES_dot_AP0=ZERO;                     P_PRES_dot_AP0=PASO_ZERO;
295                     P_PRES_dot_AP1=ZERO;                     P_PRES_dot_AP1=PASO_ZERO;
296                     P_PRES_dot_AP2=ZERO;                     P_PRES_dot_AP2=PASO_ZERO;
297                     P_PRES_dot_AP3=ZERO;                     P_PRES_dot_AP3=PASO_ZERO;
298                     P_PRES_dot_AP4=ZERO;                     P_PRES_dot_AP4=PASO_ZERO;
299                     P_PRES_dot_AP5=ZERO;                     P_PRES_dot_AP5=PASO_ZERO;
300                     #pragma ivdep                     #pragma ivdep
301                     for (z=n_start; z < n_end; ++z) {                     for (z=n_start; z < n_end; ++z) {
302                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                         R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 318  err_t Paso_Solver_GMRES( Line 318  err_t Paso_Solver_GMRES(
318                          loc_dots[6]+=P_PRES_dot_AP5;                          loc_dots[6]+=P_PRES_dot_AP5;
319                     }                     }
320                 } else {                 } else {
321                    R_PRES_dot_P=ZERO;                    R_PRES_dot_P=PASO_ZERO;
322                    P_PRES_dot_AP0=ZERO;                    P_PRES_dot_AP0=PASO_ZERO;
323                    P_PRES_dot_AP1=ZERO;                    P_PRES_dot_AP1=PASO_ZERO;
324                    P_PRES_dot_AP2=ZERO;                    P_PRES_dot_AP2=PASO_ZERO;
325                    P_PRES_dot_AP3=ZERO;                    P_PRES_dot_AP3=PASO_ZERO;
326                    P_PRES_dot_AP4=ZERO;                    P_PRES_dot_AP4=PASO_ZERO;
327                    P_PRES_dot_AP5=ZERO;                    P_PRES_dot_AP5=PASO_ZERO;
328                    P_PRES_dot_AP6=ZERO;                    P_PRES_dot_AP6=PASO_ZERO;
329                    #pragma ivdep                    #pragma ivdep
330                    for (z=n_start; z < n_end; ++z) {                    for (z=n_start; z < n_end; ++z) {
331                        R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];                        R_PRES_dot_P+=R_PRES[0][z]*P_PRES[0][z];
# Line 349  err_t Paso_Solver_GMRES( Line 349  err_t Paso_Solver_GMRES(
349                         loc_dots[7]+=P_PRES_dot_AP6;                         loc_dots[7]+=P_PRES_dot_AP6;
350                    }                    }
351                    for (i=7;i<order;++i) {                    for (i=7;i<order;++i) {
352                         P_PRES_dot_AP0=ZERO;                         P_PRES_dot_AP0=PASO_ZERO;
353                         #pragma ivdep                         #pragma ivdep
354                         for (z=n_start; z < n_end; ++z) {                         for (z=n_start; z < n_end; ++z) {
355                               P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];                               P_PRES_dot_AP0+=P_PRES[i][z]*AP[z];
# Line 378  err_t Paso_Solver_GMRES( Line 378  err_t Paso_Solver_GMRES(
378           sum_BREAKF=0.;           sum_BREAKF=0.;
379           #pragma ivdep           #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(R_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));           breakFlag=!((ABS(R_PRES_dot_AP0) > PASO_ZERO) &&  (sum_BREAKF >PASO_ZERO));
382           if (!breakFlag) {           if (!breakFlag) {
383              breakFlag=FALSE;              breakFlag=FALSE;
384              /***              /***
# Line 408  err_t Paso_Solver_GMRES( Line 408  err_t Paso_Solver_GMRES(
408               X_PRES[0]=save_XPRES;               X_PRES[0]=save_XPRES;
409               P_PRES[0]=save_P_PRES;               P_PRES[0]=save_P_PRES;
410    
411               if (ABS(Factor)<=ZERO) {               if (ABS(Factor)<=PASO_ZERO) {
412                    Factor=1.;                    Factor=1.;
413                    BREAKF[0]=ZERO;                    BREAKF[0]=PASO_ZERO;
414               } else {               } else {
415                    Factor=1./Factor;                    Factor=1./Factor;
416                    BREAKF[0]=ONE;                    BREAKF[0]=PASO_ONE;
417              }              }
418              for (i=0;i<order;++i) ALPHA[i]*=Factor;              for (i=0;i<order;++i) ALPHA[i]*=Factor;
419              /*                                                                              /*                                                                
# Line 541  err_t Paso_Solver_GMRES( Line 541  err_t Paso_Solver_GMRES(
541                    rest=n-local_n*num_threads;                    rest=n-local_n*num_threads;
542                    n_start=local_n*th+MIN(th,rest);                    n_start=local_n*th+MIN(th,rest);
543                    n_end=local_n*(th+1)+MIN(th+1,rest);                    n_end=local_n*(th+1)+MIN(th+1,rest);
544                    SC1=ZERO;                    SC1=PASO_ZERO;
545                    SC2=ZERO;                    SC2=PASO_ZERO;
546                    #pragma ivdep                    #pragma ivdep
547                    for (z=n_start; z < n_end; ++z) {                    for (z=n_start; z < n_end; ++z) {
548                         diff=R_PRES[0][z]-r[z];                         diff=R_PRES[0][z]-r[z];
# Line 563  err_t Paso_Solver_GMRES( Line 563  err_t Paso_Solver_GMRES(
563                    SC1=loc_dots[0];                    SC1=loc_dots[0];
564                    SC2=loc_dots[1];                    SC2=loc_dots[1];
565                #endif                #endif
566                gamma=(SC1<=ZERO) ? ZERO : -SC2/SC1;                gamma=(SC1<=PASO_ZERO) ? PASO_ZERO : -SC2/SC1;
567                #pragma omp parallel for private(th,z,local_n, rest, n_start, n_end,diff,L2_R)                #pragma omp parallel for private(th,z,local_n, rest, n_start, n_end,diff,L2_R)
568                for (th=0;th<num_threads;++th) {                for (th=0;th<num_threads;++th) {
569                    local_n=n/num_threads;                    local_n=n/num_threads;
570                    rest=n-local_n*num_threads;                    rest=n-local_n*num_threads;
571                    n_start=local_n*th+MIN(th,rest);                    n_start=local_n*th+MIN(th,rest);
572                    n_end=local_n*(th+1)+MIN(th+1,rest);                    n_end=local_n*(th+1)+MIN(th+1,rest);
573                    L2_R=ZERO;                    L2_R=PASO_ZERO;
574                    #pragma ivdep                    #pragma ivdep
575                    for (z=n_start; z < n_end; ++z) {                    for (z=n_start; z < n_end; ++z) {
576                        x[z]+=gamma*(X_PRES[0][z]-x[z]);                        x[z]+=gamma*(X_PRES[0][z]-x[z]);

Legend:
Removed from v.1973  
changed lines
  Added in v.1974

  ViewVC Help
Powered by ViewVC 1.1.26