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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3013 - (show annotations)
Wed Apr 28 03:47:51 2010 UTC (9 years, 1 month ago) by artak
File MIME type: text/plain
File size: 11109 byte(s)
Block separation part commented out for AMG preconditioner
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
48
49 /* call the iterative solver: */
50
51 void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
52
53 Paso_Solver_Preconditioner* prec=NULL;
54 dim_t i;
55 /*char filename[7];*/
56 if (A->solver==NULL) {
57 /* allocate structure to hold preconditioner */
58 prec=MEMALLOC(1,Paso_Solver_Preconditioner);
59
60 if (Paso_checkPtr(prec)) return;
61
62 prec->type=UNKNOWN;
63 prec->rilu=NULL;
64 prec->ilu=NULL;
65 prec->jacobi=NULL;
66 prec->gs=NULL;
67 prec->amg=NULL;
68 prec->amli=NULL;
69 prec->amgSystem=NULL;
70 prec->amliSystem=NULL;
71
72 A->solver=prec;
73 switch (options->preconditioner) {
74 default:
75 case PASO_JACOBI:
76 if (options->verbose) printf("Jacobi preconditioner is used.\n");
77 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
78 prec->type=PASO_JACOBI;
79 break;
80 case PASO_ILU0:
81 if (options->verbose) printf("ILU preconditioner is used.\n");
82 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
83 prec->type=PASO_ILU0;
84 break;
85 case PASO_RILU:
86 if (options->verbose) printf("RILU preconditioner is used.\n");
87 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
88 prec->type=PASO_RILU;
89 break;
90 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 prec->gs->sweeps=options->sweeps;
94 prec->type=PASO_GS;
95 break;
96 case PASO_AMG:
97 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 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 */
108
109 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
110 /*if (A->row_block_size==1) {
111 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
112 }
113 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 }
119 */
120
121 prec->type=PASO_AMG;
122 break;
123
124 case PASO_AMLI:
125 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
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
148 }
149 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
150 Paso_Preconditioner_free(prec);
151 A->solver=NULL;
152 }
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 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
161 dim_t i,j;
162 dim_t n=A->mainBlock->numRows;
163 double **xx;
164 double **bb;
165
166 xx=MEMALLOC(A->row_block_size,double*);
167 bb=MEMALLOC(A->row_block_size,double*);
168 if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
169
170 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 case PASO_RILU:
179 Paso_Solver_solveRILU(prec->rilu,x,b);
180 break;
181 case PASO_GS:
182 /* 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 b_new=b_new-Ax_2=b_new-AP^{-1}(2b-AP^{-1}b)
196 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
202 Paso_Solver_solveGS(prec->gs,x,b);
203 if (prec->gs->sweeps>1) {
204 double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
205 dim_t sweeps=prec->gs->sweeps;
206
207 #pragma omp parallel for private(i) schedule(static)
208 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=b[i];
209
210 while(sweeps>1) {
211 #pragma omp parallel for private(i) schedule(static)
212 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]+=b[i];
213 /* Compute the residual b=b-Ax*/
214 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A->mainBlock, x, DBLE(1), bnew);
215 /* Go round again*/
216 Paso_Solver_solveGS(prec->gs,x,bnew);
217 sweeps--;
218 }
219 MEMFREE(bnew);
220 }
221 break;
222 case PASO_AMG:
223 Paso_Solver_solveAMG(prec->amg,x,b);
224
225 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
226 /*
227 if (A->row_block_size<4) {
228 Paso_Solver_solveAMG(prec->amg,x,b);
229 }
230 else {
231
232
233 for (i=0;i<A->row_block_size;i++) {
234 xx[i]=MEMALLOC(n,double);
235 bb[i]=MEMALLOC(n,double);
236 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
237 }
238
239 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 }
244
245
246 for (i=0;i<A->row_block_size;i++) {
247 Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
248 }
249
250 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
256 for (i=0;i<A->row_block_size;i++) {
257 MEMFREE(xx[i]);
258 MEMFREE(bb[i]);
259 }
260 }
261 */
262 break;
263 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 }
305 MEMFREE(xx);
306 MEMFREE(bb);
307 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26