1 |
jfenwick |
1851 |
|
2 |
|
|
/******************************************************* |
3 |
|
|
* |
4 |
jfenwick |
2548 |
* Copyright (c) 2003-2009 by University of Queensland |
5 |
jfenwick |
1851 |
* 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 |
artak |
2475 |
/* Paso: AMG preconditioner */ |
18 |
jfenwick |
1851 |
|
19 |
|
|
/**************************************************************/ |
20 |
|
|
|
21 |
jfenwick |
2625 |
/* Author: artak@uq.edu.au */ |
22 |
jfenwick |
1851 |
|
23 |
|
|
/**************************************************************/ |
24 |
|
|
|
25 |
|
|
#include "Paso.h" |
26 |
|
|
#include "Solver.h" |
27 |
artak |
2475 |
#include "Options.h" |
28 |
jfenwick |
1851 |
#include "PasoUtil.h" |
29 |
artak |
1931 |
#include "UMFPACK.h" |
30 |
artak |
2507 |
#include "MKL.h" |
31 |
|
|
#include "SystemMatrix.h" |
32 |
phornby |
1917 |
#include "Pattern_coupling.h" |
33 |
jfenwick |
1851 |
|
34 |
|
|
/**************************************************************/ |
35 |
|
|
|
36 |
|
|
/* free all memory used by AMG */ |
37 |
|
|
|
38 |
artak |
2659 |
void Paso_Solver_AMG_System_free(Paso_Solver_AMG_System * in) { |
39 |
artak |
2662 |
dim_t i; |
40 |
artak |
2659 |
if (in!=NULL) { |
41 |
artak |
2662 |
for (i=0;i<in->block_size;++i) { |
42 |
|
|
Paso_Solver_AMG_free(in->amgblock[i]); |
43 |
artak |
2670 |
Paso_SparseMatrix_free(in->block[i]); |
44 |
artak |
2662 |
} |
45 |
artak |
2659 |
MEMFREE(in); |
46 |
|
|
} |
47 |
|
|
} |
48 |
|
|
|
49 |
|
|
|
50 |
|
|
/* free all memory used by AMG */ |
51 |
|
|
|
52 |
jfenwick |
1851 |
void Paso_Solver_AMG_free(Paso_Solver_AMG * in) { |
53 |
|
|
if (in!=NULL) { |
54 |
artak |
2828 |
|
55 |
|
|
if(in->Smoother->ID==PASO_JACOBI) |
56 |
|
|
Paso_Solver_Jacobi_free(in->Smoother->Jacobi); |
57 |
|
|
else if (in->Smoother->ID==PASO_GS) |
58 |
|
|
Paso_Solver_GS_free(in->Smoother->GS); |
59 |
|
|
MEMFREE(in->Smoother); |
60 |
|
|
|
61 |
jfenwick |
1851 |
Paso_SparseMatrix_free(in->A_FC); |
62 |
artak |
2760 |
Paso_SparseMatrix_free(in->A_FF); |
63 |
|
|
Paso_SparseMatrix_free(in->W_FC); |
64 |
jfenwick |
1851 |
Paso_SparseMatrix_free(in->A_CF); |
65 |
artak |
2760 |
Paso_SparseMatrix_free(in->P); |
66 |
|
|
Paso_SparseMatrix_free(in->R); |
67 |
artak |
2556 |
Paso_SparseMatrix_free(in->A); |
68 |
artak |
2711 |
if(in->coarsest_level==TRUE) { |
69 |
|
|
#ifdef MKL |
70 |
|
|
Paso_MKL_free1(in->AOffset1); |
71 |
|
|
Paso_SparseMatrix_free(in->AOffset1); |
72 |
|
|
#else |
73 |
|
|
#ifdef UMFPACK |
74 |
|
|
Paso_UMFPACK1_free((Paso_UMFPACK_Handler*)(in->solver)); |
75 |
|
|
#endif |
76 |
|
|
#endif |
77 |
|
|
} |
78 |
jfenwick |
1851 |
MEMFREE(in->rows_in_F); |
79 |
|
|
MEMFREE(in->rows_in_C); |
80 |
|
|
MEMFREE(in->mask_F); |
81 |
|
|
MEMFREE(in->mask_C); |
82 |
|
|
MEMFREE(in->x_F); |
83 |
|
|
MEMFREE(in->b_F); |
84 |
|
|
MEMFREE(in->x_C); |
85 |
|
|
MEMFREE(in->b_C); |
86 |
gross |
2551 |
in->solver=NULL; |
87 |
artak |
2760 |
Paso_Solver_AMG_free(in->AMG_of_Coarse); |
88 |
gross |
2551 |
MEMFREE(in->b_C); |
89 |
jfenwick |
1851 |
MEMFREE(in); |
90 |
|
|
} |
91 |
|
|
} |
92 |
|
|
|
93 |
|
|
/**************************************************************/ |
94 |
|
|
|
95 |
|
|
/* constructs the block-block factorization of |
96 |
|
|
|
97 |
|
|
[ A_FF A_FC ] |
98 |
|
|
A_p= |
99 |
|
|
[ A_CF A_FF ] |
100 |
|
|
|
101 |
|
|
to |
102 |
|
|
|
103 |
|
|
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] |
104 |
|
|
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] |
105 |
|
|
|
106 |
|
|
|
107 |
|
|
where S=A_FF-ACF*invA_FF*A_FC within the shape of S |
108 |
|
|
|
109 |
|
|
then AMG is applied to S again until S becomes empty |
110 |
|
|
|
111 |
|
|
*/ |
112 |
artak |
2520 |
Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) { |
113 |
jfenwick |
1851 |
Paso_Solver_AMG* out=NULL; |
114 |
artak |
2760 |
/* |
115 |
|
|
Paso_Pattern* temp1=NULL; |
116 |
artak |
2556 |
Paso_Pattern* temp2=NULL; |
117 |
artak |
2760 |
*/ |
118 |
artak |
2520 |
bool_t verbose=options->verbose; |
119 |
artak |
2831 |
bool_t timing=0; |
120 |
|
|
|
121 |
jfenwick |
1851 |
dim_t n=A_p->numRows; |
122 |
|
|
dim_t n_block=A_p->row_block_size; |
123 |
|
|
index_t* mis_marker=NULL; |
124 |
artak |
1881 |
index_t* counter=NULL; |
125 |
artak |
2760 |
/*index_t iPtr,*index, *where_p;*/ |
126 |
|
|
dim_t i; |
127 |
|
|
Paso_SparseMatrix * A_c=NULL; |
128 |
artak |
2765 |
double time0=0; |
129 |
artak |
2784 |
Paso_SparseMatrix * Atemp=NULL; |
130 |
|
|
double sparsity=0; |
131 |
artak |
2475 |
|
132 |
artak |
2760 |
/* |
133 |
|
|
double *temp,*temp_1; |
134 |
|
|
double S; |
135 |
|
|
index_t iptr; |
136 |
|
|
*/ |
137 |
|
|
/*char filename[8];*/ |
138 |
|
|
/*sprintf(filename,"AMGLevel%d",level); |
139 |
artak |
2726 |
|
140 |
artak |
2760 |
Paso_SparseMatrix_saveMM(A_p,filename); |
141 |
artak |
2726 |
*/ |
142 |
|
|
|
143 |
artak |
2439 |
/*Make sure we have block sizes 1*/ |
144 |
gross |
2550 |
if (A_p->col_block_size>1) { |
145 |
|
|
Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1."); |
146 |
|
|
return NULL; |
147 |
|
|
} |
148 |
|
|
if (A_p->row_block_size>1) { |
149 |
|
|
Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires row block size 1."); |
150 |
|
|
return NULL; |
151 |
|
|
} |
152 |
gross |
2551 |
out=MEMALLOC(1,Paso_Solver_AMG); |
153 |
artak |
2828 |
out->Smoother=MEMALLOC(1,Paso_Solver_Smoother); |
154 |
jfenwick |
1851 |
/* identify independend set of rows/columns */ |
155 |
|
|
mis_marker=TMPMEMALLOC(n,index_t); |
156 |
|
|
counter=TMPMEMALLOC(n,index_t); |
157 |
gross |
2551 |
if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) { |
158 |
artak |
2760 |
out->AMG_of_Coarse=NULL; |
159 |
|
|
out->A_FF=NULL; |
160 |
gross |
2551 |
out->A_FC=NULL; |
161 |
|
|
out->A_CF=NULL; |
162 |
artak |
2760 |
out->W_FC=NULL; |
163 |
|
|
out->P=NULL; |
164 |
|
|
out->R=NULL; |
165 |
gross |
2551 |
out->rows_in_F=NULL; |
166 |
|
|
out->rows_in_C=NULL; |
167 |
|
|
out->mask_F=NULL; |
168 |
|
|
out->mask_C=NULL; |
169 |
|
|
out->x_F=NULL; |
170 |
|
|
out->b_F=NULL; |
171 |
|
|
out->x_C=NULL; |
172 |
|
|
out->b_C=NULL; |
173 |
|
|
out->A=Paso_SparseMatrix_getReference(A_p); |
174 |
|
|
out->solver=NULL; |
175 |
artak |
2828 |
out->Smoother->ID=options->smoother; |
176 |
|
|
out->Smoother->Jacobi=NULL; |
177 |
|
|
out->Smoother->GS=NULL; |
178 |
gross |
2551 |
/*out->GS=Paso_Solver_getGS(A_p,verbose);*/ |
179 |
|
|
out->level=level; |
180 |
artak |
2661 |
out->n=n; |
181 |
artak |
2726 |
out->n_F=n+1; |
182 |
artak |
2661 |
out->n_block=n_block; |
183 |
artak |
2803 |
out->post_sweeps=options->post_sweeps; |
184 |
|
|
out->pre_sweeps=options->pre_sweeps; |
185 |
artak |
2556 |
|
186 |
artak |
2784 |
sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols); |
187 |
|
|
|
188 |
|
|
if (verbose) fprintf(stdout,"Stats: Sparsity of the Coarse Matrix with %d non-zeros (%d,%d) in level %d is %.6f\n",A_p->len,A_p->numRows,A_p->numCols,level,sparsity); |
189 |
|
|
|
190 |
artak |
2802 |
|
191 |
artak |
2803 |
/*if(sparsity>0.01) { |
192 |
artak |
2784 |
level=0; |
193 |
|
|
} |
194 |
artak |
2803 |
*/ |
195 |
artak |
2802 |
|
196 |
gross |
2551 |
if (level==0 || n<=options->min_coarse_matrix_size) { |
197 |
|
|
out->coarsest_level=TRUE; |
198 |
artak |
2784 |
/*out->GS=Paso_Solver_getJacobi(A_p);*/ |
199 |
|
|
|
200 |
artak |
2712 |
#ifdef MKL |
201 |
artak |
2711 |
out->AOffset1=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, out->A->pattern,1,1, FALSE); |
202 |
|
|
#pragma omp parallel for private(i) schedule(static) |
203 |
|
|
for (i=0;i<out->A->len;++i) { |
204 |
|
|
out->AOffset1->val[i]=out->A->val[i]; |
205 |
|
|
} |
206 |
artak |
2712 |
#else |
207 |
|
|
#ifdef UMFPACK |
208 |
|
|
#else |
209 |
artak |
2828 |
if (options->smoother == PASO_JACOBI) |
210 |
|
|
out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p); |
211 |
|
|
else if (options->smoother == PASO_GS) |
212 |
artak |
2831 |
out->Smoother->GS=Paso_Solver_getGS(A_p,verbose); |
213 |
gross |
2551 |
#endif |
214 |
|
|
#endif |
215 |
artak |
2784 |
|
216 |
gross |
2551 |
} else { |
217 |
|
|
out->coarsest_level=FALSE; |
218 |
artak |
2828 |
|
219 |
|
|
if (options->smoother == PASO_JACOBI) |
220 |
|
|
out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p); |
221 |
|
|
else if (options->smoother == PASO_GS) |
222 |
|
|
out->Smoother->GS=Paso_Solver_getGS(A_p,verbose); |
223 |
artak |
2442 |
|
224 |
gross |
2551 |
/* identify independend set of rows/columns */ |
225 |
|
|
#pragma omp parallel for private(i) schedule(static) |
226 |
|
|
for (i=0;i<n;++i) mis_marker[i]=-1; |
227 |
|
|
|
228 |
artak |
2765 |
/*mesuring coarsening time */ |
229 |
|
|
time0=Paso_timer(); |
230 |
|
|
|
231 |
gross |
2551 |
if (options->coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) { |
232 |
artak |
2652 |
Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold); |
233 |
gross |
2551 |
} |
234 |
|
|
else if (options->coarsening_method == PASO_RUGE_STUEBEN_COARSENING) { |
235 |
|
|
Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold); |
236 |
|
|
} |
237 |
|
|
else if (options->coarsening_method == PASO_AGGREGATION_COARSENING) { |
238 |
|
|
Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold); |
239 |
|
|
} |
240 |
artak |
2816 |
else if (options->coarsening_method == PASO_STANDARD_COARSENING) { |
241 |
|
|
Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold); |
242 |
|
|
} |
243 |
gross |
2551 |
else { |
244 |
|
|
/*Default coarseneing*/ |
245 |
artak |
2831 |
/*Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold);*/ |
246 |
|
|
/*Paso_Pattern_Read("RS.spl",n,mis_marker);*/ |
247 |
artak |
2659 |
/*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/ |
248 |
artak |
2828 |
/*Paso_Pattern_greedy_Agg(A_p,mis_marker,options->coarsening_threshold);*/ |
249 |
artak |
2831 |
Paso_Pattern_greedy(A_p->pattern,mis_marker); |
250 |
artak |
2726 |
/*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/ |
251 |
artak |
2760 |
|
252 |
gross |
2551 |
} |
253 |
artak |
2765 |
|
254 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: Profilining for level %d:\n",level); |
255 |
artak |
2767 |
|
256 |
artak |
2765 |
time0=Paso_timer()-time0; |
257 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: Coarsening: %e\n",time0); |
258 |
artak |
2760 |
|
259 |
artak |
2686 |
#pragma omp parallel for private(i) schedule(static) |
260 |
|
|
for (i = 0; i < n; ++i) counter[i]=mis_marker[i]; |
261 |
artak |
2726 |
|
262 |
artak |
2686 |
out->n_F=Paso_Util_cumsum(n,counter); |
263 |
|
|
|
264 |
artak |
2726 |
if (out->n_F==0) { |
265 |
artak |
2686 |
out->coarsest_level=TRUE; |
266 |
artak |
2726 |
level=0; |
267 |
artak |
2686 |
if (verbose) { |
268 |
artak |
2726 |
/*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/ |
269 |
|
|
printf("AMG coarsening does not eliminate any of the unknowns, switching to Jacobi preconditioner.\n"); |
270 |
artak |
2686 |
} |
271 |
artak |
2726 |
} |
272 |
|
|
else if (out->n_F==n) { |
273 |
|
|
out->coarsest_level=TRUE; |
274 |
|
|
level=0; |
275 |
|
|
if (verbose) { |
276 |
|
|
/*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/ |
277 |
|
|
printf("AMG coarsening eliminates all of the unknowns, switching to Jacobi preconditioner.\n"); |
278 |
|
|
|
279 |
|
|
} |
280 |
artak |
2686 |
} else { |
281 |
artak |
2475 |
|
282 |
artak |
2686 |
if (Paso_noError()) { |
283 |
|
|
|
284 |
|
|
/*#pragma omp parallel for private(i) schedule(static) |
285 |
|
|
for (i = 0; i < n; ++i) counter[i]=mis_marker[i]; |
286 |
|
|
out->n_F=Paso_Util_cumsum(n,counter); |
287 |
|
|
*/ |
288 |
artak |
2760 |
|
289 |
artak |
2686 |
out->mask_F=MEMALLOC(n,index_t); |
290 |
|
|
out->rows_in_F=MEMALLOC(out->n_F,index_t); |
291 |
artak |
2760 |
if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->rows_in_F) ) ) { |
292 |
artak |
2686 |
/* creates an index for F from mask */ |
293 |
|
|
#pragma omp parallel for private(i) schedule(static) |
294 |
|
|
for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1; |
295 |
|
|
#pragma omp parallel for private(i) schedule(static) |
296 |
|
|
for (i = 0; i < n; ++i) { |
297 |
|
|
if (mis_marker[i]) { |
298 |
|
|
out->rows_in_F[counter[i]]=i; |
299 |
|
|
out->mask_F[i]=counter[i]; |
300 |
|
|
} else { |
301 |
|
|
out->mask_F[i]=-1; |
302 |
|
|
} |
303 |
|
|
} |
304 |
|
|
|
305 |
jfenwick |
1851 |
} |
306 |
|
|
} |
307 |
artak |
2831 |
|
308 |
|
|
/* if(level==1) { |
309 |
artak |
2802 |
printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F); |
310 |
|
|
for (i = 0; i < n; ++i) { |
311 |
artak |
2831 |
printf("##%d %d\n",i,!mis_marker[i]); |
312 |
artak |
2802 |
} |
313 |
artak |
2816 |
} |
314 |
|
|
*/ |
315 |
|
|
|
316 |
artak |
2686 |
/*check whether coarsening process actually makes sense to continue. |
317 |
|
|
if coarse matrix at least smaller by 30% then continue, otherwise we stop.*/ |
318 |
|
|
if ((out->n_F*100/n)<30) { |
319 |
|
|
level=1; |
320 |
artak |
1881 |
} |
321 |
gross |
2551 |
|
322 |
artak |
2686 |
if ( Paso_noError()) { |
323 |
|
|
out->n_C=n-out->n_F; |
324 |
|
|
out->rows_in_C=MEMALLOC(out->n_C,index_t); |
325 |
|
|
out->mask_C=MEMALLOC(n,index_t); |
326 |
|
|
if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) { |
327 |
|
|
/* creates an index for C from mask */ |
328 |
|
|
#pragma omp parallel for private(i) schedule(static) |
329 |
|
|
for (i = 0; i < n; ++i) counter[i]=! mis_marker[i]; |
330 |
|
|
Paso_Util_cumsum(n,counter); |
331 |
|
|
#pragma omp parallel for private(i) schedule(static) |
332 |
|
|
for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1; |
333 |
|
|
#pragma omp parallel for private(i) schedule(static) |
334 |
|
|
for (i = 0; i < n; ++i) { |
335 |
|
|
if (! mis_marker[i]) { |
336 |
|
|
out->rows_in_C[counter[i]]=i; |
337 |
|
|
out->mask_C[i]=counter[i]; |
338 |
|
|
} else { |
339 |
|
|
out->mask_C[i]=-1; |
340 |
|
|
} |
341 |
|
|
} |
342 |
|
|
} |
343 |
|
|
} |
344 |
|
|
if ( Paso_noError()) { |
345 |
artak |
2760 |
/* get A_FF block: */ |
346 |
|
|
/* |
347 |
|
|
out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F); |
348 |
artak |
2686 |
out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F); |
349 |
|
|
out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C); |
350 |
artak |
2760 |
*/ |
351 |
|
|
|
352 |
|
|
/*Compute W_FC*/ |
353 |
|
|
/*initialy W_FC=A_FC*/ |
354 |
|
|
out->W_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C); |
355 |
|
|
|
356 |
|
|
/*sprintf(filename,"W_FCbefore_%d",level); |
357 |
|
|
Paso_SparseMatrix_saveMM(out->W_FC,filename); |
358 |
|
|
*/ |
359 |
artak |
2802 |
/* for (i = 0; i < n; ++i) { |
360 |
|
|
printf("##mis_marker[%d]=%d\n",i,mis_marker[i]); |
361 |
|
|
} |
362 |
|
|
*/ |
363 |
artak |
2765 |
time0=Paso_timer(); |
364 |
artak |
2760 |
Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker); |
365 |
artak |
2765 |
time0=Paso_timer()-time0; |
366 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0); |
367 |
artak |
2765 |
|
368 |
artak |
2802 |
|
369 |
|
|
/*sprintf(filename,"W_FCafter_%d",level); |
370 |
artak |
2760 |
Paso_SparseMatrix_saveMM(out->W_FC,filename); |
371 |
|
|
*/ |
372 |
artak |
2802 |
|
373 |
artak |
2760 |
/* get Prolongation and Restriction */ |
374 |
artak |
2765 |
time0=Paso_timer(); |
375 |
artak |
2760 |
out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker); |
376 |
artak |
2765 |
time0=Paso_timer()-time0; |
377 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0); |
378 |
|
|
/*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/ |
379 |
artak |
2760 |
|
380 |
artak |
2802 |
|
381 |
artak |
2831 |
/*sprintf(filename,"P_%d",level); |
382 |
artak |
2760 |
Paso_SparseMatrix_saveMM(out->P,filename); |
383 |
|
|
*/ |
384 |
|
|
|
385 |
artak |
2765 |
time0=Paso_timer(); |
386 |
artak |
2760 |
out->R=Paso_SparseMatrix_getRestriction(out->P); |
387 |
artak |
2765 |
time0=Paso_timer()-time0; |
388 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0); |
389 |
|
|
/*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/ |
390 |
artak |
2765 |
|
391 |
artak |
2831 |
/* |
392 |
|
|
sprintf(filename,"R_%d",level); |
393 |
artak |
2760 |
Paso_SparseMatrix_saveMM(out->R,filename); |
394 |
|
|
*/ |
395 |
|
|
|
396 |
gross |
2551 |
} |
397 |
artak |
2686 |
if ( Paso_noError()) { |
398 |
artak |
2765 |
|
399 |
|
|
time0=Paso_timer(); |
400 |
artak |
2784 |
|
401 |
|
|
Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P); |
402 |
artak |
2831 |
|
403 |
artak |
2784 |
A_c=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp); |
404 |
|
|
|
405 |
artak |
2831 |
/*A_c=Paso_SparseMatrix_loadMM_toCSR("A_C1.mtx");*/ |
406 |
|
|
|
407 |
artak |
2784 |
Paso_SparseMatrix_free(Atemp); |
408 |
|
|
|
409 |
|
|
/*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/ |
410 |
artak |
2765 |
time0=Paso_timer()-time0; |
411 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0); |
412 |
artak |
2765 |
|
413 |
artak |
2784 |
|
414 |
artak |
2760 |
/*Paso_Solver_getCoarseMatrix(A_c, A_p,out->R,out->P);*/ |
415 |
artak |
2767 |
|
416 |
artak |
2831 |
/* |
417 |
|
|
sprintf(filename,"A_C_%d",level); |
418 |
artak |
2760 |
Paso_SparseMatrix_saveMM(A_c,filename); |
419 |
artak |
2767 |
*/ |
420 |
|
|
|
421 |
artak |
2760 |
out->AMG_of_Coarse=Paso_Solver_getAMG(A_c,level-1,options); |
422 |
gross |
2551 |
} |
423 |
artak |
2760 |
|
424 |
artak |
2686 |
/* allocate work arrays for AMG application */ |
425 |
|
|
if (Paso_noError()) { |
426 |
artak |
2767 |
/* |
427 |
|
|
out->x_F=MEMALLOC(n_block*out->n_F,double); |
428 |
artak |
2686 |
out->b_F=MEMALLOC(n_block*out->n_F,double); |
429 |
artak |
2767 |
*/ |
430 |
artak |
2686 |
out->x_C=MEMALLOC(n_block*out->n_C,double); |
431 |
|
|
out->b_C=MEMALLOC(n_block*out->n_C,double); |
432 |
|
|
|
433 |
artak |
2767 |
/*if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {*/ |
434 |
|
|
if ( ! ( Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) { |
435 |
|
|
|
436 |
|
|
/* |
437 |
|
|
#pragma omp parallel for private(i) schedule(static) |
438 |
artak |
2686 |
for (i = 0; i < out->n_F; ++i) { |
439 |
|
|
out->x_F[i]=0.; |
440 |
|
|
out->b_F[i]=0.; |
441 |
|
|
} |
442 |
artak |
2767 |
*/ |
443 |
|
|
|
444 |
artak |
2686 |
#pragma omp parallel for private(i) schedule(static) |
445 |
|
|
for (i = 0; i < out->n_C; ++i) { |
446 |
|
|
out->x_C[i]=0.; |
447 |
|
|
out->b_C[i]=0.; |
448 |
|
|
} |
449 |
|
|
} |
450 |
|
|
} |
451 |
artak |
2760 |
Paso_SparseMatrix_free(A_c); |
452 |
artak |
2686 |
} |
453 |
jfenwick |
1851 |
} |
454 |
|
|
} |
455 |
|
|
TMPMEMFREE(mis_marker); |
456 |
|
|
TMPMEMFREE(counter); |
457 |
artak |
2686 |
|
458 |
jfenwick |
1851 |
if (Paso_noError()) { |
459 |
artak |
2524 |
if (verbose && level>0 && !out->coarsest_level) { |
460 |
gross |
2551 |
printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C); |
461 |
jfenwick |
1851 |
} |
462 |
|
|
return out; |
463 |
|
|
} else { |
464 |
|
|
Paso_Solver_AMG_free(out); |
465 |
|
|
return NULL; |
466 |
|
|
} |
467 |
|
|
} |
468 |
|
|
|
469 |
|
|
/**************************************************************/ |
470 |
|
|
|
471 |
|
|
/* apply AMG precondition b-> x |
472 |
|
|
|
473 |
|
|
in fact it solves |
474 |
|
|
|
475 |
artak |
2381 |
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FC ] [ x_F ] = [b_F] |
476 |
jfenwick |
1851 |
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C] |
477 |
|
|
|
478 |
|
|
in the form |
479 |
|
|
|
480 |
|
|
b->[b_F,b_C] |
481 |
|
|
x_F=invA_FF*b_F |
482 |
|
|
b_C=b_C-A_CF*x_F |
483 |
|
|
x_C=AMG(b_C) |
484 |
|
|
b_F=b_F-A_FC*x_C |
485 |
|
|
x_F=invA_FF*b_F |
486 |
|
|
x<-[x_F,x_C] |
487 |
|
|
|
488 |
|
|
should be called within a parallel region |
489 |
|
|
barrier synconization should be performed to make sure that the input vector available |
490 |
|
|
|
491 |
|
|
*/ |
492 |
|
|
|
493 |
|
|
void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) { |
494 |
artak |
2107 |
dim_t i; |
495 |
artak |
2507 |
double time0=0; |
496 |
gross |
2551 |
double *r=NULL, *x0=NULL; |
497 |
artak |
2831 |
bool_t timing=0; |
498 |
artak |
2784 |
|
499 |
artak |
2803 |
dim_t post_sweeps=amg->post_sweeps; |
500 |
|
|
dim_t pre_sweeps=amg->pre_sweeps; |
501 |
|
|
|
502 |
artak |
2699 |
#ifdef UMFPACK |
503 |
|
|
Paso_UMFPACK_Handler * ptr=NULL; |
504 |
artak |
2507 |
#endif |
505 |
artak |
2784 |
|
506 |
gross |
2551 |
r=MEMALLOC(amg->n,double); |
507 |
|
|
x0=MEMALLOC(amg->n,double); |
508 |
artak |
2507 |
|
509 |
artak |
2524 |
if (amg->coarsest_level) { |
510 |
artak |
2507 |
|
511 |
|
|
time0=Paso_timer(); |
512 |
artak |
2699 |
/*If all unknown are eliminated then Jacobi is the best preconditioner*/ |
513 |
artak |
2831 |
|
514 |
artak |
2726 |
if (amg->n_F==0 || amg->n_F==amg->n) { |
515 |
artak |
2828 |
if(amg->Smoother->ID==PASO_JACOBI) |
516 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b); |
517 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
518 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x,b); |
519 |
artak |
2699 |
} |
520 |
|
|
else { |
521 |
artak |
2712 |
#ifdef MKL |
522 |
|
|
Paso_MKL1(amg->AOffset1,x,b,verbose); |
523 |
|
|
#else |
524 |
|
|
#ifdef UMFPACK |
525 |
gross |
2551 |
ptr=(Paso_UMFPACK_Handler *)(amg->solver); |
526 |
artak |
2831 |
Paso_UMFPACK1(&ptr,amg->A,x,b,timing); |
527 |
gross |
2551 |
amg->solver=(void*) ptr; |
528 |
artak |
2712 |
#else |
529 |
artak |
2828 |
if(amg->Smoother->ID==PASO_JACOBI) |
530 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b); |
531 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
532 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x,b); |
533 |
artak |
2699 |
#endif |
534 |
artak |
2507 |
#endif |
535 |
artak |
2699 |
} |
536 |
artak |
2784 |
|
537 |
artak |
2475 |
time0=Paso_timer()-time0; |
538 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0); |
539 |
artak |
2507 |
|
540 |
jfenwick |
1851 |
} else { |
541 |
artak |
1902 |
/* presmoothing */ |
542 |
artak |
2507 |
time0=Paso_timer(); |
543 |
artak |
2828 |
if(amg->Smoother->ID==PASO_JACOBI) |
544 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b); |
545 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
546 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x,b); |
547 |
artak |
2803 |
|
548 |
|
|
/***********/ |
549 |
artak |
2816 |
if (pre_sweeps>1) { |
550 |
|
|
#pragma omp parallel for private(i) schedule(static) |
551 |
|
|
for (i=0;i<amg->n;++i) r[i]=b[i]; |
552 |
|
|
} |
553 |
artak |
2803 |
|
554 |
|
|
while(pre_sweeps>1) { |
555 |
|
|
#pragma omp parallel for private(i) schedule(static) |
556 |
|
|
for (i=0;i<amg->n;++i) r[i]+=b[i]; |
557 |
|
|
|
558 |
artak |
2816 |
/* Compute the residual r=r-Ax*/ |
559 |
artak |
2803 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
560 |
|
|
/* Go round again*/ |
561 |
artak |
2828 |
|
562 |
|
|
if(amg->Smoother->ID==PASO_JACOBI) |
563 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,r); |
564 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
565 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x,r); |
566 |
|
|
|
567 |
artak |
2803 |
pre_sweeps-=1; |
568 |
|
|
} |
569 |
|
|
/***********/ |
570 |
|
|
|
571 |
artak |
2507 |
time0=Paso_timer()-time0; |
572 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0); |
573 |
artak |
2803 |
/* end of presmoothing */ |
574 |
artak |
1939 |
|
575 |
artak |
2507 |
time0=Paso_timer(); |
576 |
artak |
1890 |
#pragma omp parallel for private(i) schedule(static) |
577 |
|
|
for (i=0;i<amg->n;++i) r[i]=b[i]; |
578 |
|
|
|
579 |
|
|
/*r=b-Ax*/ |
580 |
|
|
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
581 |
artak |
2760 |
|
582 |
artak |
2831 |
/* b_c = R*r */ |
583 |
artak |
2760 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C); |
584 |
artak |
1954 |
|
585 |
artak |
2507 |
time0=Paso_timer()-time0; |
586 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0); |
587 |
artak |
2507 |
|
588 |
jfenwick |
1851 |
/* x_C=AMG(b_C) */ |
589 |
artak |
2760 |
Paso_Solver_solveAMG(amg->AMG_of_Coarse,amg->x_C,amg->b_C); |
590 |
artak |
1890 |
|
591 |
artak |
2507 |
time0=Paso_timer(); |
592 |
artak |
2726 |
|
593 |
artak |
2831 |
/* x_0 = P*x_c */ |
594 |
artak |
2760 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0); |
595 |
|
|
|
596 |
|
|
/* x=x+x0 */ |
597 |
artak |
1954 |
#pragma omp parallel for private(i) schedule(static) |
598 |
artak |
2760 |
for (i=0;i<amg->n;++i) x[i]+=x0[i]; |
599 |
artak |
1954 |
|
600 |
artak |
2767 |
/*postsmoothing*/ |
601 |
artak |
2803 |
|
602 |
artak |
2767 |
time0=Paso_timer(); |
603 |
|
|
#pragma omp parallel for private(i) schedule(static) |
604 |
|
|
for (i=0;i<amg->n;++i) r[i]=b[i]; |
605 |
|
|
|
606 |
|
|
/*r=b-Ax */ |
607 |
|
|
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
608 |
artak |
2828 |
if(amg->Smoother->ID==PASO_JACOBI) |
609 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r); |
610 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
611 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x0,r); |
612 |
artak |
2767 |
|
613 |
|
|
#pragma omp parallel for private(i) schedule(static) |
614 |
|
|
for (i=0;i<amg->n;++i) { |
615 |
|
|
x[i]+=x0[i]; |
616 |
|
|
} |
617 |
artak |
2831 |
|
618 |
artak |
2803 |
/***************/ |
619 |
|
|
while(post_sweeps>1) { |
620 |
|
|
|
621 |
|
|
#pragma omp parallel for private(i) schedule(static) |
622 |
|
|
for (i=0;i<amg->n;++i) r[i]=b[i]; |
623 |
|
|
|
624 |
|
|
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
625 |
artak |
2828 |
|
626 |
|
|
if(amg->Smoother->ID==PASO_JACOBI) |
627 |
|
|
Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r); |
628 |
|
|
else if (amg->Smoother->ID==PASO_GS) |
629 |
|
|
Paso_Solver_solveGS(amg->Smoother->GS,x0,r); |
630 |
|
|
|
631 |
artak |
2803 |
#pragma omp parallel for private(i) schedule(static) |
632 |
|
|
for (i=0;i<amg->n;++i) { |
633 |
|
|
x[i]+=x0[i]; |
634 |
|
|
} |
635 |
|
|
post_sweeps-=1; |
636 |
|
|
} |
637 |
|
|
/**************/ |
638 |
artak |
2767 |
|
639 |
|
|
time0=Paso_timer()-time0; |
640 |
artak |
2831 |
if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0); |
641 |
artak |
2767 |
|
642 |
|
|
/*end of postsmoothing*/ |
643 |
artak |
2107 |
|
644 |
artak |
2760 |
} |
645 |
artak |
1890 |
MEMFREE(r); |
646 |
artak |
1954 |
MEMFREE(x0); |
647 |
artak |
2760 |
|
648 |
gross |
2551 |
return; |
649 |
jfenwick |
1851 |
} |