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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26