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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26