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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2659 - (show annotations)
Thu Sep 10 03:54:50 2009 UTC (10 years ago) by artak
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 11575 byte(s)
first attempt of AMG on blocks
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 MEMFREE(in);
43 }
44 }
45 /* call the iterative solver: */
46
47 void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48
49 Paso_Solver_Preconditioner* prec=NULL;
50 if (A->solver==NULL) {
51 /* allocate structure to hold preconditioner */
52 prec=MEMALLOC(1,Paso_Solver_Preconditioner);
53 prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
54 if (Paso_checkPtr(prec)) return;
55 prec->type=UNKNOWN;
56 prec->rilu=NULL;
57 prec->ilu=NULL;
58 prec->jacobi=NULL;
59 prec->gs=NULL;
60 prec->amg=NULL;
61 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 A->solver=prec;
68 switch (options->preconditioner) {
69 default:
70 case PASO_JACOBI:
71 if (options->verbose) printf("Jacobi preconditioner is used.\n");
72 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
73 prec->type=PASO_JACOBI;
74 break;
75 case PASO_ILU0:
76 if (options->verbose) printf("ILU preconditioner is used.\n");
77 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
78 prec->type=PASO_ILU0;
79 break;
80 case PASO_RILU:
81 if (options->verbose) printf("RILU preconditioner is used.\n");
82 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
83 prec->type=PASO_RILU;
84 break;
85 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 prec->gs->sweeps=options->sweeps;
89 prec->type=PASO_GS;
90 break;
91 case PASO_AMG:
92 if (options->verbose) printf("AMG preconditioner is used.\n");
93
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 prec->type=PASO_AMG;
123 break;
124
125 }
126 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
127 Paso_Preconditioner_free(prec);
128 A->solver=NULL;
129 }
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 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
138
139 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 case PASO_RILU:
148 Paso_Solver_solveRILU(prec->rilu,x,b);
149 break;
150 case PASO_GS:
151 /* 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 Paso_Solver_solveGS(prec->gs,x,b);
171 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 break;
195 case PASO_AMG:
196 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 }
210
211 #pragma omp parallel for private(i) schedule(static)
212 for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
213 b2[i]=b[2*i+1];
214 }
215
216 Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
217 Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
218
219 #pragma omp parallel for private(i) schedule(static)
220 for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
221 x[2*i]=x1[i];
222 }
223
224 #pragma omp parallel for private(i) schedule(static)
225 for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
226 x[2*i+1]=x2[i];
227 }
228
229 MEMFREE(x1);
230 MEMFREE(x2);
231 MEMFREE(b1);
232 MEMFREE(b2);
233
234 }
235 else if (A->row_block_size==3) {
236 dim_t i;
237
238 double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
239 double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
240 double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
241 double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
242 double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
243 double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
244
245 #pragma omp parallel for private(i) schedule(static)
246 for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
247 b1[i]=b[3*i];
248 }
249
250 #pragma omp parallel for private(i) schedule(static)
251 for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
252 b2[i]=b[3*i+1];
253 }
254
255 #pragma omp parallel for private(i) schedule(static)
256 for (i=0;i<prec->amgSystem->amgblock3->n;i++) {
257 b3[i]=b[3*i+2];
258 }
259
260
261 Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
262 Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
263 Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);
264
265 #pragma omp parallel for private(i) schedule(static)
266 for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
267 x[3*i]=x1[i];
268 }
269
270 #pragma omp parallel for private(i) schedule(static)
271 for (i=0;i<prec->amgSystem->amgblock2->n;i++) {
272 x[3*i+1]=x2[i];
273 }
274
275 #pragma omp parallel for private(i) schedule(static)
276 for (i=0;i<prec->amgSystem->amgblock3->n;i++) {
277 x[3*i+2]=x3[i];
278 }
279
280 MEMFREE(x1);
281 MEMFREE(x2);
282 MEMFREE(x3);
283 MEMFREE(b1);
284 MEMFREE(b2);
285 MEMFREE(b3);
286 }
287 break;
288 }
289
290 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26