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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3098 - (show annotations)
Fri Aug 20 08:08:27 2010 UTC (9 years, 1 month ago) by gross
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 9853 byte(s)
fix to the GaussSeidel
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 #include "PasoUtil.h"
30
31 /***********************************************************************************/
32
33 /* free space */
34
35 void Paso_Preconditioner_free(Paso_Solver_Preconditioner* in) {
36 if (in!=NULL) {
37 Paso_Solver_ILU_free(in->ilu);
38 Paso_Solver_RILU_free(in->rilu);
39 Paso_Solver_Jacobi_free(in->jacobi);
40 Paso_Solver_GS_free(in->gs);
41 Paso_Solver_AMG_free(in->amg);
42 Paso_Solver_AMG_System_free(in->amgSystem);
43 Paso_Solver_AMLI_free(in->amli);
44 Paso_Solver_AMLI_System_free(in->amliSystem);
45 MEMFREE(in);
46 }
47 }
48
49
50 /* call the iterative solver: */
51
52 void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
53
54 Paso_Solver_Preconditioner* prec=NULL;
55 dim_t i;
56 /*char filename[7];*/
57 if (A->solver==NULL) {
58 /* allocate structure to hold preconditioner */
59 prec=MEMALLOC(1,Paso_Solver_Preconditioner);
60
61 if (Paso_checkPtr(prec)) return;
62
63 prec->type=UNKNOWN;
64 prec->rilu=NULL;
65 prec->ilu=NULL;
66 prec->jacobi=NULL;
67 prec->gs=NULL;
68 prec->amg=NULL;
69 prec->amli=NULL;
70 prec->amgSystem=NULL;
71 prec->amliSystem=NULL;
72 if (options->verbose && options->use_local_preconditioner) printf("Apply preconditioner locally only.\n");
73
74 A->solver=prec;
75 switch (options->preconditioner) {
76 default:
77 case PASO_JACOBI:
78 if (options->verbose) printf("Jacobi preconditioner is used.\n");
79 prec->jacobi=Paso_Solver_getJacobi(A);
80 prec->type=PASO_JACOBI;
81 break;
82 case PASO_GS:
83 if (options->verbose) printf("Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
84 prec->gs=Paso_Solver_getGS(A, options->sweeps, options->use_local_preconditioner, options->verbose);
85 prec->type=PASO_GS;
86 break;
87
88 /***************************************************************************************/
89 case PASO_ILU0:
90 if (options->verbose) printf("ILU preconditioner is used.\n");
91 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
92 prec->type=PASO_ILU0;
93 Paso_MPIInfo_noError(A->mpi_info);
94 break;
95 case PASO_RILU:
96 if (options->verbose) printf("RILU preconditioner is used.\n");
97 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
98 Paso_MPIInfo_noError(A->mpi_info);
99 prec->type=PASO_RILU;
100 break;
101
102 case PASO_AMG:
103 if (options->verbose) printf("AMG preconditioner is used.\n");
104 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
105 Paso_MPIInfo_noError(A->mpi_info);
106 /*
107 prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
108 if (Paso_checkPtr(prec->amgSystem)) return;
109 prec->amgSystem->block_size=A->row_block_size;
110 for (i=0;i<A->row_block_size;++i) {
111 prec->amgSystem->amgblock[i]=NULL;
112 prec->amgSystem->block[i]=NULL;
113 }
114 */
115
116 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
117 /*if (A->row_block_size==1) {
118 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
119 }
120 else {
121 for (i=0;i<A->row_block_size;++i) {
122 prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
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
133 prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);
134 if (Paso_checkPtr(prec->amliSystem)) return;
135
136 prec->amliSystem->block_size=A->row_block_size;
137
138 for (i=0;i<A->row_block_size;++i) {
139 prec->amliSystem->amliblock[i]=NULL;
140 prec->amliSystem->block[i]=NULL;
141 }
142
143 if (options->verbose) printf("AMLI preconditioner is used.\n");
144
145 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
146 if (A->row_block_size==1) {
147 prec->amli=Paso_Solver_getAMLI(A->mainBlock,options->level_max,options);
148 }
149 else {
150 for (i=0;i<A->row_block_size;++i) {
151 prec->amliSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
152 prec->amliSystem->amliblock[i]=Paso_Solver_getAMLI(prec->amliSystem->block[i],options->level_max,options);
153 }
154 }
155 prec->type=PASO_AMLI;
156 break;
157
158 }
159 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
160 Paso_Preconditioner_free(prec);
161 A->solver=NULL;
162 }
163 }
164 }
165
166 /* applies the preconditioner */
167 /* has to be called within a parallel reqion */
168 /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
169 void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
170 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
171 dim_t i,j;
172 dim_t n=A->mainBlock->numRows;
173
174
175
176
177 switch (prec->type) {
178 default:
179 case PASO_JACOBI:
180 Paso_Solver_solveJacobi(A, prec->jacobi,x,b);
181 break;
182
183 case PASO_GS:
184 Paso_Solver_solveGS(A, prec->gs,x,b);
185 break;
186 case PASO_ILU0:
187 Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
188 break;
189 case PASO_RILU:
190 Paso_Solver_solveRILU(prec->rilu,x,b);
191 break;
192
193 case PASO_AMG:
194 Paso_Solver_solveAMG(prec->amg,x,b);
195
196 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
197 /*
198 if (A->row_block_size<4) {
199 Paso_Solver_solveAMG(prec->amg,x,b);
200 }
201 else {
202
203
204 for (i=0;i<A->row_block_size;i++) {
205 xx[i]=MEMALLOC(n,double);
206 bb[i]=MEMALLOC(n,double);
207 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
208 }
209
210 for (i=0;i<n;i++) {
211 for (j=0;j<A->row_block_size;j++) {
212 bb[j][i]=b[A->row_block_size*i+j];
213 }
214 }
215
216
217 for (i=0;i<A->row_block_size;i++) {
218 Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
219 }
220
221 for (i=0;i<n;i++) {
222 for (j=0;j<A->row_block_size;j++) {
223 x[A->row_block_size*i+j]=xx[j][i];
224 }
225 }
226
227 for (i=0;i<A->row_block_size;i++) {
228 MEMFREE(xx[i]);
229 MEMFREE(bb[i]);
230 }
231 }
232 */
233 break;
234 case PASO_AMLI:
235
236 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
237 if (A->row_block_size==1) {
238 Paso_Solver_solveAMLI(prec->amli,x,b);
239 }
240 else {
241 double **xx;
242 double **bb;
243 xx=MEMALLOC(A->row_block_size,double*);
244 bb=MEMALLOC(A->row_block_size,double*);
245 if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
246 for (i=0;i<A->row_block_size;i++) {
247 xx[i]=MEMALLOC(n,double);
248 bb[i]=MEMALLOC(n,double);
249 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
250 }
251
252 /*#pragma omp parallel for private(i,j) schedule(static)*/
253 for (i=0;i<n;i++) {
254 for (j=0;j<A->row_block_size;j++) {
255 bb[j][i]=b[A->row_block_size*i+j];
256 }
257 }
258
259
260 for (i=0;i<A->row_block_size;i++) {
261 Paso_Solver_solveAMLI(prec->amliSystem->amliblock[i],xx[i],bb[i]);
262 }
263
264 /*#pragma omp parallel for private(i,j) schedule(static)*/
265 for (i=0;i<n;i++) {
266 for (j=0;j<A->row_block_size;j++) {
267 x[A->row_block_size*i+j]=xx[j][i];
268 }
269 }
270
271 for (i=0;i<A->row_block_size;i++) {
272 MEMFREE(xx[i]);
273 MEMFREE(bb[i]);
274 }
275 MEMFREE(xx);
276 MEMFREE(bb);
277 }
278 break;
279 }
280
281 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26