/[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 2705 by artak, Thu Oct 1 23:21:10 2009 UTC revision 2760 by artak, Thu Nov 19 05:22:45 2009 UTC
# Line 39  void Paso_Preconditioner_free(Paso_Solve Line 39  void Paso_Preconditioner_free(Paso_Solve
39        Paso_Solver_GS_free(in->gs);        Paso_Solver_GS_free(in->gs);
40        Paso_Solver_AMG_free(in->amg);        Paso_Solver_AMG_free(in->amg);
41        Paso_Solver_AMG_System_free(in->amgSystem);        Paso_Solver_AMG_System_free(in->amgSystem);
42          Paso_Solver_AMLI_free(in->amli);
43          Paso_Solver_AMLI_System_free(in->amliSystem);
44        MEMFREE(in);        MEMFREE(in);
45      }      }
46  }  }
# Line 52  void Paso_Solver_setPreconditioner(Paso_ Line 54  void Paso_Solver_setPreconditioner(Paso_
54          /* allocate structure to hold preconditioner */          /* allocate structure to hold preconditioner */
55          prec=MEMALLOC(1,Paso_Solver_Preconditioner);          prec=MEMALLOC(1,Paso_Solver_Preconditioner);
56          prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);          prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
57            prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);
58          if (Paso_checkPtr(prec)) return;          if (Paso_checkPtr(prec)) return;
59            
60            if (Paso_checkPtr(prec->amliSystem)) return;
61            
62          if (Paso_checkPtr(prec->amgSystem)) return;          if (Paso_checkPtr(prec->amgSystem)) return;
63            
64          prec->type=UNKNOWN;          prec->type=UNKNOWN;
65          prec->rilu=NULL;          prec->rilu=NULL;
66          prec->ilu=NULL;          prec->ilu=NULL;
67          prec->jacobi=NULL;          prec->jacobi=NULL;
68          prec->gs=NULL;          prec->gs=NULL;
69          prec->amg=NULL;          prec->amg=NULL;
70            prec->amli=NULL;
71                    
72          prec->amgSystem->block_size=A->row_block_size;          prec->amgSystem->block_size=A->row_block_size;
73          for (i=0;i<A->row_block_size;++i) {          for (i=0;i<A->row_block_size;++i) {
# Line 67  void Paso_Solver_setPreconditioner(Paso_ Line 75  void Paso_Solver_setPreconditioner(Paso_
75            prec->amgSystem->block[i]=NULL;            prec->amgSystem->block[i]=NULL;
76          }          }
77                    
78            prec->amliSystem->block_size=A->row_block_size;
79            for (i=0;i<A->row_block_size;++i) {
80              prec->amliSystem->amliblock[i]=NULL;
81              prec->amliSystem->block[i]=NULL;
82            }
83    
84          A->solver=prec;          A->solver=prec;
85          switch (options->preconditioner) {          switch (options->preconditioner) {
# Line 107  void Paso_Solver_setPreconditioner(Paso_ Line 120  void Paso_Solver_setPreconditioner(Paso_
120                }                }
121                prec->type=PASO_AMG;                prec->type=PASO_AMG;
122                break;                break;
123                case PASO_AMLI:
124                  if (options->verbose) printf("AMLI preconditioner is used.\n");
125                  
126                  /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
127                  if (A->row_block_size==1) {
128                    prec->amli=Paso_Solver_getAMLI(A->mainBlock,options->level_max,options);  
129                  }
130                  else {
131                    for (i=0;i<A->row_block_size;++i) {
132                    prec->amliSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
133                    prec->amliSystem->amliblock[i]=Paso_Solver_getAMLI(prec->amliSystem->block[i],options->level_max,options);
134                    }
135                  }
136                  prec->type=PASO_AMLI;
137                  break;
138    
139          }          }
140          if (! Paso_MPIInfo_noError(A->mpi_info ) ){          if (! Paso_MPIInfo_noError(A->mpi_info ) ){
# Line 211  void Paso_Solver_solvePreconditioner(Pas Line 239  void Paso_Solver_solvePreconditioner(Pas
239                  }                  }
240                                                                
241                  /*#pragma omp parallel for private(i,j) schedule(static)*/                  /*#pragma omp parallel for private(i,j) schedule(static)*/
242                    for (i=0;i<n;i++) {
243                        for (j=0;j<A->row_block_size;j++) {
244                        x[A->row_block_size*i+j]=xx[j][i];
245                        }
246                     }
247                    
248                    for (i=0;i<A->row_block_size;i++) {
249                    MEMFREE(xx[i]);
250                    MEMFREE(bb[i]);
251                    }
252                  
253                }
254            break;
255            case PASO_AMLI:
256                
257                /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
258                if (A->row_block_size==1) {
259                    Paso_Solver_solveAMLI(prec->amli,x,b);
260                }
261                else {
262              
263                    
264                     for (i=0;i<A->row_block_size;i++) {
265                        xx[i]=MEMALLOC(n,double);
266                        bb[i]=MEMALLOC(n,double);
267                        if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
268                    }
269                    
270                    /*#pragma omp parallel for private(i,j) schedule(static)*/
271                    for (i=0;i<n;i++) {
272                        for (j=0;j<A->row_block_size;j++) {
273                         bb[j][i]=b[A->row_block_size*i+j];
274                        }
275                     }
276                    
277                    
278                    for (i=0;i<A->row_block_size;i++) {
279                    Paso_Solver_solveAMLI(prec->amliSystem->amliblock[i],xx[i],bb[i]);
280                    }
281                                  
282                    /*#pragma omp parallel for private(i,j) schedule(static)*/
283                  for (i=0;i<n;i++) {                  for (i=0;i<n;i++) {
284                      for (j=0;j<A->row_block_size;j++) {                      for (j=0;j<A->row_block_size;j++) {
285                      x[A->row_block_size*i+j]=xx[j][i];                      x[A->row_block_size*i+j]=xx[j][i];

Legend:
Removed from v.2705  
changed lines
  Added in v.2760

  ViewVC Help
Powered by ViewVC 1.1.26