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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3094 - (hide annotations)
Fri Aug 13 08:38:06 2010 UTC (9 years, 1 month ago) by gross
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 9833 byte(s)
The MPI and sequational GAUSS_SEIDEL have been merged.
The couring and main diagonal pointer is now manged by the patternm which means that they are calculated once only even if the preconditioner is deleted.



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 lgao 3051 #include "PasoUtil.h"
30 jgs 150
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 gross 430 Paso_Solver_RILU_free(in->rilu);
39 jgs 150 Paso_Solver_Jacobi_free(in->jacobi);
40 gross 3094 Paso_Solver_GS_free(in->gs);
41 artak 1842 Paso_Solver_AMG_free(in->amg);
42 artak 2659 Paso_Solver_AMG_System_free(in->amgSystem);
43 artak 2760 Paso_Solver_AMLI_free(in->amli);
44     Paso_Solver_AMLI_System_free(in->amliSystem);
45 jgs 150 MEMFREE(in);
46     }
47     }
48 gross 3005
49    
50 jgs 150 /* call the iterative solver: */
51    
52     void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
53 artak 2659
54 jgs 150 Paso_Solver_Preconditioner* prec=NULL;
55 artak 2662 dim_t i;
56 artak 2992 /*char filename[7];*/
57 gross 425 if (A->solver==NULL) {
58 jgs 150 /* allocate structure to hold preconditioner */
59     prec=MEMALLOC(1,Paso_Solver_Preconditioner);
60 gross 3005
61 jgs 150 if (Paso_checkPtr(prec)) return;
62 artak 2760
63 jgs 150 prec->type=UNKNOWN;
64 gross 430 prec->rilu=NULL;
65 jgs 150 prec->ilu=NULL;
66     prec->jacobi=NULL;
67 artak 1819 prec->gs=NULL;
68 artak 1842 prec->amg=NULL;
69 artak 2760 prec->amli=NULL;
70 gross 3005 prec->amgSystem=NULL;
71     prec->amliSystem=NULL;
72 gross 3094 if (options->verbose && options->use_local_preconditioner) printf("Apply preconditioner locally only.\n");
73 artak 2662
74 gross 425 A->solver=prec;
75 jgs 150 switch (options->preconditioner) {
76     default:
77     case PASO_JACOBI:
78     if (options->verbose) printf("Jacobi preconditioner is used.\n");
79 gross 3094 prec->jacobi=Paso_Solver_getJacobi(A);
80 jgs 150 prec->type=PASO_JACOBI;
81     break;
82 gross 3094 case PASO_GS:
83     if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
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 jgs 150 case PASO_ILU0:
90     if (options->verbose) printf("ILU preconditioner is used.\n");
91 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
92 jgs 150 prec->type=PASO_ILU0;
93 gross 3094 Paso_MPIInfo_noError(A->mpi_info);
94 jgs 150 break;
95 gross 430 case PASO_RILU:
96     if (options->verbose) printf("RILU preconditioner is used.\n");
97 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
98 gross 3094 Paso_MPIInfo_noError(A->mpi_info);
99 gross 430 prec->type=PASO_RILU;
100     break;
101 gross 3094
102 artak 1842 case PASO_AMG:
103 artak 3013 if (options->verbose) printf("AMG preconditioner is used.\n");
104 gross 3094 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
105     Paso_MPIInfo_noError(A->mpi_info);
106 artak 3013 /*
107     prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
108 gross 3005 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 artak 3013 */
115    
116 artak 2662 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
117 artak 3013 /*if (A->row_block_size==1) {
118 artak 2659 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
119     }
120 artak 2662 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 artak 2659 }
126 artak 3013 */
127 artak 2992
128 artak 1842 prec->type=PASO_AMG;
129     break;
130 artak 2992
131 artak 2760 case PASO_AMLI:
132 artak 3017
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 gross 3005 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 artak 2760
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 artak 1842
158 jgs 150 }
159 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
160 artak 2659 Paso_Preconditioner_free(prec);
161     A->solver=NULL;
162 jgs 150 }
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 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
171 artak 2662 dim_t i,j;
172 artak 2670 dim_t n=A->mainBlock->numRows;
173 gross 3094
174 artak 2992
175 artak 2662
176 gross 3094
177 jgs 150 switch (prec->type) {
178     default:
179     case PASO_JACOBI:
180 gross 3094 Paso_Solver_solveJacobi(A, prec->jacobi,x,b);
181 jgs 150 break;
182 gross 3094
183     case PASO_GS:
184     Paso_Solver_solveGS(A, prec->gs,x,b);
185     break;
186 jgs 150 case PASO_ILU0:
187 gross 3094 Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
188 jgs 150 break;
189 gross 430 case PASO_RILU:
190 gross 3094 Paso_Solver_solveRILU(prec->rilu,x,b);
191 gross 430 break;
192 gross 3094
193 artak 2662 case PASO_AMG:
194 artak 3013 Paso_Solver_solveAMG(prec->amg,x,b);
195 artak 2992
196 artak 2662 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
197 artak 3013 /*
198 artak 3010 if (A->row_block_size<4) {
199 artak 2659 Paso_Solver_solveAMG(prec->amg,x,b);
200     }
201 artak 2662 else {
202    
203 artak 2659
204 artak 2670 for (i=0;i<A->row_block_size;i++) {
205 artak 2662 xx[i]=MEMALLOC(n,double);
206     bb[i]=MEMALLOC(n,double);
207 artak 2670 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
208 artak 2659 }
209 artak 2670
210 artak 2662 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 artak 2660 }
215 artak 2670
216    
217 artak 2662 for (i=0;i<A->row_block_size;i++) {
218     Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
219 artak 2659 }
220 artak 2670
221 artak 2662 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 artak 2670
227 artak 2662 for (i=0;i<A->row_block_size;i++) {
228     MEMFREE(xx[i]);
229     MEMFREE(bb[i]);
230     }
231 artak 2659 }
232 artak 3013 */
233 artak 2659 break;
234 artak 2760 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 gross 3094 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 artak 2760 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 gross 3094 MEMFREE(xx);
276     MEMFREE(bb);
277 artak 2760 }
278     break;
279 jgs 150 }
280 gross 3094
281 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26