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

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

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

revision 1841 by artak, Wed Oct 1 05:56:05 2008 UTC revision 1842 by artak, Fri Oct 3 04:18:26 2008 UTC
# Line 37  void Paso_Preconditioner_free(Paso_Solve Line 37  void Paso_Preconditioner_free(Paso_Solve
37        Paso_Solver_RILU_free(in->rilu);        Paso_Solver_RILU_free(in->rilu);
38        Paso_Solver_Jacobi_free(in->jacobi);        Paso_Solver_Jacobi_free(in->jacobi);
39        Paso_Solver_GS_free(in->gs);        Paso_Solver_GS_free(in->gs);
40          Paso_Solver_AMG_free(in->amg);
41        MEMFREE(in);        MEMFREE(in);
42      }      }
43  }  }
# Line 53  void Paso_Solver_setPreconditioner(Paso_ Line 54  void Paso_Solver_setPreconditioner(Paso_
54          prec->ilu=NULL;          prec->ilu=NULL;
55          prec->jacobi=NULL;          prec->jacobi=NULL;
56          prec->gs=NULL;          prec->gs=NULL;
57            prec->amg=NULL;
58          A->solver=prec;          A->solver=prec;
59          switch (options->preconditioner) {          switch (options->preconditioner) {
60             default:             default:
# Line 77  void Paso_Solver_setPreconditioner(Paso_ Line 79  void Paso_Solver_setPreconditioner(Paso_
79                prec->gs->sweeps=options->sweeps;                prec->gs->sweeps=options->sweeps;
80                prec->type=PASO_GS;                prec->type=PASO_GS;
81                break;                break;
82                case PASO_AMG:
83                  if (options->verbose) printf("AMG preconditioner is used.\n");
84                  prec->amg=Paso_Solver_getAMG(A->mainBlock,options->verbose);
85                  prec->type=PASO_AMG;
86                  break;
87    
88          }          }
89          if (! Paso_MPIInfo_noError(A->mpi_info ) ){          if (! Paso_MPIInfo_noError(A->mpi_info ) ){
90             Paso_Preconditioner_free(prec);             Paso_Preconditioner_free(prec);
# Line 90  void Paso_Solver_setPreconditioner(Paso_ Line 98  void Paso_Solver_setPreconditioner(Paso_
98  /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */  /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
99  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
100      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
101        
102        
103      switch (prec->type) {      switch (prec->type) {
104          default:          default:
105          case PASO_JACOBI:          case PASO_JACOBI:
# Line 102  void Paso_Solver_solvePreconditioner(Pas Line 112  void Paso_Solver_solvePreconditioner(Pas
112             Paso_Solver_solveRILU(prec->rilu,x,b);             Paso_Solver_solveRILU(prec->rilu,x,b);
113             break;             break;
114           case PASO_GS:           case PASO_GS:
115                /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
116                  We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
117                  x_0=0
118                  x_1=P^{-1}(b)
119                  
120                  b_old=b
121                  
122                  b_new=b_old+b
123                  x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
124                     =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
125                  b_old=b_new
126                  
127                  b_new=b_old+b
128                  x_3=....=P^{-1}(b_new-AP^{-1}b_old)
129                  b_old=b_new
130                  
131                  So for n- steps we will use loop, but we always make sure that every value calculated only once!
132                */
133    
134             Paso_Solver_solveGS(prec->gs,x,b);             Paso_Solver_solveGS(prec->gs,x,b);
135               if (prec->gs->sweeps>1) {
136               double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
137               double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
138               dim_t i;
139               #pragma omp parallel for private(i) schedule(static)
140               for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
141              
142               while(prec->gs->sweeps>1) {
143                   #pragma omp parallel for private(i) schedule(static)
144                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
145                    /* Compute the residual b=b-Ax*/
146                   Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
147                   /* Go round again*/
148                   Paso_Solver_solveGS(prec->gs,x,bnew);
149                   #pragma omp parallel for private(i) schedule(static)
150                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
151                   prec->gs->sweeps=prec->gs->sweeps-1;
152               }
153              
154               MEMFREE(bold);
155               MEMFREE(bnew);
156              
157               }
158               /* prec->gs->sweeps=prec->gs->sweeps-1;*/
159              
160               break;
161             case PASO_AMG:
162               Paso_Solver_solveAMG(prec->amg,x,b);
163             break;             break;
164      }      }
165    
166  }  }

Legend:
Removed from v.1841  
changed lines
  Added in v.1842

  ViewVC Help
Powered by ViewVC 1.1.26