/[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 2658 by jfenwick, Tue Aug 18 01:25:18 2009 UTC revision 2659 by artak, Thu Sep 10 03:54:50 2009 UTC
# Line 38  void Paso_Preconditioner_free(Paso_Solve Line 38  void Paso_Preconditioner_free(Paso_Solve
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);        Paso_Solver_AMG_free(in->amg);
41          Paso_Solver_AMG_System_free(in->amgSystem);
42        MEMFREE(in);        MEMFREE(in);
43      }      }
44  }  }
45  /*  call the iterative solver: */  /*  call the iterative solver: */
46    
47  void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {  void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48    
49      Paso_Solver_Preconditioner* prec=NULL;      Paso_Solver_Preconditioner* prec=NULL;
50      if (A->solver==NULL) {      if (A->solver==NULL) {
51          /* allocate structure to hold preconditioner */          /* allocate structure to hold preconditioner */
52          prec=MEMALLOC(1,Paso_Solver_Preconditioner);          prec=MEMALLOC(1,Paso_Solver_Preconditioner);
53            prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
54          if (Paso_checkPtr(prec)) return;          if (Paso_checkPtr(prec)) return;
55          prec->type=UNKNOWN;          prec->type=UNKNOWN;
56          prec->rilu=NULL;          prec->rilu=NULL;
# Line 55  void Paso_Solver_setPreconditioner(Paso_ Line 58  void Paso_Solver_setPreconditioner(Paso_
58          prec->jacobi=NULL;          prec->jacobi=NULL;
59          prec->gs=NULL;          prec->gs=NULL;
60          prec->amg=NULL;          prec->amg=NULL;
61            prec->amgSystem->amgblock1=NULL;
62            prec->amgSystem->amgblock2=NULL;
63            prec->amgSystem->amgblock3=NULL;
64            prec->amgSystem->block1=NULL;
65            prec->amgSystem->block2=NULL;
66            prec->amgSystem->block3=NULL;
67          A->solver=prec;          A->solver=prec;
68          switch (options->preconditioner) {          switch (options->preconditioner) {
69             default:             default:
# Line 81  void Paso_Solver_setPreconditioner(Paso_ Line 90  void Paso_Solver_setPreconditioner(Paso_
90                break;                break;
91              case PASO_AMG:              case PASO_AMG:
92                if (options->verbose) printf("AMG preconditioner is used.\n");                if (options->verbose) printf("AMG preconditioner is used.\n");
93                prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);                
94                  if (A->row_block_size==1) {
95                    prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);  
96                  }
97                  else if (A->row_block_size==2) {
98                    
99                    prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
100                    prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);
101    
102                    prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
103                    prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
104                    
105                  }
106                  else if (A->row_block_size==3)
107                  {
108    
109                    prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
110                    prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);
111                    prec->amgSystem->block3=Paso_SparseMatrix_getBlock(A->mainBlock,3);
112                    
113                    /*
114                    Paso_SparseMatrix_saveMM(prec->amgSystem->block1,"amgtemp1.mat");
115                    Paso_SparseMatrix_saveMM(prec->amgSystem->block2,"amgtemp2.mat");
116                    Paso_SparseMatrix_saveMM(prec->amgSystem->block3,"amgtemp3.mat");
117                    */
118                    prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
119                    prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
120                    prec->amgSystem->amgblock3=Paso_Solver_getAMG(prec->amgSystem->block3,options->level_max,options);
121                  }
122                prec->type=PASO_AMG;                prec->type=PASO_AMG;
123                break;                break;
124    
125          }          }
126          if (! Paso_MPIInfo_noError(A->mpi_info ) ){          if (! Paso_MPIInfo_noError(A->mpi_info ) ){
127             Paso_Preconditioner_free(prec);              Paso_Preconditioner_free(prec);
128             A->solver=NULL;              A->solver=NULL;
129          }          }
130      }      }
131  }  }
# Line 156  void Paso_Solver_solvePreconditioner(Pas Line 193  void Paso_Solver_solvePreconditioner(Pas
193             }             }
194             break;             break;
195           case PASO_AMG:           case PASO_AMG:
196             Paso_Solver_solveAMG(prec->amg,x,b);              if (A->row_block_size==1) {
197             break;                  Paso_Solver_solveAMG(prec->amg,x,b);
198                }
199                else if (A->row_block_size==2) {
200                    dim_t i;
201                    double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
202                    double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
203                    double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
204                    double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
205                    
206                    #pragma omp parallel for private(i) schedule(static)
207                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
208                        b1[i]=b[2*i];
209                    }
210                    
211                    #pragma omp parallel for private(i) schedule(static)
212                    for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
213                        b2[i]=b[2*i+1];
214                    }
215                    
216                    Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
217                    Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
218                    
219                    #pragma omp parallel for private(i) schedule(static)
220                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
221                        x[2*i]=x1[i];
222                    }
223                    
224                    #pragma omp parallel for private(i) schedule(static)
225                    for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
226                        x[2*i+1]=x2[i];
227                    }
228                                                    
229                    MEMFREE(x1);
230                    MEMFREE(x2);
231                    MEMFREE(b1);
232                    MEMFREE(b2);
233    
234                }
235                else if (A->row_block_size==3) {
236                    dim_t i;
237                    
238                    double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
239                    double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
240                    double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
241                    double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
242                    double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
243                    double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
244                    
245                    #pragma omp parallel for private(i) schedule(static)
246                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
247                        b1[i]=b[3*i];
248                     }
249                    
250                    #pragma omp parallel for private(i) schedule(static)
251                    for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
252                        b2[i]=b[3*i+1];
253                    }
254                    
255                    #pragma omp parallel for private(i) schedule(static)
256                    for (i=0;i<prec->amgSystem->amgblock3->n;i++) {
257                        b3[i]=b[3*i+2];
258                    }
259                    
260                    
261                    Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
262                    Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
263                    Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);
264                    
265                    #pragma omp parallel for private(i) schedule(static)
266                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
267                        x[3*i]=x1[i];
268                    }
269                    
270                    #pragma omp parallel for private(i) schedule(static)
271                    for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
272                        x[3*i+1]=x2[i];
273                    }
274                    
275                    #pragma omp parallel for private(i) schedule(static)
276                    for (i=0;i<prec->amgSystem->amgblock3->n;i++) {
277                        x[3*i+2]=x3[i];
278                    }
279                    
280                    MEMFREE(x1);
281                    MEMFREE(x2);
282                    MEMFREE(x3);
283                    MEMFREE(b1);
284                    MEMFREE(b2);
285                    MEMFREE(b3);
286                }
287            break;
288      }      }
289    
290  }  }

Legend:
Removed from v.2658  
changed lines
  Added in v.2659

  ViewVC Help
Powered by ViewVC 1.1.26