/[escript]/trunk/paso/src/Solver_preconditioner.c
ViewVC logotype

Diff of /trunk/paso/src/Solver_preconditioner.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3093 by lgao, Tue Jun 29 00:45:49 2010 UTC revision 3094 by gross, Fri Aug 13 08:38:06 2010 UTC
# Line 37  void Paso_Preconditioner_free(Paso_Solve Line 37  void Paso_Preconditioner_free(Paso_Solve
37        Paso_Solver_ILU_free(in->ilu);        Paso_Solver_ILU_free(in->ilu);
38        Paso_Solver_RILU_free(in->rilu);        Paso_Solver_RILU_free(in->rilu);
39        Paso_Solver_Jacobi_free(in->jacobi);        Paso_Solver_Jacobi_free(in->jacobi);
40        if (in->type==PASO_GAUSS_SEIDEL_MPI)        Paso_Solver_GS_free(in->gs);
         Paso_Solver_GSMPI_free(in->gs);  
       else  
         Paso_Solver_GS_free(in->gs);  
41        Paso_Solver_AMG_free(in->amg);        Paso_Solver_AMG_free(in->amg);
42        Paso_Solver_AMG_System_free(in->amgSystem);        Paso_Solver_AMG_System_free(in->amgSystem);
43        Paso_Solver_AMLI_free(in->amli);        Paso_Solver_AMLI_free(in->amli);
# Line 72  void Paso_Solver_setPreconditioner(Paso_ Line 69  void Paso_Solver_setPreconditioner(Paso_
69          prec->amli=NULL;          prec->amli=NULL;
70          prec->amgSystem=NULL;          prec->amgSystem=NULL;
71          prec->amliSystem=NULL;          prec->amliSystem=NULL;
72        if (options->verbose && options->use_local_preconditioner) printf("Apply preconditioner locally only.\n");
73    
74          A->solver=prec;          A->solver=prec;
75          switch (options->preconditioner) {          switch (options->preconditioner) {
76             default:             default:
77             case PASO_JACOBI:             case PASO_JACOBI:
78                if (options->verbose) printf("Jacobi preconditioner is used.\n");                if (options->verbose) printf("Jacobi preconditioner is used.\n");
79                prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);                prec->jacobi=Paso_Solver_getJacobi(A);
80                prec->type=PASO_JACOBI;                prec->type=PASO_JACOBI;
81                break;                break;
82           case PASO_GS:
83              if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
84              prec->gs=Paso_Solver_getGS(A, options->sweeps, options->use_local_preconditioner, options->verbose);
85              prec->type=PASO_GS;
86              break;
87              
88           /***************************************************************************************/  
89             case PASO_ILU0:             case PASO_ILU0:
90                if (options->verbose) printf("ILU preconditioner is used.\n");                if (options->verbose) printf("ILU preconditioner is used.\n");
91                prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);                prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
92                prec->type=PASO_ILU0;                prec->type=PASO_ILU0;
93              Paso_MPIInfo_noError(A->mpi_info);
94                break;                break;
95             case PASO_RILU:             case PASO_RILU:
96                if (options->verbose) printf("RILU preconditioner is used.\n");                if (options->verbose) printf("RILU preconditioner is used.\n");
97                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
98              Paso_MPIInfo_noError(A->mpi_info);
99                prec->type=PASO_RILU;                prec->type=PASO_RILU;
100                break;                break;
101              case PASO_GS:  
               if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");  
               prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);  
               prec->gs->sweeps=options->sweeps;  
               prec->type=PASO_GS;  
               break;  
             case PASO_GAUSS_SEIDEL_MPI:  
               if (options->verbose) printf("MPI versioned Gauss-Seidel preconditioner is used.\n");  
               prec->gs=Paso_Solver_getGSMPI(A->mainBlock,options->verbose);  
               prec->gs->sweeps=options->sweeps;  
               prec->type=PASO_GAUSS_SEIDEL_MPI;  
               break;  
