1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2008 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: AMG preconditioner with reordering */ |
18 |
|
19 |
/**************************************************************/ |
20 |
|
21 |
/* Author: artak@access.edu.au */ |
22 |
|
23 |
/**************************************************************/ |
24 |
|
25 |
#include "Paso.h" |
26 |
#include "Solver.h" |
27 |
#include "PasoUtil.h" |
28 |
#include "UMFPACK.h" |
29 |
#include "Pattern_coupling.h" |
30 |
|
31 |
/**************************************************************/ |
32 |
|
33 |
/* free all memory used by AMG */ |
34 |
|
35 |
void Paso_Solver_AMG_free(Paso_Solver_AMG * in) { |
36 |
if (in!=NULL) { |
37 |
Paso_Solver_AMG_free(in->AMG_of_Schur); |
38 |
Paso_Solver_GS_free(in->GS); |
39 |
MEMFREE(in->inv_A_FF); |
40 |
MEMFREE(in->A_FF_pivot); |
41 |
Paso_SparseMatrix_free(in->A_FC); |
42 |
Paso_SparseMatrix_free(in->A_CF); |
43 |
MEMFREE(in->rows_in_F); |
44 |
MEMFREE(in->rows_in_C); |
45 |
MEMFREE(in->mask_F); |
46 |
MEMFREE(in->mask_C); |
47 |
MEMFREE(in->x_F); |
48 |
MEMFREE(in->b_F); |
49 |
MEMFREE(in->x_C); |
50 |
MEMFREE(in->b_C); |
51 |
MEMFREE(in); |
52 |
} |
53 |
} |
54 |
|
55 |
/**************************************************************/ |
56 |
|
57 |
/* constructs the block-block factorization of |
58 |
|
59 |
[ A_FF A_FC ] |
60 |
A_p= |
61 |
[ A_CF A_FF ] |
62 |
|
63 |
to |
64 |
|
65 |
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] |
66 |
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] |
67 |
|
68 |
|
69 |
where S=A_FF-ACF*invA_FF*A_FC within the shape of S |
70 |
|
71 |
then AMG is applied to S again until S becomes empty |
72 |
|
73 |
*/ |
74 |
Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level) { |
75 |
Paso_Solver_AMG* out=NULL; |
76 |
dim_t n=A_p->numRows; |
77 |
dim_t n_block=A_p->row_block_size; |
78 |
index_t* mis_marker=NULL; |
79 |
index_t* counter=NULL; |
80 |
index_t iPtr,*index, *where_p, iPtr_s; |
81 |
dim_t i,k,j,j0; |
82 |
Paso_SparseMatrix * schur=NULL; |
83 |
Paso_SparseMatrix * schur_withFillIn=NULL; |
84 |
double time0,time1,time2,S; |
85 |
|
86 |
/* identify independend set of rows/columns */ |
87 |
mis_marker=TMPMEMALLOC(n,index_t); |
88 |
counter=TMPMEMALLOC(n,index_t); |
89 |
out=MEMALLOC(1,Paso_Solver_AMG); |
90 |
out->AMG_of_Schur=NULL; |
91 |
out->inv_A_FF=NULL; |
92 |
out->A_FF_pivot=NULL; |
93 |
out->A_FC=NULL; |
94 |
out->A_CF=NULL; |
95 |
out->rows_in_F=NULL; |
96 |
out->rows_in_C=NULL; |
97 |
out->mask_F=NULL; |
98 |
out->mask_C=NULL; |
99 |
out->x_F=NULL; |
100 |
out->b_F=NULL; |
101 |
out->x_C=NULL; |
102 |
out->b_C=NULL; |
103 |
out->A=Paso_SparseMatrix_getReference(A_p); |
104 |
out->GS=Paso_Solver_getGS(A_p,verbose); |
105 |
out->GS->sweeps=2; |
106 |
out->level=level; |
107 |
|
108 |
if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) { |
109 |
/* identify independend set of rows/columns */ |
110 |
time0=Paso_timer(); |
111 |
#pragma omp parallel for private(i) schedule(static) |
112 |
for (i=0;i<n;++i) mis_marker[i]=-1; |
113 |
Paso_Pattern_RS(A_p,mis_marker,1/n); |
114 |
/*Paso_Pattern_coup(A_p,mis_marker,1/n);*/ |
115 |
time2=Paso_timer()-time0; |
116 |
if (Paso_noError()) { |
117 |
#pragma omp parallel for private(i) schedule(static) |
118 |
for (i = 0; i < n; ++i) counter[i]=mis_marker[i]; |
119 |
out->n=n; |
120 |
out->n_block=n_block; |
121 |
out->n_F=Paso_Util_cumsum(n,counter); |
122 |
out->mask_F=MEMALLOC(n,index_t); |
123 |
out->rows_in_F=MEMALLOC(out->n_F,index_t); |
124 |
out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double); |
125 |
out->A_FF_pivot=NULL; /* later use for block size>3 */ |
126 |
if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) { |
127 |
#pragma omp parallel |
128 |
{ |
129 |
/* creates an index for F from mask */ |
130 |
#pragma omp for private(i) schedule(static) |
131 |
for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1; |
132 |
#pragma omp for private(i) schedule(static) |
133 |
for (i = 0; i < n; ++i) { |
134 |
if (mis_marker[i]) { |
135 |
out->rows_in_F[counter[i]]=i; |
136 |
out->mask_F[i]=counter[i]; |
137 |
} else { |
138 |
out->mask_F[i]=-1; |
139 |
} |
140 |
} |
141 |
/* Compute row-sum for getting rs(A_FF)*/ |
142 |
#pragma omp for private(i,iPtr) schedule(static) |
143 |
for (i = 0; i < out->n_F; ++i) { |
144 |
out->inv_A_FF[i]=0; |
145 |
for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) { |
146 |
out->inv_A_FF[i]+=A_p->val[iPtr]; |
147 |
} |
148 |
} |
149 |
|
150 |
#pragma omp for private(i, where_p,iPtr,index) schedule(static) |
151 |
for (i = 0; i < out->n_F; i++) { |
152 |
/* find main diagonal */ |
153 |
iPtr=A_p->pattern->ptr[out->rows_in_F[i]]; |
154 |
index=&(A_p->pattern->index[iPtr]); |
155 |
where_p=(index_t*)bsearch(&out->rows_in_F[i], |
156 |
index, |
157 |
A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]], |
158 |
sizeof(index_t), |
159 |
Paso_comparIndex); |
160 |
if (where_p==NULL) { |
161 |
Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing."); |
162 |
} else { |
163 |
iPtr+=(index_t)(where_p-index); |
164 |
/* get inverse of A_FF block: */ |
165 |
S=out->inv_A_FF[i]; |
166 |
if (ABS(A_p->val[iPtr])>0.) { |
167 |
if(ABS(S)>0.) |
168 |
out->inv_A_FF[i]=1./S; |
169 |
/* else |
170 |
{ |
171 |
fprintf(stderr,"ROWSUM OF ROW %d is 0\n",out->rows_in_F[i]); |
172 |
S=0; |
173 |
for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) { |
174 |
if(A_p->val[iPtr]!=0) { |
175 |
fprintf(stderr,"A_p[%d,%d]=%f ",out->rows_in_F[i],A_p->pattern->index[iPtr],A_p->val[iPtr]); |
176 |
S+=A_p->val[iPtr]; |
177 |
} |
178 |
} |
179 |
fprintf(stderr,"\n SUMMMMMMM %f\n",S); |
180 |
} |
181 |
*/ } else { |
182 |
Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getAMG: Break-down in AMG decomposition: non-regular main diagonal block."); |
183 |
} |
184 |
} |
185 |
} |
186 |
} /* end parallel region */ |
187 |
|
188 |
if( Paso_noError()) { |
189 |
/* if there are no nodes in the coarse level there is no more work to do */ |
190 |
out->n_C=n-out->n_F; |
191 |
/*if (level>0) {*/ |
192 |
if (out->n_C>10) { |
193 |
out->rows_in_C=MEMALLOC(out->n_C,index_t); |
194 |
out->mask_C=MEMALLOC(n,index_t); |
195 |
if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) { |
196 |
/* creates an index for C from mask */ |
197 |
#pragma omp parallel for private(i) schedule(static) |
198 |
for (i = 0; i < n; ++i) counter[i]=! mis_marker[i]; |
199 |
Paso_Util_cumsum(n,counter); |
200 |
#pragma omp parallel |
201 |
{ |
202 |
#pragma omp for private(i) schedule(static) |
203 |
for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1; |
204 |
#pragma omp for private(i) schedule(static) |
205 |
for (i = 0; i < n; ++i) { |
206 |
if (! mis_marker[i]) { |
207 |
out->rows_in_C[counter[i]]=i; |
208 |
out->mask_C[i]=counter[i]; |
209 |
} else { |
210 |
out->mask_C[i]=-1; |
211 |
} |
212 |
} |
213 |
} /* end parallel region */ |
214 |
/* get A_CF block: */ |
215 |
out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F); |
216 |
if (Paso_noError()) { |
217 |
/* get A_FC block: */ |
218 |
out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C); |
219 |
/* get A_CC block: */ |
220 |
if (Paso_noError()) { |
221 |
schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C); |
222 |
|
223 |
/*find the pattern of the schur complement with fill in*/ |
224 |
schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1); |
225 |
|
226 |
/* copy values over*/ |
227 |
#pragma omp for private(i,iPtr,iPtr_s,j,j0) schedule(static) |
228 |
for (i = 0; i < schur_withFillIn->numRows; ++i) { |
229 |
for (iPtr=schur_withFillIn->pattern->ptr[i];iPtr<schur_withFillIn->pattern->ptr[i + 1]; ++iPtr) { |
230 |
j=schur_withFillIn->pattern->index[iPtr]; |
231 |
schur_withFillIn->val[iPtr]=0.; |
232 |
for (iPtr_s=schur->pattern->ptr[i];iPtr_s<schur->pattern->ptr[i + 1]; ++iPtr_s){ |
233 |
j0=schur->pattern->index[iPtr_s]; |
234 |
if (j==j0) { |
235 |
schur_withFillIn->val[iPtr]=schur->val[iPtr_s]; |
236 |
break; |
237 |
} |
238 |
} |
239 |
} |
240 |
} |
241 |
time0=Paso_timer()-time0; |
242 |
|
243 |
if (Paso_noError()) { |
244 |
time1=Paso_timer(); |
245 |
/* update A_CC block to get Schur complement and then apply AMG to it */ |
246 |
Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC); |
247 |
time1=Paso_timer()-time1; |
248 |
out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level+1); |
249 |
|
250 |
/*Paso_SparseMatrix_free(schur);*/ |
251 |
/* Paso_SparseMatrix_free(schur_withFillIn);*/ |
252 |
} |
253 |
/* allocate work arrays for AMG application */ |
254 |
if (Paso_noError()) { |
255 |
out->x_F=MEMALLOC(n_block*out->n_F,double); |
256 |
out->b_F=MEMALLOC(n_block*out->n_F,double); |
257 |
out->x_C=MEMALLOC(n_block*out->n_C,double); |
258 |
out->b_C=MEMALLOC(n_block*out->n_C,double); |
259 |
|
260 |
if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) { |
261 |
#pragma omp parallel |
262 |
{ |
263 |
#pragma omp for private(i,k) schedule(static) |
264 |
for (i = 0; i < out->n_F; ++i) { |
265 |
for (k=0; k<n_block;++k) { |
266 |
out->x_F[i*n_block+k]=0.; |
267 |
out->b_F[i*n_block+k]=0.; |
268 |
} |
269 |
} |
270 |
#pragma omp for private(i,k) schedule(static) |
271 |
for (i = 0; i < out->n_C; ++i) { |
272 |
for (k=0; k<n_block;++k) { |
273 |
out->x_C[i*n_block+k]=0.; |
274 |
out->b_C[i*n_block+k]=0.; |
275 |
} |
276 |
} |
277 |
} /* end parallel region */ |
278 |
} |
279 |
} |
280 |
} |
281 |
} |
282 |
} |
283 |
} |
284 |
} |
285 |
} |
286 |
} |
287 |
} |
288 |
TMPMEMFREE(mis_marker); |
289 |
TMPMEMFREE(counter); |
290 |
if (Paso_noError()) { |
291 |
if (verbose) { |
292 |
printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F); |
293 |
if (out->n_C>0) { |
294 |
printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1); |
295 |
} else { |
296 |
printf("timing: AMG: MIS: %e\n",time2); |
297 |
} |
298 |
} |
299 |
return out; |
300 |
} else { |
301 |
Paso_Solver_AMG_free(out); |
302 |
return NULL; |
303 |
} |
304 |
} |
305 |
|
306 |
/**************************************************************/ |
307 |
|
308 |
/* apply AMG precondition b-> x |
309 |
|
310 |
in fact it solves |
311 |
|
312 |
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] [ x_F ] = [b_F] |
313 |
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C] |
314 |
|
315 |
in the form |
316 |
|
317 |
b->[b_F,b_C] |
318 |
x_F=invA_FF*b_F |
319 |
b_C=b_C-A_CF*x_F |
320 |
x_C=AMG(b_C) |
321 |
b_F=b_F-A_FC*x_C |
322 |
x_F=invA_FF*b_F |
323 |
x<-[x_F,x_C] |
324 |
|
325 |
should be called within a parallel region |
326 |
barrier synconization should be performed to make sure that the input vector available |
327 |
|
328 |
*/ |
329 |
|
330 |
void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) { |
331 |
dim_t i,k,oldsweeps; |
332 |
dim_t n_block=amg->n_block; |
333 |
double *r=MEMALLOC(amg->n,double); |
334 |
/*Paso_Solver_GS* GS=NULL;*/ |
335 |
double *bold=MEMALLOC(amg->n*amg->n_block,double); |
336 |
double *bnew=MEMALLOC(amg->n*amg->n_block,double); |
337 |
|
338 |
/*if (amg->level==0) {*/ |
339 |
if (amg->n_C<=10) { |
340 |
Paso_UMFPACK1(amg->A,x,b,0); |
341 |
} else { |
342 |
/* presmoothing */ |
343 |
Paso_Solver_solveGS(amg->GS,x,b); |
344 |
oldsweeps=amg->GS->sweeps; |
345 |
if (amg->GS->sweeps>1) { |
346 |
|
347 |
#pragma omp parallel for private(i) schedule(static) |
348 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=b[i]; |
349 |
|
350 |
while(amg->GS->sweeps>1) { |
351 |
#pragma omp parallel for private(i) schedule(static) |
352 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bnew[i]=bold[i]+b[i]; |
353 |
/* Compute the residual b=b-Ax*/ |
354 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew); |
355 |
/* Go round again*/ |
356 |
Paso_Solver_solveGS(amg->GS,x,bnew); |
357 |
#pragma omp parallel for private(i) schedule(static) |
358 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=bnew[i]; |
359 |
amg->GS->sweeps=amg->GS->sweeps-1; |
360 |
} |
361 |
} |
362 |
amg->GS->sweeps=oldsweeps; |
363 |
/* end of presmoothing */ |
364 |
|
365 |
#pragma omp parallel for private(i) schedule(static) |
366 |
for (i=0;i<amg->n;++i) r[i]=b[i]; |
367 |
|
368 |
/*r=b-Ax*/ |
369 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
370 |
|
371 |
/* b->[b_F,b_C] */ |
372 |
if (n_block==1) { |
373 |
#pragma omp parallel for private(i) schedule(static) |
374 |
for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]]; |
375 |
#pragma omp parallel for private(i) schedule(static) |
376 |
for (i=0;i<amg->n_C;++i) amg->b_C[i]=r[amg->rows_in_C[i]]; |
377 |
} else { |
378 |
#pragma omp parallel for private(i,k) schedule(static) |
379 |
for (i=0;i<amg->n_F;++i) |
380 |
for (k=0;k<n_block;k++) amg->b_F[amg->n_block*i+k]=r[n_block*amg->rows_in_F[i]+k]; |
381 |
#pragma omp parallel for private(i,k) schedule(static) |
382 |
for (i=0;i<amg->n_C;++i) |
383 |
for (k=0;k<n_block;k++) amg->b_C[amg->n_block*i+k]=r[n_block*amg->rows_in_C[i]+k]; |
384 |
} |
385 |
|
386 |
/* x_F=invA_FF*b_F */ |
387 |
Paso_Solver_applyBlockDiagonalMatrix(n_block,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F); |
388 |
|
389 |
/* b_C=b_C-A_CF*x_F */ |
390 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_CF,amg->x_F,1.,amg->b_C); |
391 |
|
392 |
/* x_C=AMG(b_C) */ |
393 |
Paso_Solver_solveAMG(amg->AMG_of_Schur,amg->x_C,amg->b_C); |
394 |
|
395 |
/* b_F=b_F-A_FC*x_C */ |
396 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F); |
397 |
/* x_F=invA_FF*b_F */ |
398 |
Paso_Solver_applyBlockDiagonalMatrix(n_block,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F); |
399 |
/* x<-[x_F,x_C] */ |
400 |
|
401 |
if (n_block==1) { |
402 |
#pragma omp parallel for private(i) schedule(static) |
403 |
for (i=0;i<amg->n;++i) { |
404 |
if (amg->mask_C[i]>-1) { |
405 |
x[i]+=amg->x_C[amg->mask_C[i]]; |
406 |
} else { |
407 |
x[i]+=amg->x_F[amg->mask_F[i]]; |
408 |
} |
409 |
} |
410 |
} else { |
411 |
#pragma omp parallel for private(i,k) schedule(static) |
412 |
for (i=0;i<amg->n;++i) { |
413 |
if (amg->mask_C[i]>-1) { |
414 |
for (k=0;k<n_block;k++) x[n_block*i+k]+=amg->x_C[n_block*amg->mask_C[i]+k]; |
415 |
} else { |
416 |
for (k=0;k<n_block;k++) x[n_block*i+k]+=amg->x_F[n_block*amg->mask_F[i]+k]; |
417 |
} |
418 |
} |
419 |
} |
420 |
|
421 |
/*postsmoothing*/ |
422 |
Paso_Solver_solveGS1(amg->GS,x,b); |
423 |
if (amg->GS->sweeps>1) { |
424 |
#pragma omp parallel for private(i) schedule(static) |
425 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=b[i]; |
426 |
|
427 |
while(amg->GS->sweeps>1) { |
428 |
#pragma omp parallel for private(i) schedule(static) |
429 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bnew[i]=bold[i]+b[i]; |
430 |
/* Compute the residual b=b-Ax*/ |
431 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew); |
432 |
/* Go round again*/ |
433 |
Paso_Solver_solveGS1(amg->GS,x,bnew); |
434 |
#pragma omp parallel for private(i) schedule(static) |
435 |
for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=bnew[i]; |
436 |
amg->GS->sweeps=amg->GS->sweeps-1; |
437 |
} |
438 |
} |
439 |
/*end of postsmoothing*/ |
440 |
|
441 |
/*Paso_Solver_GS_free(GS);*/ |
442 |
} |
443 |
MEMFREE(r); |
444 |
MEMFREE(bold); |
445 |
MEMFREE(bnew); |
446 |
return; |
447 |
} |
448 |
|
449 |
/* |
450 |
* $Log$ |
451 |
* |
452 |
*/ |