/[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 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 2178 by artak, Thu Dec 18 00:08:58 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 37  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        MEMFREE(in);        MEMFREE(in);
42      }      }
43  }  }
# Line 52  void Paso_Solver_setPreconditioner(Paso_ Line 53  void Paso_Solver_setPreconditioner(Paso_
53          prec->rilu=NULL;          prec->rilu=NULL;
54          prec->ilu=NULL;          prec->ilu=NULL;
55          prec->jacobi=NULL;          prec->jacobi=NULL;
56            prec->gs=NULL;
57            prec->amg=NULL;
58          A->solver=prec;          A->solver=prec;
59          switch (options->preconditioner) {          switch (options->preconditioner) {
60             default:             default:
# Line 70  void Paso_Solver_setPreconditioner(Paso_ Line 73  void Paso_Solver_setPreconditioner(Paso_
73                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
74                prec->type=PASO_RILU;                prec->type=PASO_RILU;
75                break;                break;
76                case PASO_GS:
77                  if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
78                  prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
79                  prec->gs->sweeps=options->sweeps;
80                  prec->type=PASO_GS;
81                  break;
82                case PASO_AMG:
83                  if (options->verbose) printf("AMG preconditioner is used.\n");
84                  prec->amg=Paso_Solver_getAMG(A->mainBlock,options->verbose,0,options->couplingParam);
85                  prec->type=PASO_AMG;
86                  break;
87    
88          }          }
89          if (! Paso_MPIInfo_noError(A->mpi_info ) ){          if (! Paso_MPIInfo_noError(A->mpi_info ) ){
90             Paso_Preconditioner_free(prec);             Paso_Preconditioner_free(prec);
# Line 83  void Paso_Solver_setPreconditioner(Paso_ Line 98  void Paso_Solver_setPreconditioner(Paso_
98  /* 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 */
99  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
100      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
101      #pragma omp barrier      
102      switch (prec->type) {      switch (prec->type) {
103          default:          default:
104          case PASO_JACOBI:          case PASO_JACOBI:
# Line 95  void Paso_Solver_solvePreconditioner(Pas Line 110  void Paso_Solver_solvePreconditioner(Pas
110          case PASO_RILU:          case PASO_RILU:
111             Paso_Solver_solveRILU(prec->rilu,x,b);             Paso_Solver_solveRILU(prec->rilu,x,b);
112             break;             break;
113             case PASO_GS:
114                /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
115                  We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
116                  x_0=0
117                  x_1=P^{-1}(b)
118                  
119                  b_old=b
120                  
121                  b_new=b_old+b
122                  x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
123                     =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
124                  b_old=b_new
125                  
126                  b_new=b_old+b
127                  x_3=....=P^{-1}(b_new-AP^{-1}b_old)
128                  b_old=b_new
129                  
130                  So for n- steps we will use loop, but we always make sure that every value calculated only once!
131                */
132    
133               Paso_Solver_solveGS(prec->gs,x,b);
134               if (prec->gs->sweeps>1) {
135               double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
136               double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
137               dim_t i;
138               #pragma omp parallel for private(i) schedule(static)
139               for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
140              
141               while(prec->gs->sweeps>1) {
142                   #pragma omp parallel for private(i) schedule(static)
143                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
144                    /* Compute the residual b=b-Ax*/
145                   Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
146                   /* Go round again*/
147                   Paso_Solver_solveGS(prec->gs,x,bnew);
148                   #pragma omp parallel for private(i) schedule(static)
149                   for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
150                   prec->gs->sweeps=prec->gs->sweeps-1;
151               }
152              
153               MEMFREE(bold);
154               MEMFREE(bnew);
155              
156               }
157               break;
158             case PASO_AMG:
159               Paso_Solver_solveAMG(prec->amg,x,b);
160               break;
161      }      }
162    
163  }  }

Legend:
Removed from v.1388  
changed lines
  Added in v.2178

  ViewVC Help
Powered by ViewVC 1.1.26