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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 7 months ago) by jfenwick
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 11254 byte(s)
Don't panic.
Updating copyright stamps

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 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 artak 2760 Paso_Solver_AMLI_free(in->amli);
43     Paso_Solver_AMLI_System_free(in->amliSystem);
44 jgs 150 MEMFREE(in);
45     }
46     }
47     /* call the iterative solver: */
48    
49     void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
50 artak 2659
51 jgs 150 Paso_Solver_Preconditioner* prec=NULL;
52 artak 2662 dim_t i;
53 gross 425 if (A->solver==NULL) {
54 jgs 150 /* allocate structure to hold preconditioner */
55     prec=MEMALLOC(1,Paso_Solver_Preconditioner);
56 artak 2659 prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
57 artak 2760 prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);
58 jgs 150 if (Paso_checkPtr(prec)) return;
59 artak 2760
60     if (Paso_checkPtr(prec->amliSystem)) return;
61    
62 artak 2670 if (Paso_checkPtr(prec->amgSystem)) return;
63 artak 2760
64 jgs 150 prec->type=UNKNOWN;
65 gross 430 prec->rilu=NULL;
66 jgs 150 prec->ilu=NULL;
67     prec->jacobi=NULL;
68 artak 1819 prec->gs=NULL;
69 artak 1842 prec->amg=NULL;
70 artak 2760 prec->amli=NULL;
71 artak 2662
72     prec->amgSystem->block_size=A->row_block_size;
73     for (i=0;i<A->row_block_size;++i) {
74     prec->amgSystem->amgblock[i]=NULL;
75     prec->amgSystem->block[i]=NULL;
76     }
77    
78 artak 2760 prec->amliSystem->block_size=A->row_block_size;
79     for (i=0;i<A->row_block_size;++i) {
80     prec->amliSystem->amliblock[i]=NULL;
81     prec->amliSystem->block[i]=NULL;
82     }
83 artak 2662
84 gross 425 A->solver=prec;
85 jgs 150 switch (options->preconditioner) {
86     default:
87     case PASO_JACOBI:
88     if (options->verbose) printf("Jacobi preconditioner is used.\n");
89 ksteube 1312 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
90 jgs 150 prec->type=PASO_JACOBI;
91     break;
92     case PASO_ILU0:
93     if (options->verbose) printf("ILU preconditioner is used.\n");
94 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
95 jgs 150 prec->type=PASO_ILU0;
96     break;
97 gross 430 case PASO_RILU:
98     if (options->verbose) printf("RILU preconditioner is used.\n");
99 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
100 gross 430 prec->type=PASO_RILU;
101     break;
102 artak 1819 case PASO_GS:
103     if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
104     prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
105 artak 1823 prec->gs->sweeps=options->sweeps;
106 artak 1819 prec->type=PASO_GS;
107     break;
108 artak 1842 case PASO_AMG:
109     if (options->verbose) printf("AMG preconditioner is used.\n");
110 artak 2659
111 artak 2662 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
112 artak 2659 if (A->row_block_size==1) {
113     prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
114     }
115 artak 2662 else {
116     for (i=0;i<A->row_block_size;++i) {
117     prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
118     prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options);
119     }
120 artak 2659 }
121 artak 1842 prec->type=PASO_AMG;
122     break;
123 artak 2760 case PASO_AMLI:
124     if (options->verbose) printf("AMLI preconditioner is used.\n");
125    
126     /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
127     if (A->row_block_size==1) {
128     prec->amli=Paso_Solver_getAMLI(A->mainBlock,options->level_max,options);
129     }
130     else {
131     for (i=0;i<A->row_block_size;++i) {
132     prec->amliSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
133     prec->amliSystem->amliblock[i]=Paso_Solver_getAMLI(prec->amliSystem->block[i],options->level_max,options);
134     }
135     }
136     prec->type=PASO_AMLI;
137     break;
138 artak 1842
139 jgs 150 }
140 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
141 artak 2659 Paso_Preconditioner_free(prec);
142     A->solver=NULL;
143 jgs 150 }
144     }
145     }
146    
147     /* applies the preconditioner */
148     /* has to be called within a parallel reqion */
149     /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
150     void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
151 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
152 artak 2662 dim_t i,j;
153 artak 2670 dim_t n=A->mainBlock->numRows;
154 artak 2662 double **xx;
155     double **bb;
156 artak 2670 xx=MEMALLOC(A->row_block_size,double*);
157     bb=MEMALLOC(A->row_block_size,double*);
158     if (Paso_checkPtr(xx) && Paso_checkPtr(bb)) return;
159 artak 2662
160 jgs 150 switch (prec->type) {
161     default:
162     case PASO_JACOBI:
163     Paso_Solver_solveJacobi(prec->jacobi,x,b);
164     break;
165     case PASO_ILU0:
166     Paso_Solver_solveILU(prec->ilu,x,b);
167     break;
168 gross 430 case PASO_RILU:
169     Paso_Solver_solveRILU(prec->rilu,x,b);
170     break;
171 artak 2662 case PASO_GS:
172 artak 1842 /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
173     We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
174     x_0=0
175     x_1=P^{-1}(b)
176    
177     b_old=b
178    
179     b_new=b_old+b
180     x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
181     =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
182     b_old=b_new
183    
184     b_new=b_old+b
185     x_3=....=P^{-1}(b_new-AP^{-1}b_old)
186     b_old=b_new
187    
188     So for n- steps we will use loop, but we always make sure that every value calculated only once!
189     */
190    
191 artak 1819 Paso_Solver_solveGS(prec->gs,x,b);
192 artak 1842 if (prec->gs->sweeps>1) {
193 artak 2662 double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
194     double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
195 artak 1842
196 artak 2662 #pragma omp parallel for private(i) schedule(static)
197     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
198 artak 1842
199 artak 2662 while(prec->gs->sweeps>1) {
200     #pragma omp parallel for private(i) schedule(static)
201     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
202     /* Compute the residual b=b-Ax*/
203     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
204     /* Go round again*/
205     Paso_Solver_solveGS(prec->gs,x,bnew);
206     #pragma omp parallel for private(i) schedule(static)
207     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
208     prec->gs->sweeps=prec->gs->sweeps-1;
209     }
210     MEMFREE(bold);
211     MEMFREE(bnew);
212 artak 1842 }
213 artak 1819 break;
214 artak 2662 case PASO_AMG:
215    
216     /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
217 artak 2659 if (A->row_block_size==1) {
218     Paso_Solver_solveAMG(prec->amg,x,b);
219     }
220 artak 2662 else {
221    
222 artak 2659
223 artak 2670 for (i=0;i<A->row_block_size;i++) {
224 artak 2662 xx[i]=MEMALLOC(n,double);
225     bb[i]=MEMALLOC(n,double);
226 artak 2670 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
227 artak 2659 }
228 artak 2670
229 artak 2705 /*#pragma omp parallel for private(i,j) schedule(static)*/
230 artak 2662 for (i=0;i<n;i++) {
231     for (j=0;j<A->row_block_size;j++) {
232     bb[j][i]=b[A->row_block_size*i+j];
233     }
234 artak 2660 }
235 artak 2670
236    
237 artak 2662 for (i=0;i<A->row_block_size;i++) {
238     Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
239 artak 2659 }
240 artak 2670
241 artak 2705 /*#pragma omp parallel for private(i,j) schedule(static)*/
242 artak 2662 for (i=0;i<n;i++) {
243     for (j=0;j<A->row_block_size;j++) {
244     x[A->row_block_size*i+j]=xx[j][i];
245     }
246     }
247 artak 2670
248 artak 2662 for (i=0;i<A->row_block_size;i++) {
249     MEMFREE(xx[i]);
250     MEMFREE(bb[i]);
251     }
252 artak 2670
253 artak 2659 }
254     break;
255 artak 2760 case PASO_AMLI:
256    
257     /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
258     if (A->row_block_size==1) {
259     Paso_Solver_solveAMLI(prec->amli,x,b);
260     }
261     else {
262    
263    
264     for (i=0;i<A->row_block_size;i++) {
265     xx[i]=MEMALLOC(n,double);
266     bb[i]=MEMALLOC(n,double);
267     if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
268     }
269    
270     /*#pragma omp parallel for private(i,j) schedule(static)*/
271     for (i=0;i<n;i++) {
272     for (j=0;j<A->row_block_size;j++) {
273     bb[j][i]=b[A->row_block_size*i+j];
274     }
275     }
276    
277    
278     for (i=0;i<A->row_block_size;i++) {
279     Paso_Solver_solveAMLI(prec->amliSystem->amliblock[i],xx[i],bb[i]);
280     }
281    
282     /*#pragma omp parallel for private(i,j) schedule(static)*/
283     for (i=0;i<n;i++) {
284     for (j=0;j<A->row_block_size;j++) {
285     x[A->row_block_size*i+j]=xx[j][i];
286     }
287     }
288    
289     for (i=0;i<A->row_block_size;i++) {
290     MEMFREE(xx[i]);
291     MEMFREE(bb[i]);
292     }
293    
294     }
295     break;
296 jgs 150 }
297 artak 2670 MEMFREE(xx);
298     MEMFREE(bb);
299 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26