/[escript]/trunk/paso/src/Solver_preconditioner.c
ViewVC logotype

Annotation of /trunk/paso/src/Solver_preconditioner.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2475 - (hide annotations)
Wed Jun 17 01:48:46 2009 UTC (10 years, 4 months ago) by artak
File MIME type: text/plain
File size: 5940 byte(s)
AMG now takes into account new SolverOptions class, namely one can select coarsening algorithms from python. Not all options are included into AMG, this will be implemented later on.
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 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 dhawcroft 631
14 ksteube 1811
15 jgs 150 /**************************************************************/
16    
17     /* Paso: SystemMatrix: sets-up the preconditioner */
18    
19     /**************************************************************/
20    
21     /* Copyrights by ACcESS Australia 2003/04 */
22     /* Author: gross@access.edu.au */
23    
24     /**************************************************************/
25    
26 gross 700 #include "Paso.h"
27     #include "SystemMatrix.h"
28 jgs 150 #include "Solver.h"
29    
30     /***********************************************************************************/
31    
32     /* free space */
33    
34     void Paso_Preconditioner_free(Paso_Solver_Preconditioner* in) {
35     if (in!=NULL) {
36     Paso_Solver_ILU_free(in->ilu);
37 gross 430 Paso_Solver_RILU_free(in->rilu);
38 jgs 150 Paso_Solver_Jacobi_free(in->jacobi);
39 artak 1819 Paso_Solver_GS_free(in->gs);
40 artak 1842 Paso_Solver_AMG_free(in->amg);
41 jgs 150 MEMFREE(in);
42     }
43     }
44     /* call the iterative solver: */
45    
46     void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
47     Paso_Solver_Preconditioner* prec=NULL;
48 gross 425 if (A->solver==NULL) {
49 jgs 150 /* allocate structure to hold preconditioner */
50     prec=MEMALLOC(1,Paso_Solver_Preconditioner);
51     if (Paso_checkPtr(prec)) return;
52     prec->type=UNKNOWN;
53 gross 430 prec->rilu=NULL;
54 jgs 150 prec->ilu=NULL;
55     prec->jacobi=NULL;
56 artak 1819 prec->gs=NULL;
57 artak 1842 prec->amg=NULL;
58 gross 425 A->solver=prec;
59 jgs 150 switch (options->preconditioner) {
60     default:
61     case PASO_JACOBI:
62     if (options->verbose) printf("Jacobi preconditioner is used.\n");
63 ksteube 1312 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
64 jgs 150 prec->type=PASO_JACOBI;
65     break;
66     case PASO_ILU0:
67     if (options->verbose) printf("ILU preconditioner is used.\n");
68 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
69 jgs 150 prec->type=PASO_ILU0;
70     break;
71 gross 430 case PASO_RILU:
72     if (options->verbose) printf("RILU preconditioner is used.\n");
73 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
74 gross 430 prec->type=PASO_RILU;
75     break;
76 artak 1819 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 artak 1823 prec->gs->sweeps=options->sweeps;
80 artak 1819 prec->type=PASO_GS;
81     break;
82 artak 1842 case PASO_AMG:
83     if (options->verbose) printf("AMG preconditioner is used.\n");
84 artak 2475 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->verbose,options->level_max,options->coarsening_threshold,options->coarsening_method);
85 artak 1842 prec->type=PASO_AMG;
86     break;
87    
88 jgs 150 }
89 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
90 jgs 150 Paso_Preconditioner_free(prec);
91 gross 425 A->solver=NULL;
92 jgs 150 }
93     }
94     }
95    
96     /* applies the preconditioner */
97     /* has to be called within a parallel reqion */
98     /* 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){
100 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
101 artak 1842
102 jgs 150 switch (prec->type) {
103     default:
104     case PASO_JACOBI:
105     Paso_Solver_solveJacobi(prec->jacobi,x,b);
106     break;
107     case PASO_ILU0:
108     Paso_Solver_solveILU(prec->ilu,x,b);
109     break;
110 gross 430 case PASO_RILU:
111     Paso_Solver_solveRILU(prec->rilu,x,b);
112     break;
113 artak 1819 case PASO_GS:
114 artak 1842 /* 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 artak 1819 Paso_Solver_solveGS(prec->gs,x,b);
134 artak 1842 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 artak 1819 break;
158 artak 1842 case PASO_AMG:
159     Paso_Solver_solveAMG(prec->amg,x,b);
160     break;
161 jgs 150 }
162 artak 1842
163 jgs 150 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26