/[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

trunk/paso/src/Solvers/Solver_preconditioner.c revision 633 by dhawcroft, Thu Mar 23 05:37:00 2006 UTC trunk/paso/src/Solver_preconditioner.c revision 2661 by artak, Fri Sep 11 00:59:59 2009 UTC
# Line 1  Line 1 
 /* $Id$ */  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2009 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
 /*  
 ********************************************************************************  
 *               Copyright   2006 by ACcESS MNRF                                *  
 *                                                                              *  
 *                 http://www.access.edu.au                                     *  
 *           Primary Business: Queensland, Australia                            *  
 *     Licensed under the Open Software License version 3.0             *  
 *        http://www.opensource.org/licenses/osl-3.0.php                        *  
 ********************************************************************************  
 */  
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 19  Line 19 
19  /**************************************************************/  /**************************************************************/
20    
21  /* Copyrights by ACcESS Australia 2003/04 */  /* Copyrights by ACcESS Australia 2003/04 */
22  /* Author: gross@access.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
23    
24  /**************************************************************/  /**************************************************************/
25    
# Line 36  void Paso_Preconditioner_free(Paso_Solve Line 36  void Paso_Preconditioner_free(Paso_Solve
36        Paso_Solver_ILU_free(in->ilu);        Paso_Solver_ILU_free(in->ilu);
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);
40          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;
57          prec->ilu=NULL;          prec->ilu=NULL;
58          prec->jacobi=NULL;          prec->jacobi=NULL;
59            prec->gs=NULL;
60            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:
70             case PASO_JACOBI:             case PASO_JACOBI:
71                if (options->verbose) printf("Jacobi preconditioner is used.\n");                if (options->verbose) printf("Jacobi preconditioner is used.\n");
72                prec->jacobi=Paso_Solver_getJacobi(A);                prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
73                prec->type=PASO_JACOBI;                prec->type=PASO_JACOBI;
74                break;                break;
75             case PASO_ILU0:             case PASO_ILU0:
76                if (options->verbose) printf("ILU preconditioner is used.\n");                if (options->verbose) printf("ILU preconditioner is used.\n");
77                prec->ilu=Paso_Solver_getILU(A,options->verbose);                prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
78                prec->type=PASO_ILU0;                prec->type=PASO_ILU0;
79                break;                break;
80             case PASO_RILU:             case PASO_RILU:
81                if (options->verbose) printf("RILU preconditioner is used.\n");                if (options->verbose) printf("RILU preconditioner is used.\n");
82                prec->rilu=Paso_Solver_getRILU(A,options->verbose);                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
83                prec->type=PASO_RILU;                prec->type=PASO_RILU;
84                break;                break;
85                case PASO_GS:
86                  if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
87                  prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
88                  prec->gs->sweeps=options->sweeps;
89                  prec->type=PASO_GS;
90                  break;
91                case PASO_AMG:
92                  if (options->verbose) printf("AMG preconditioner is used.\n");
93                  
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;
123                  break;
124    
125          }          }
126          if (! Paso_noError()) {          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 82  void Paso_Solver_setPreconditioner(Paso_ Line 135  void Paso_Solver_setPreconditioner(Paso_
135  /* 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 */
136  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
137      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
138      #pragma omp barrier      
139      switch (prec->type) {      switch (prec->type) {
140          default:          default:
141          case PASO_JACOBI:          case PASO_JACOBI:
# Line 94  void Paso_Solver_solvePreconditioner(Pas Line 147  void Paso_Solver_solvePreconditioner(Pas
147          case PASO_RILU:          case PASO_RILU:
148             Paso_Solver_solveRILU(prec->rilu,x,b);             Paso_Solver_solveRILU(prec->rilu,x,b);
149             break;             break;
150             case PASO_GS:
151                /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
152                  We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
153                  x_0=0
154                  x_1=P^{-1}(b)
155                  
156                  b_old=b
157                  
158                  b_new=b_old+b
159                  x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
160                     =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
161                  b_old=b_new
162                  
163                  b_new=b_old+b
164                  x_3=....=P^{-1}(b_new-AP^{-1}b_old)
165                  b_old=b_new
166                  
167                  So for n- steps we will use loop, but we always make sure that every value calculated only once!
168                */
169    
170               Paso_Solver_solveGS(prec->gs,x,b);
171               if (prec->gs->sweeps>1) {
172               double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
173               double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
174               dim_t i;
175               #pragma omp parallel for private(i) schedule(static)
176               for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
177              
178               while(prec->gs->sweeps>1) {
179                   #pragma omp parallel for private(i) schedule(static)
180                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
181                    /* Compute the residual b=b-Ax*/
182                   Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
183                   /* Go round again*/
184                   Paso_Solver_solveGS(prec->gs,x,bnew);
185                   #pragma omp parallel for private(i) schedule(static)
186                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
187                   prec->gs->sweeps=prec->gs->sweeps-1;
188               }
189              
190               MEMFREE(bold);
191               MEMFREE(bnew);
192              
193               }
194               break;
195             case PASO_AMG:
196                if (A->row_block_size==1) {
197                    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                        b2[i]=b[2*i+1];
210                    }
211                    
212                    
213                    Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
214                    Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
215                    
216                    #pragma omp parallel for private(i) schedule(static)
217                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
218                        x[2*i]=x1[i];
219                        x[2*i+1]=x2[i];
220                    }
221                    
222                    MEMFREE(x1);
223                    MEMFREE(x2);
224                    MEMFREE(b1);
225                    MEMFREE(b2);
226    
227                }
228                else if (A->row_block_size==3) {
229                    dim_t i;
230                    
231                    double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
232                    double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
233                    double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
234                    double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
235                    double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
236                    double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
237                    
238                    #pragma omp parallel for private(i) schedule(static)
239                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
240                        b1[i]=b[3*i];
241                        b2[i]=b[3*i+1];
242                        b3[i]=b[3*i+2];
243                     }
244                    
245                    Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
246                    Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
247                    Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);
248                    
249                    #pragma omp parallel for private(i) schedule(static)
250                    for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
251                        x[3*i]=x1[i];
252                        x[3*i+1]=x2[i];
253                        x[3*i+2]=x3[i];
254                    }
255                    
256                    
257                    MEMFREE(x1);
258                    MEMFREE(x2);
259                    MEMFREE(x3);
260                    MEMFREE(b1);
261                    MEMFREE(b2);
262                    MEMFREE(b3);
263                }
264            break;
265      }      }
 }  
266    
267  /*  }
  * $Log$  
  * Revision 1.2  2005/09/15 03:44:40  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.1.2.2  2005/09/07 00:59:09  gross  
  * some inconsistent renaming fixed to make the linking work.  
  *  
  * Revision 1.1.2.1  2005/09/05 06:29:50  gross  
  * These files have been extracted from finley to define a stand alone libray for iterative  
  * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  
  * has not been tested yet.  
  *  
  *  
  */  

Legend:
Removed from v.633  
changed lines
  Added in v.2661

  ViewVC Help
Powered by ViewVC 1.1.26