102              case PASO_AMG:              case PASO_AMG:
103                if (options->verbose) printf("AMG preconditioner is used.\n");                if (options->verbose) printf("AMG preconditioner is used.\n");
104                prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);                  prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
105              Paso_MPIInfo_noError(A->mpi_info);
106            /*            /*
107                prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);                prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
108            if (Paso_checkPtr(prec->amgSystem)) return;            if (Paso_checkPtr(prec->amgSystem)) return;
# Line 173  void Paso_Solver_solvePreconditioner(Pas Line 170  void Paso_Solver_solvePreconditioner(Pas
170      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
171      dim_t i,j;      dim_t i,j;
172      dim_t n=A->mainBlock->numRows;      dim_t n=A->mainBlock->numRows;
173      double **xx;  
     double **bb;  
174            
175      xx=MEMALLOC(A->row_block_size,double*);  
     bb=MEMALLOC(A->row_block_size,double*);  
     if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;  
176    
177      switch (prec->type) {      switch (prec->type) {
178          default:          default:
179          case PASO_JACOBI:          case PASO_JACOBI:
180             Paso_Solver_solveJacobi(prec->jacobi,x,b);             Paso_Solver_solveJacobi(A, prec->jacobi,x,b);
181             break;             break;
182          
183        case PASO_GS:
184           Paso_Solver_solveGS(A, prec->gs,x,b);
185           break;      
186          case PASO_ILU0:          case PASO_ILU0:
187             Paso_Solver_solveILU(prec->ilu,x,b);             Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
188             break;             break;
189          case PASO_RILU:          case PASO_RILU:
190             Paso_Solver_solveRILU(prec->rilu,x,b);         Paso_Solver_solveRILU(prec->rilu,x,b);
            break;  
         case PASO_GS:  
             /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.  
               We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:  
               x_0=0  
               x_1=P^{-1}(b)  
                 
               b_old=b  
                 
               b_new=b_old+b  
               x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))  
                  =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)  
               b_old=b_new  
                 
               b_new=b_old+b  
               b_new=b_new-Ax_2=b_new-AP^{-1}(2b-AP^{-1}b)  
               x_3=....=P^{-1}(b_new-AP^{-1}b_old)  
               b_old=b_new  
                 
               So for n- steps we will use loop, but we always make sure that every value calculated only once!  
             */  
               
            Paso_Solver_solveGS(prec->gs,x,b);  
            if (prec->gs->sweeps>1) {  
                 double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);  
                 dim_t sweeps=prec->gs->sweeps;  
             
                 #pragma omp parallel for private(i) schedule(static)  
                 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=b[i];  
             
                 while(sweeps>1) {  
                    #pragma omp parallel for private(i) schedule(static)  
                    for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]+=b[i];  
                     /* Compute the residual b=b-Ax*/  
                    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A->mainBlock, x, DBLE(1), bnew);  
                    /* Go round again*/  
                    Paso_Solver_solveGS(prec->gs,x,bnew);  
                    sweeps--;  
                 }  
                 MEMFREE(bnew);  
            }  
            break;  
         case PASO_GAUSS_SEIDEL_MPI:  
            Paso_Solver_solveGSMPI(A,prec->gs,x,b);  
191             break;             break;
192    
193          case PASO_AMG:          case PASO_AMG:
194              Paso_Solver_solveAMG(prec->amg,x,b);              Paso_Solver_solveAMG(prec->amg,x,b);
195    
# Line 283  void Paso_Solver_solvePreconditioner(Pas Line 238  void Paso_Solver_solvePreconditioner(Pas
238                  Paso_Solver_solveAMLI(prec->amli,x,b);                  Paso_Solver_solveAMLI(prec->amli,x,b);
239              }              }
240              else {              else {
241                         double **xx;
242                               double **bb;
243               xx=MEMALLOC(A->row_block_size,double*);
244               bb=MEMALLOC(A->row_block_size,double*);
245               if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
246                   for (i=0;i<A->row_block_size;i++) {                   for (i=0;i<A->row_block_size;i++) {
247                      xx[i]=MEMALLOC(n,double);                      xx[i]=MEMALLOC(n,double);
248                      bb[i]=MEMALLOC(n,double);                      bb[i]=MEMALLOC(n,double);
# Line 314  void Paso_Solver_solvePreconditioner(Pas Line 272  void Paso_Solver_solvePreconditioner(Pas
272                  MEMFREE(xx[i]);                  MEMFREE(xx[i]);
273                  MEMFREE(bb[i]);                  MEMFREE(bb[i]);
274                  }                  }
275                                  MEMFREE(xx);
276            MEMFREE(bb);
277              }              }
278          break;          break;
279      }      }
280      MEMFREE(xx);  
     MEMFREE(bb);  
281  }  }

Legend:
Removed from v.3093  
changed lines
  Added in v.3094

  ViewVC Help
Powered by ViewVC 1.1.26