/[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 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    
4    /*
5    ********************************************************************************
6    *               Copyright   2006 by ACcESS MNRF                                *
7    *                                                                              *
8    *                 http://www.access.edu.au                                     *
9    *           Primary Business: Queensland, Australia                            *
10    *     Licensed under the Open Software License version 3.0             *
11    *        http://www.opensource.org/licenses/osl-3.0.php                        *
12    ********************************************************************************
13    */
14    
15  /*  /*
16  *  Purpose  *  Purpose
17  *  =======  *  =======
# Line 56  err_t Paso_Solver_GMRES( Line 68  err_t Paso_Solver_GMRES(
68      double * r,      double * r,
69      double * x,      double * x,
70      dim_t *iter,      dim_t *iter,
71      double * tolerance,dim_t Length_of_recursion,dim_t restart) {      double * tolerance,dim_t Length_of_recursion,dim_t restart,
72        Paso_Performance* pp) {
73    
74    /* Local variables */    /* Local variables */
75    
76    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;
77    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;
78    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;
79    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;
80    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;
81    dim_t maxit,Num_iter_global,num_iter_restart,num_iter;    dim_t maxit,Num_iter_global,num_iter_restart,num_iter;
82    dim_t i,z,order;    dim_t i,z,order,n, Length_of_mem;
83    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE;
84    err_t Status=SOLVER_NO_ERROR;    err_t Status=SOLVER_NO_ERROR;
85    
86      
87    /* adapt original routine parameters */    /* adapt original routine parameters */
88    
89    dim_t n=A->num_cols * A-> col_block_size;    n=A->num_cols * A-> col_block_size;
90    dim_t Length_of_mem=MAX(Length_of_recursion,0)+1;    Length_of_mem=MAX(Length_of_recursion,0)+1;
91    
92    
93    /*     Test the input parameters. */    /*     Test the input parameters. */
94    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 97  err_t Paso_Solver_GMRES(
97    }    }
98    
99    /*     allocate memory: */    /*     allocate memory: */
100      X_PRES=TMPMEMALLOC(Length_of_mem,double*);
101      R_PRES=TMPMEMALLOC(Length_of_mem,double*);
102      P_PRES=TMPMEMALLOC(Length_of_mem,double*);
103      P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double);
104      R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double);
105      BREAKF=TMPMEMALLOC(Length_of_mem,double);
106      ALPHA=TMPMEMALLOC(Length_of_mem,double);
107    AP=TMPMEMALLOC(n,double);    AP=TMPMEMALLOC(n,double);
108    if (AP==NULL) Status=SOLVER_MEMORY_ERROR;    if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL ||
109    for (i=0;i<Length_of_mem;i++) {           P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL) {
110      X_PRES[i]=TMPMEMALLOC(n,double);       Status=SOLVER_MEMORY_ERROR;
111      R_PRES[i]=TMPMEMALLOC(n,double);    } else {
112      P_PRES[i]=TMPMEMALLOC(n,double);       for (i=0;i<Length_of_mem;i++) {
113      if (X_PRES[i]==NULL || R_PRES[i]==NULL ||  P_PRES[i]==NULL) Status=SOLVER_MEMORY_ERROR;         X_PRES[i]=TMPMEMALLOC(n,double);
114           R_PRES[i]=TMPMEMALLOC(n,double);
115           P_PRES[i]=TMPMEMALLOC(n,double);
116           if (X_PRES[i]==NULL || R_PRES[i]==NULL ||  P_PRES[i]==NULL) Status=SOLVER_MEMORY_ERROR;
117         }
118    }    }
119    if ( Status ==SOLVER_NO_ERROR ) {    if ( Status ==SOLVER_NO_ERROR ) {
120    
# Line 136  err_t Paso_Solver_GMRES( Line 161  err_t Paso_Solver_GMRES(
161           /***                                                                           /***                                                                
162           *** calculate new search direction P from R_PRES           *** calculate new search direction P from R_PRES
163           ***/           ***/
164             #pragma omp barrier
165           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);
166           /***                                                                           /***                                                                
167           *** apply A to P to get AP           *** apply A to P to get AP
168           ***/           ***/
169             #pragma omp barrier
170       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &P_PRES[0][0],ZERO, &AP[0]);
171           /***                                                                           /***                                                                
172           ***** calculation of the norm of R and the scalar products of                 ***** calculation of the norm of R and the scalar products of      
# Line 352  err_t Paso_Solver_GMRES( Line 379  err_t Paso_Solver_GMRES(
379           /* this fixes a problem with the intel compiler */           /* this fixes a problem with the intel compiler */
380           #pragma omp master           #pragma omp master
381           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];
          #pragma omp barrier  
382           /***   if sum_BREAKF is equal to zero a breakdown occurs.           /***   if sum_BREAKF is equal to zero a breakdown occurs.
383            ***   iteration procedure can be continued but R_PRES is not the            ***   iteration procedure can be continued but R_PRES is not the
384            ***   residual of X_PRES approximation.            ***   residual of X_PRES approximation.
385            ***/            ***/
386             #pragma omp barrier
387           sum_BREAKF=0.;           sum_BREAKF=0.;
388           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
389           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));           breakFlag=!((ABS(P_PRES_dot_AP0) > ZERO) &&  (sum_BREAKF >ZERO));
# Line 552  err_t Paso_Solver_GMRES( Line 579  err_t Paso_Solver_GMRES(
579        }        }
580          }          }
581      }  /* end of parallel region */      }  /* end of parallel region */
     TMPMEMFREE(AP);  
582      for (i=0;i<Length_of_recursion;i++) {      for (i=0;i<Length_of_recursion;i++) {
583         TMPMEMFREE(X_PRES[i]);         TMPMEMFREE(X_PRES[i]);
584         TMPMEMFREE(R_PRES[i]);         TMPMEMFREE(R_PRES[i]);
585         TMPMEMFREE(P_PRES[i]);         TMPMEMFREE(P_PRES[i]);
586      }      }
587        TMPMEMFREE(AP);
588        TMPMEMFREE(X_PRES);
589        TMPMEMFREE(R_PRES);
590        TMPMEMFREE(P_PRES);
591        TMPMEMFREE(P_PRES_dot_AP);
592        TMPMEMFREE(R_PRES_dot_P_PRES);
593        TMPMEMFREE(BREAKF);
594        TMPMEMFREE(ALPHA);
595      *iter=Num_iter_global;      *iter=Num_iter_global;
596      *tolerance=Norm_of_residual_global;      *tolerance=Norm_of_residual_global;
597    }    }
598    return Status;    return Status;
599  }  }
   

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

  ViewVC Help
Powered by ViewVC 1.1.26