/[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 2991 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 2992 by artak, Tue Mar 23 04:16:28 2010 UTC
# Line 50  void Paso_Solver_setPreconditioner(Paso_ Line 50  void Paso_Solver_setPreconditioner(Paso_
50    
51      Paso_Solver_Preconditioner* prec=NULL;      Paso_Solver_Preconditioner* prec=NULL;
52      dim_t i;      dim_t i;
53        /*char filename[7];*/
54      if (A->solver==NULL) {      if (A->solver==NULL) {
55          /* allocate structure to hold preconditioner */          /* allocate structure to hold preconditioner */
56          prec=MEMALLOC(1,Paso_Solver_Preconditioner);          prec=MEMALLOC(1,Paso_Solver_Preconditioner);
# Line 106  void Paso_Solver_setPreconditioner(Paso_ Line 107  void Paso_Solver_setPreconditioner(Paso_
107                prec->type=PASO_GS;                prec->type=PASO_GS;
108                break;                break;
109              case PASO_AMG:              case PASO_AMG:
110                    
111                if (options->verbose) printf("AMG preconditioner is used.\n");                if (options->verbose) printf("AMG preconditioner is used.\n");
112                                
113                /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/                /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
# Line 115  void Paso_Solver_setPreconditioner(Paso_ Line 117  void Paso_Solver_setPreconditioner(Paso_
117                else {                else {
118                  for (i=0;i<A->row_block_size;++i) {                  for (i=0;i<A->row_block_size;++i) {
119                  prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);                  prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
120                    /*sprintf(filename,"ABlock%d",i);
121                    Paso_SparseMatrix_saveMM(prec->amgSystem->block[i],filename);
122                    */
123                  prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options);                  prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options);
124                  }                  }
125                }                }
126                  
127                              
128                prec->type=PASO_AMG;                prec->type=PASO_AMG;
129                break;                break;
130    
131              case PASO_AMLI:              case PASO_AMLI:
132                if (options->verbose) printf("AMLI preconditioner is used.\n");                if (options->verbose) printf("AMLI preconditioner is used.\n");
133                                
# Line 153  void Paso_Solver_solvePreconditioner(Pas Line 161  void Paso_Solver_solvePreconditioner(Pas
161      dim_t n=A->mainBlock->numRows;      dim_t n=A->mainBlock->numRows;
162      double **xx;      double **xx;
163      double **bb;      double **bb;
164        
165      xx=MEMALLOC(A->row_block_size,double*);      xx=MEMALLOC(A->row_block_size,double*);
166      bb=MEMALLOC(A->row_block_size,double*);      bb=MEMALLOC(A->row_block_size,double*);
167      if (Paso_checkPtr(xx) && Paso_checkPtr(bb)) return;      if (Paso_checkPtr(xx) && Paso_checkPtr(bb)) return;
# Line 182  void Paso_Solver_solvePreconditioner(Pas Line 191  void Paso_Solver_solvePreconditioner(Pas
191                b_old=b_new                b_old=b_new
192                                
193                b_new=b_old+b                b_new=b_old+b
194                  b_new=b_new-Ax_2=b_new-AP^{-1}(2b-AP^{-1}b)
195                x_3=....=P^{-1}(b_new-AP^{-1}b_old)                x_3=....=P^{-1}(b_new-AP^{-1}b_old)
196                b_old=b_new                b_old=b_new
197                                
198                So for n- steps we will use loop, but we always make sure that every value calculated only once!                So for n- steps we will use loop, but we always make sure that every value calculated only once!
199              */              */
200                
201             Paso_Solver_solveGS(prec->gs,x,b);             Paso_Solver_solveGS(prec->gs,x,b);
202             if (prec->gs->sweeps>1) {             if (prec->gs->sweeps>1) {
203                  double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);                  double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
204                  double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);                  double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
205                    dim_t sweeps=prec->gs->sweeps;
206                        
207                  #pragma omp parallel for private(i) schedule(static)                  #pragma omp parallel for private(i) schedule(static)
208                  for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];                  for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
209                        
210                  while(prec->gs->sweeps>1) {                  while(sweeps>1) {
211                     #pragma omp parallel for private(i) schedule(static)                     #pragma omp parallel for private(i) schedule(static)
212                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
213                      /* Compute the residual b=b-Ax*/                      /* Compute the residual b=b-Ax*/
# Line 205  void Paso_Solver_solvePreconditioner(Pas Line 216  void Paso_Solver_solvePreconditioner(Pas
216                     Paso_Solver_solveGS(prec->gs,x,bnew);                     Paso_Solver_solveGS(prec->gs,x,bnew);
217                     #pragma omp parallel for private(i) schedule(static)                     #pragma omp parallel for private(i) schedule(static)
218                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
219                     prec->gs->sweeps=prec->gs->sweeps-1;                     sweeps=sweeps-1;
220                  }                  }
221                  MEMFREE(bold);                  MEMFREE(bold);
222                  MEMFREE(bnew);                  MEMFREE(bnew);
223             }             }
224             break;             break;
225          case PASO_AMG:          case PASO_AMG:
226                
227              /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/              /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
228              if (A->row_block_size==1) {              if (A->row_block_size==1) {
229                  Paso_Solver_solveAMG(prec->amg,x,b);                  Paso_Solver_solveAMG(prec->amg,x,b);

Legend:
Removed from v.2991  
changed lines
  Added in v.2992

  ViewVC Help
Powered by ViewVC 1.1.26