/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver_ILU.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver_ILU.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 115 - (show annotations)
Fri Mar 4 07:12:47 2005 UTC (14 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 15868 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: ILU preconditioner with reordering */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003,2004,2005 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Solver.h"
15 #include "Util.h"
16
17 /**************************************************************/
18
19 /* free all memory used by ILU */
20
21 void Finley_Solver_ILU_free(Finley_Solver_ILU * in) {
22 if (in!=NULL) {
23 Finley_Solver_ILU_free(in->ILU_of_Schur);
24 MEMFREE(in->inv_A_FF);
25 MEMFREE(in->A_FF_pivot);
26 Finley_SystemMatrix_dealloc(in->A_FC);
27 Finley_SystemMatrix_dealloc(in->A_CF);
28 MEMFREE(in->rows_in_F);
29 MEMFREE(in->rows_in_C);
30 MEMFREE(in->mask_F);
31 MEMFREE(in->mask_C);
32 MEMFREE(in->x_F);
33 MEMFREE(in->b_F);
34 MEMFREE(in->x_C);
35 MEMFREE(in->b_C);
36 MEMFREE(in);
37 }
38 }
39
40 /**************************************************************/
41
42 /* constructs the block-block factorization of
43
44 [ A_FF A_FC ]
45 A_p=
46 [ A_CF A_FF ]
47
48 to
49
50 [ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ]
51 [ A_CF*invA_FF I ] [ 0 S ] [ 0 I ]
52
53
54 where S=A_FF-ACF*invA_FF*A_FC within the shape of S
55
56 then ILU is applied to S again until S becomes empty
57
58 */
59 Finley_Solver_ILU* Finley_Solver_getILU(Finley_SystemMatrix * A_p,int verbose) {
60 Finley_Solver_ILU* out=NULL;
61 maybelong n=A_p->num_rows;
62 maybelong n_block=A_p->row_block_size;
63 maybelong* mis_marker=NULL;
64 maybelong* counter=NULL;
65 maybelong i,iPtr, *index, *where_p,k;
66 Finley_SystemMatrix * schur=NULL;
67 double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;
68
69
70 /* identify independend set of rows/columns */
71 mis_marker=TMPMEMALLOC(n,maybelong);
72 counter=TMPMEMALLOC(n,maybelong);
73 out=MEMALLOC(1,Finley_Solver_ILU);
74 out->ILU_of_Schur=NULL;
75 out->inv_A_FF=NULL;
76 out->A_FF_pivot=NULL;
77 out->A_FC=NULL;
78 out->A_CF=NULL;
79 out->rows_in_F=NULL;
80 out->rows_in_C=NULL;
81 out->mask_F=NULL;
82 out->mask_C=NULL;
83 out->x_F=NULL;
84 out->b_F=NULL;
85 out->x_C=NULL;
86 out->b_C=NULL;
87
88 if ( !(Finley_checkPtr(mis_marker) || Finley_checkPtr(out) || Finley_checkPtr(counter) ) ) {
89 /* identify independend set of rows/columns */
90 time0=Finley_timer();
91 Finley_SystemMatrixPattern_mis(A_p->pattern,mis_marker);
92 time2=Finley_timer()-time0;
93 if (Finley_ErrorCode==NO_ERROR) {
94 #pragma omp parallel for private(i) schedule(static)
95 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
96 out->n=n;
97 out->n_block=n_block;
98 out->n_F=Finley_Util_cumsum(n,counter);
99 if (verbose) {
100 printf("ILU: number of vertices to be eliminated = %d\n",out->n_F);
101 printf("ILU: number of vertices in coarse level = %d\n",n-out->n_F);
102 }
103 out->mask_F=MEMALLOC(n,maybelong);
104 out->rows_in_F=MEMALLOC(out->n_F,maybelong);
105 out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);
106 out->A_FF_pivot=NULL; /* later use for block size>3 */
107 if (! (Finley_checkPtr(out->mask_F) || Finley_checkPtr(out->inv_A_FF) || Finley_checkPtr(out->rows_in_F) ) ) {
108 #pragma omp parallel
109 {
110 /* creates an index for F from mask */
111 #pragma omp for private(i) schedule(static)
112 for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
113 #pragma omp for private(i) schedule(static)
114 for (i = 0; i < n; ++i) {
115 if (mis_marker[i]) {
116 out->rows_in_F[counter[i]]=i;
117 out->mask_F[i]=counter[i];
118 } else {
119 out->mask_F[i]=-1;
120 }
121 }
122 #pragma omp for private(i, where_p,iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D,index) schedule(static)
123 for (i = 0; i < out->n_F; i++) {
124 /* find main diagonal */
125 iPtr=A_p->pattern->ptr[out->rows_in_F[i]];
126 index=&(A_p->pattern->index[iPtr]);
127 where_p=(maybelong*)bsearch(&out->rows_in_F[i],
128 index,
129 A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]],
130 sizeof(maybelong),
131 Finley_comparIndex);
132 if (where_p==NULL) {
133 Finley_ErrorCode = VALUE_ERROR;
134 sprintf(Finley_ErrorMsg, "no main diagonal in row %d",i);
135 } else {
136 iPtr+=(maybelong)(where_p-index);
137 /* get inverse of A_FF block: */
138 if (n_block==1) {
139 if (ABS(A_p->val[iPtr])>0.) {
140 out->inv_A_FF[i]=1./A_p->val[iPtr];
141 } else {
142 Finley_ErrorCode = ZERO_DIVISION_ERROR;
143 sprintf(Finley_ErrorMsg, "Break-down in ILU decomposition: non-regular main diagonal block in row %d",i);
144 }
145 } else if (n_block==2) {
146 A11=A_p->val[iPtr*4];
147 A12=A_p->val[iPtr*4+2];
148 A21=A_p->val[iPtr*4+1];
149 A22=A_p->val[iPtr*4+3];
150 D = A11*A22-A12*A21;
151 if (ABS(D) > 0 ){
152 D=1./D;
153 out->inv_A_FF[i*4]= A22*D;
154 out->inv_A_FF[i*4+1]=-A21*D;
155 out->inv_A_FF[i*4+2]=-A12*D;
156 out->inv_A_FF[i*4+3]= A11*D;
157 } else {
158 Finley_ErrorCode = ZERO_DIVISION_ERROR;
159 sprintf(Finley_ErrorMsg, "Break-down in ILU decomposition: non-regular main diagonal block in row %d",i);
160 }
161 } else if (n_block==3) {
162 A11=A_p->val[iPtr*9 ];
163 A21=A_p->val[iPtr*9+1];
164 A31=A_p->val[iPtr*9+2];
165 A12=A_p->val[iPtr*9+3];
166 A22=A_p->val[iPtr*9+4];
167 A32=A_p->val[iPtr*9+5];
168 A13=A_p->val[iPtr*9+6];
169 A23=A_p->val[iPtr*9+7];
170 A33=A_p->val[iPtr*9+8];
171 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
172 if (ABS(D) > 0 ){
173 D=1./D;
174 out->inv_A_FF[i*9 ]=(A22*A33-A23*A32)*D;
175 out->inv_A_FF[i*9+1]=(A31*A23-A21*A33)*D;
176 out->inv_A_FF[i*9+2]=(A21*A32-A31*A22)*D;
177 out->inv_A_FF[i*9+3]=(A13*A32-A12*A33)*D;
178 out->inv_A_FF[i*9+4]=(A11*A33-A31*A13)*D;
179 out->inv_A_FF[i*9+5]=(A12*A31-A11*A32)*D;
180 out->inv_A_FF[i*9+6]=(A12*A23-A13*A22)*D;
181 out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;
182 out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;
183 } else {
184 Finley_ErrorCode = ZERO_DIVISION_ERROR;
185 sprintf(Finley_ErrorMsg, "Break-down in ILU decomposition: non-regular main diagonal block in row %d",i);
186 }
187 }
188 }
189 }
190 } /* end parallel region */
191
192 /* if there are no nodes in the coarse level there is no more work to do */
193 out->n_C=n-out->n_F;
194 if (out->n_C>0) {
195 out->rows_in_C=MEMALLOC(out->n_C,maybelong);
196 out->mask_C=MEMALLOC(n,maybelong);
197 if (! (Finley_checkPtr(out->mask_C) || Finley_checkPtr(out->rows_in_C) ) ) {
198 /* creates an index for C from mask */
199 #pragma omp parallel for private(i) schedule(static)
200 for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
201 Finley_Util_cumsum(n,counter);
202 #pragma omp parallel
203 {
204 #pragma omp for private(i) schedule(static)
205 for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;
206 #pragma omp for private(i) schedule(static)
207 for (i = 0; i < n; ++i) {
208 if (! mis_marker[i]) {
209 out->rows_in_C[counter[i]]=i;
210 out->mask_C[i]=counter[i];
211 } else {
212 out->mask_C[i]=-1;
213 }
214 }
215 } /* end parallel region */
216 /* get A_CF block: */
217 out->A_CF=Finley_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_F);
218 if (Finley_ErrorCode==NO_ERROR) {
219 /* get A_FC block: */
220 out->A_FC=Finley_SystemMatrix_getSubmatrix(A_p,out->n_F,out->rows_in_F,out->mask_C);
221 /* get A_FF block: */
222 if (Finley_ErrorCode==NO_ERROR) {
223 schur=Finley_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_C);
224 time0=Finley_timer()-time0;
225 if (Finley_ErrorCode==NO_ERROR) {
226 time1=Finley_timer();
227 /* update A_CC block to get Schur complement and then apply ILU to it */
228 Finley_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
229 time1=Finley_timer()-time1;
230 out->ILU_of_Schur=Finley_Solver_getILU(schur,verbose);
231 Finley_SystemMatrix_dealloc(schur);
232 }
233 /* allocate work arrays for ILU application */
234 if (Finley_ErrorCode==NO_ERROR) {
235 out->x_F=MEMALLOC(n_block*out->n_F,double);
236 out->b_F=MEMALLOC(n_block*out->n_F,double);
237 out->x_C=MEMALLOC(n_block*out->n_C,double);
238 out->b_C=MEMALLOC(n_block*out->n_C,double);
239 if (! (Finley_checkPtr(out->x_F) || Finley_checkPtr(out->b_F) || Finley_checkPtr(out->x_C) || Finley_checkPtr(out->b_C) ) ) {
240 #pragma omp parallel
241 {
242 #pragma omp for private(i,k) schedule(static)
243 for (i = 0; i < out->n_F; ++i) {
244 for (k=0; k<n_block;++k) {
245 out->x_F[i*n_block+k]=0.;
246 out->b_F[i*n_block+k]=0.;
247 }
248 }
249 #pragma omp for private(i,k) schedule(static)
250 for (i = 0; i < out->n_C; ++i) {
251 for (k=0; k<n_block;++k) {
252 out->x_C[i*n_block+k]=0.;
253 out->b_C[i*n_block+k]=0.;
254 }
255 }
256 } /* end parallel region */
257 }
258 }
259 }
260 }
261 }
262 }
263 }
264 }
265 }
266 TMPMEMFREE(mis_marker);
267 TMPMEMFREE(counter);
268 if (Finley_ErrorCode==NO_ERROR) {
269 if (verbose) {
270 /* if (out->n_C>0) {
271 printf("timing: ILU: reordering : %e\n",time0);
272 printf("timing: ILU: elemination : %e\n",time1);
273 }
274 printf("timing: ILU: mis %e\n",time2); */
275 }
276 return out;
277 } else {
278 Finley_Solver_ILU_free(out);
279 return NULL;
280 }
281 }
282
283 /**************************************************************/
284
285 /* apply ILU precondition b-> x
286
287 in fact it solves
288
289 [ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] [ x_F ] = [b_F]
290 [ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C]
291
292 in the form
293
294 b->[b_F,b_C]
295 x_F=invA_FF*b_F
296 b_C=b_C-A_CF*x_F
297 x_C=ILU(b_C)
298 b_F=b_F-A_FC*x_C
299 x_F=invA_FF*b_F
300 x<-[x_F,x_C]
301
302 should be called within a parallel region
303 barrier synconization should be performed to make sure that the input vector available
304
305 */
306
307 void Finley_Solver_solveILU(Finley_Solver_ILU * ilu, double * x, double * b) {
308 maybelong i,k;
309 maybelong n_block=ilu->n_block;
310
311 if (ilu->n_C==0) {
312 /* x=invA_FF*b */
313 Finley_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,x,b);
314 } else {
315 /* b->[b_F,b_C] */
316 if (n_block==1) {
317 #pragma omp for private(i) schedule(static)
318 for (i=0;i<ilu->n_F;++i) ilu->b_F[i]=b[ilu->rows_in_F[i]];
319 #pragma omp for private(i) schedule(static)
320 for (i=0;i<ilu->n_C;++i) ilu->b_C[i]=b[ilu->rows_in_C[i]];
321 } else {
322 #pragma omp for private(i,k) schedule(static)
323 for (i=0;i<ilu->n_F;++i)
324 for (k=0;k<n_block;k++) ilu->b_F[ilu->n_block*i+k]=b[n_block*ilu->rows_in_F[i]+k];
325 #pragma omp for private(i,k) schedule(static)
326 for (i=0;i<ilu->n_C;++i)
327 for (k=0;k<n_block;k++) ilu->b_C[ilu->n_block*i+k]=b[n_block*ilu->rows_in_C[i]+k];
328 }
329 /* x_F=invA_FF*b_F */
330 Finley_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,ilu->x_F,ilu->b_F);
331 /* b_C=b_C-A_CF*x_F */
332 Finley_RawScaledSystemMatrixVector(-1.,ilu->A_CF,ilu->x_F,1.,ilu->b_C);
333 /* x_C=ILU(b_C) */
334 Finley_Solver_solveILU(ilu->ILU_of_Schur,ilu->x_C,ilu->b_C);
335 /* b_F=b_F-A_FC*x_C */
336 Finley_RawScaledSystemMatrixVector(-1.,ilu->A_FC,ilu->x_C,1.,ilu->b_F);
337 /* x_F=invA_FF*b_F */
338 Finley_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,ilu->x_F,ilu->b_F);
339 /* x<-[x_F,x_C] */
340 if (n_block==1) {
341 #pragma omp for private(i) schedule(static)
342 for (i=0;i<ilu->n;++i) {
343 if (ilu->mask_C[i]>-1) {
344 x[i]=ilu->x_C[ilu->mask_C[i]];
345 } else {
346 x[i]=ilu->x_F[ilu->mask_F[i]];
347 }
348 }
349 } else {
350 #pragma omp for private(i,k) schedule(static)
351 for (i=0;i<ilu->n;++i) {
352 if (ilu->mask_C[i]>-1) {
353 for (k=0;k<n_block;k++) x[n_block*i+k]=ilu->x_C[n_block*ilu->mask_C[i]+k];
354 } else {
355 for (k=0;k<n_block;k++) x[n_block*i+k]=ilu->x_F[n_block*ilu->mask_F[i]+k];
356 }
357 }
358 }
359 /* all done */
360 }
361 return;
362 }
363 /*
364 * $Log$
365 * Revision 1.2 2005/03/04 07:12:47 jgs
366 * *** empty log message ***
367 *
368 * Revision 1.1.2.3 2005/03/04 05:27:08 gross
369 * bug in SystemPattern fixed.
370 *
371 * Revision 1.1.2.2 2005/03/04 05:09:46 gross
372 * a missing file from the ILU reimplementation
373 *
374 * Revision 1.1.2.1 2005/03/02 23:35:07 gross
375 * reimplementation of the ILU in Finley. block size>1 still needs some testing
376 *
377 *
378 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26