/[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 2669 by artak, Thu Sep 17 00:10:00 2009 UTC revision 2670 by artak, Thu Sep 17 06:17:57 2009 UTC
# Line 53  void Paso_Solver_setPreconditioner(Paso_ Line 53  void Paso_Solver_setPreconditioner(Paso_
53          prec=MEMALLOC(1,Paso_Solver_Preconditioner);          prec=MEMALLOC(1,Paso_Solver_Preconditioner);
54          prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);          prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
55          if (Paso_checkPtr(prec)) return;          if (Paso_checkPtr(prec)) return;
56            if (Paso_checkPtr(prec->amgSystem)) return;
57          prec->type=UNKNOWN;          prec->type=UNKNOWN;
58          prec->rilu=NULL;          prec->rilu=NULL;
59          prec->ilu=NULL;          prec->ilu=NULL;
# Line 104  void Paso_Solver_setPreconditioner(Paso_ Line 105  void Paso_Solver_setPreconditioner(Paso_
105                  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);
106                  }                  }
107                }                }
   
108                prec->type=PASO_AMG;                prec->type=PASO_AMG;
109                break;                break;
110    
# Line 122  void Paso_Solver_setPreconditioner(Paso_ Line 122  void Paso_Solver_setPreconditioner(Paso_
122  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
123      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
124      dim_t i,j;      dim_t i,j;
125      dim_t n=Paso_SystemMatrix_getGlobalNumRows(A);      dim_t n=A->mainBlock->numRows;
126      double **xx;      double **xx;
127      double **bb;      double **bb;
128        xx=MEMALLOC(A->row_block_size,double*);
129        bb=MEMALLOC(A->row_block_size,double*);
130        if (Paso_checkPtr(xx) && Paso_checkPtr(bb)) return;
131    
132      switch (prec->type) {      switch (prec->type) {
133          default:          default:
# Line 188  void Paso_Solver_solvePreconditioner(Pas Line 191  void Paso_Solver_solvePreconditioner(Pas
191              }              }
192              else {              else {
193                        
                 xx=MEMALLOC(A->row_block_size,double*);  
                 bb=MEMALLOC(A->row_block_size,double*);  
194                                    
195                  for (i=0;i<A->row_block_size;i++) {                   for (i=0;i<A->row_block_size;i++) {
196                      xx[i]=MEMALLOC(n,double);                      xx[i]=MEMALLOC(n,double);
197                      bb[i]=MEMALLOC(n,double);                      bb[i]=MEMALLOC(n,double);
198                        if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
199                  }                  }
200                    
201                  #pragma omp parallel for private(i,j) schedule(static)                  #pragma omp parallel for private(i,j) schedule(static)
202                  for (i=0;i<n;i++) {                  for (i=0;i<n;i++) {
203                      for (j=0;j<A->row_block_size;j++) {                      for (j=0;j<A->row_block_size;j++) {
204                       bb[j][i]=b[A->row_block_size*i+j];                       bb[j][i]=b[A->row_block_size*i+j];
205                      }                      }
206                   }                   }
207                                                    
208                    
209                  for (i=0;i<A->row_block_size;i++) {                  for (i=0;i<A->row_block_size;i++) {
210                  Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);                  Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
211                  }                  }
212                                  
213                  #pragma omp parallel for private(i,j) schedule(static)                  #pragma omp parallel for private(i,j) schedule(static)
214                  for (i=0;i<n;i++) {                  for (i=0;i<n;i++) {
215                      for (j=0;j<A->row_block_size;j++) {                      for (j=0;j<A->row_block_size;j++) {
216                      x[A->row_block_size*i+j]=xx[j][i];                      x[A->row_block_size*i+j]=xx[j][i];
217                      }                      }
218                   }                   }
219                                    
220                  for (i=0;i<A->row_block_size;i++) {                  for (i=0;i<A->row_block_size;i++) {
221                  MEMFREE(xx[i]);                  MEMFREE(xx[i]);
222                  MEMFREE(bb[i]);                  MEMFREE(bb[i]);
223                  }                  }
224                                  
                 MEMFREE(xx);  
                 MEMFREE(bb);  
225              }              }
               
226          break;          break;
227      }      }
228        MEMFREE(xx);
229        MEMFREE(bb);
230  }  }

Legend:
Removed from v.2669  
changed lines
  Added in v.2670

  ViewVC Help
Powered by ViewVC 1.1.26