/[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 temp_trunk_copy/paso/src/GMRES.c revision 1384 by phornby, Fri Jan 11 02:29:38 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 136  err_t Paso_Solver_GMRES( Line 160  err_t Paso_Solver_GMRES(
160           /***                                                                           /***                                                                
161           *** calculate new search direction P from R_PRES           *** calculate new search direction P from R_PRES
162           ***/           ***/
163             #pragma omp barrier
164           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);           Paso_Solver_solvePreconditioner(A,&P_PRES[0][0], &R_PRES[0][0]);
165           /***                                                                           /***                                                                
166           *** apply A to P to get AP           *** apply A to P to get AP
167           ***/           ***/
168             #pragma omp barrier
169       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]);
170           /***                                                                           /***                                                                
171           ***** 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 378  err_t Paso_Solver_GMRES(
378           /* this fixes a problem with the intel compiler */           /* this fixes a problem with the intel compiler */
379           #pragma omp master           #pragma omp master
380           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];           P_PRES_dot_AP0=R_PRES_dot_P_PRES[0];
          #pragma omp barrier  
381           /***   if sum_BREAKF is equal to zero a breakdown occurs.           /***   if sum_BREAKF is equal to zero a breakdown occurs.
382            ***   iteration procedure can be continued but R_PRES is not the            ***   iteration procedure can be continued but R_PRES is not the
383            ***   residual of X_PRES approximation.            ***   residual of X_PRES approximation.
384            ***/            ***/
385             #pragma omp barrier
386           sum_BREAKF=0.;           sum_BREAKF=0.;
387           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];           for (i=0;i<order;++i) sum_BREAKF +=BREAKF[i];
388           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 578  err_t Paso_Solver_GMRES(
578        }        }
579          }          }
580      }  /* end of parallel region */      }  /* end of parallel region */
     TMPMEMFREE(AP);  
581      for (i=0;i<Length_of_recursion;i++) {      for (i=0;i<Length_of_recursion;i++) {
582         TMPMEMFREE(X_PRES[i]);         TMPMEMFREE(X_PRES[i]);
583         TMPMEMFREE(R_PRES[i]);         TMPMEMFREE(R_PRES[i]);
584         TMPMEMFREE(P_PRES[i]);         TMPMEMFREE(P_PRES[i]);
585      }      }
586        TMPMEMFREE(AP);
587        TMPMEMFREE(X_PRES);
588        TMPMEMFREE(R_PRES);
589        TMPMEMFREE(P_PRES);
590        TMPMEMFREE(P_PRES_dot_AP);
591        TMPMEMFREE(R_PRES_dot_P_PRES);
592        TMPMEMFREE(BREAKF);
593        TMPMEMFREE(ALPHA);
594      *iter=Num_iter_global;      *iter=Num_iter_global;
595      *tolerance=Norm_of_residual_global;      *tolerance=Norm_of_residual_global;
596    }    }
597    return Status;    return Status;
598  }  }
   

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

  ViewVC Help
Powered by ViewVC 1.1.26