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