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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26