/[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 2661 - (hide annotations)
Fri Sep 11 00:59:59 2009 UTC (10 years, 1 month ago) by artak
File MIME type: text/plain
File size: 10552 byte(s)
Finially unit tests for AMG preconditioner for single equations as well as for systems. Minor bug in AMG is also fixed
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
23 jgs 150
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 artak 2659 Paso_Solver_AMG_System_free(in->amgSystem);
42 jgs 150 MEMFREE(in);
43     }
44     }
45     /* call the iterative solver: */
46    
47     void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48 artak 2659
49 jgs 150 Paso_Solver_Preconditioner* prec=NULL;
50 gross 425 if (A->solver==NULL) {
51 jgs 150 /* allocate structure to hold preconditioner */
52     prec=MEMALLOC(1,Paso_Solver_Preconditioner);
53 artak 2659 prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
54 jgs 150 if (Paso_checkPtr(prec)) return;
55     prec->type=UNKNOWN;
56 gross 430 prec->rilu=NULL;
57 jgs 150 prec->ilu=NULL;
58     prec->jacobi=NULL;
59 artak 1819 prec->gs=NULL;
60 artak 1842 prec->amg=NULL;
61 artak 2659 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 gross 425 A->solver=prec;
68 jgs 150 switch (options->preconditioner) {
69     default:
70     case PASO_JACOBI:
71     if (options->verbose) printf("Jacobi preconditioner is used.\n");
72 ksteube 1312 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
73 jgs 150 prec->type=PASO_JACOBI;
74     break;
75     case PASO_ILU0:
76     if (options->verbose) printf("ILU preconditioner is used.\n");
77 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
78 jgs 150 prec->type=PASO_ILU0;
79     break;
80 gross 430 case PASO_RILU:
81     if (options->verbose) printf("RILU preconditioner is used.\n");
82 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
83 gross 430 prec->type=PASO_RILU;
84     break;
85 artak 1819 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 artak 1823 prec->gs->sweeps=options->sweeps;
89 artak 1819 prec->type=PASO_GS;
90     break;
91 artak 1842 case PASO_AMG:
92     if (options->verbose) printf("AMG preconditioner is used.\n");
93 artak 2659
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 artak 1842 prec->type=PASO_AMG;
123     break;
124    
125 jgs 150 }
126 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
127 artak 2659 Paso_Preconditioner_free(prec);
128     A->solver=NULL;
129 jgs 150 }
130     }
131     }
132    
133     /* applies the preconditioner */
134     /* has to be called within a parallel reqion */
135     /* 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){
137 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
138 artak 1842
139 jgs 150 switch (prec->type) {
140     default:
141     case PASO_JACOBI:
142     Paso_Solver_solveJacobi(prec->jacobi,x,b);
143     break;
144     case PASO_ILU0:
145     Paso_Solver_solveILU(prec->ilu,x,b);
146     break;
147 gross 430 case PASO_RILU:
148     Paso_Solver_solveRILU(prec->rilu,x,b);
149     break;
150 artak 1819 case PASO_GS:
151 artak 1842 /* 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 artak 1819 Paso_Solver_solveGS(prec->gs,x,b);
171 artak 1842 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 artak 1819 break;
195 artak 1842 case PASO_AMG:
196 artak 2659 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 artak 2660
213 artak 2659 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 artak 2660 x[2*i+1]=x2[i];
220 artak 2659 }
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 artak 2660 }
244 artak 2659
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 artak 2661
257 artak 2659 MEMFREE(x1);
258     MEMFREE(x2);
259     MEMFREE(x3);
260     MEMFREE(b1);
261     MEMFREE(b2);
262     MEMFREE(b3);
263     }
264     break;
265 jgs 150 }
266 artak 1842
267 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26