1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2009 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: SparseMatrix_AMGcomponents */ |
18 |
|
19 |
/**************************************************************/ |
20 |
|
21 |
/* Author: Artak Amirbekyan, artak@uq.edu.au */ |
22 |
|
23 |
/**************************************************************/ |
24 |
|
25 |
#include "Paso.h" |
26 |
#include "SparseMatrix.h" |
27 |
#include "PasoUtil.h" |
28 |
|
29 |
/************************************************************** |
30 |
|
31 |
Methods nessecary for AMG preconditioner |
32 |
|
33 |
*/ |
34 |
|
35 |
/* Prolongation matrix is constracted by W_FC as P<- [W_FC] |
36 |
[I_CC] |
37 |
*/ |
38 |
|
39 |
Paso_SparseMatrix* Paso_SparseMatrix_getProlongation(Paso_SparseMatrix* W, index_t* mis_marker){ |
40 |
|
41 |
Paso_Pattern *outpattern=NULL; |
42 |
Paso_SparseMatrix *out=NULL; |
43 |
index_t iptr,wptr; |
44 |
|
45 |
Paso_IndexList* index_list=NULL; |
46 |
dim_t n=W->numRows+W->numCols; |
47 |
dim_t i,j,k=0; |
48 |
|
49 |
index_list=TMPMEMALLOC(n,Paso_IndexList); |
50 |
if (! Paso_checkPtr(index_list)) { |
51 |
#pragma omp parallel for private(i) schedule(static) |
52 |
for(i=0;i<n;++i) { |
53 |
index_list[i].extension=NULL; |
54 |
index_list[i].n=0; |
55 |
} |
56 |
} |
57 |
|
58 |
for (i=0;i<n;++i) { |
59 |
if (mis_marker[i]) { |
60 |
for (iptr=W->pattern->ptr[k];iptr<W->pattern->ptr[k+1]; ++iptr) { |
61 |
j=W->pattern->index[iptr]; |
62 |
Paso_IndexList_insertIndex(&(index_list[i]),j); |
63 |
} |
64 |
k++; |
65 |
} |
66 |
else { |
67 |
Paso_IndexList_insertIndex(&(index_list[i]),i-k); |
68 |
} |
69 |
} |
70 |
|
71 |
outpattern=Paso_IndexList_createPattern(0, n,index_list,0,W->numCols,0); |
72 |
out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE); |
73 |
|
74 |
k=0; |
75 |
for (i=0;i<out->numRows;++i) { |
76 |
if (mis_marker[i]) { |
77 |
wptr=W->pattern->ptr[k]; |
78 |
for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) { |
79 |
out->val[iptr]=W->val[wptr]; |
80 |
wptr++; |
81 |
} |
82 |
k++; |
83 |
} |
84 |
else { |
85 |
for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr){ |
86 |
out->val[iptr]=1; |
87 |
} |
88 |
} |
89 |
} |
90 |
|
91 |
/* clean up */ |
92 |
if (index_list!=NULL) { |
93 |
#pragma omp parallel for private(i) schedule(static) |
94 |
for(i=0;i<n;++i) Paso_IndexList_free(index_list[i].extension); |
95 |
} |
96 |
TMPMEMFREE(index_list); |
97 |
Paso_Pattern_free(outpattern); |
98 |
return out; |
99 |
} |
100 |
|
101 |
|
102 |
/* Restriction matrix R=P^T */ |
103 |
|
104 |
Paso_SparseMatrix* Paso_SparseMatrix_getRestriction(Paso_SparseMatrix* P){ |
105 |
|
106 |
Paso_Pattern *outpattern=NULL; |
107 |
Paso_SparseMatrix *out=NULL; |
108 |
|
109 |
Paso_IndexList* index_list=NULL; |
110 |
dim_t C=P->numCols; |
111 |
dim_t F=P->numRows-C; |
112 |
dim_t n=C+F; |
113 |
dim_t i,j,k=0; |
114 |
index_t iptr,jptr; |
115 |
|
116 |
index_list=TMPMEMALLOC(C,Paso_IndexList); |
117 |
if (! Paso_checkPtr(index_list)) { |
118 |
#pragma omp parallel for private(i) schedule(static) |
119 |
for(i=0;i<C;++i) { |
120 |
index_list[i].extension=NULL; |
121 |
index_list[i].n=0; |
122 |
} |
123 |
} |
124 |
|
125 |
|
126 |
/*#pragma omp parallel for private(i,iptr,j) schedule(static)*/ |
127 |
for (i=0;i<n;++i) { |
128 |
for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) { |
129 |
j=P->pattern->index[iptr]; |
130 |
Paso_IndexList_insertIndex(&(index_list[j]),i); |
131 |
} |
132 |
} |
133 |
|
134 |
outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0); |
135 |
out=Paso_SparseMatrix_alloc(P->type,outpattern,1,1,TRUE); |
136 |
|
137 |
for (i=0;i<out->numRows;++i) { |
138 |
for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) { |
139 |
j=out->pattern->index[iptr]; |
140 |
/*This for later will be replaced by bsearch!!*/ |
141 |
for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) { |
142 |
k=P->pattern->index[jptr]; |
143 |
if(k==i) { |
144 |
out->val[iptr]=P->val[jptr]; |
145 |
} |
146 |
} |
147 |
} |
148 |
} |
149 |
|
150 |
/* clean up */ |
151 |
if (index_list!=NULL) { |
152 |
#pragma omp parallel for private(i) schedule(static) |
153 |
for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension); |
154 |
} |
155 |
TMPMEMFREE(index_list); |
156 |
Paso_Pattern_free(outpattern); |
157 |
return out; |
158 |
} |
159 |
|
160 |
|
161 |
void Paso_SparseMatrix_updateWeights(Paso_SparseMatrix* A,Paso_SparseMatrix* W_FC, index_t* mis_marker){ |
162 |
|
163 |
double *alpha; |
164 |
double *beta; |
165 |
double *alpha_den; |
166 |
double *beta_den; |
167 |
double a_ii=0; |
168 |
|
169 |
dim_t n=A->numRows; |
170 |
dim_t F=W_FC->numRows; |
171 |
dim_t i,j,k; |
172 |
|
173 |
index_t iPtr,*index, *where_p; |
174 |
|
175 |
alpha=TMPMEMALLOC(F,double); |
176 |
beta=TMPMEMALLOC(F,double); |
177 |
alpha_den=TMPMEMALLOC(F,double); |
178 |
beta_den=TMPMEMALLOC(F,double); |
179 |
|
180 |
k=0; |
181 |
for (i = 0; i < n; ++i) { |
182 |
if(mis_marker[i]) { |
183 |
iPtr=A->pattern->ptr[i]; |
184 |
index=&(A->pattern->index[iPtr]); |
185 |
where_p=(index_t*)bsearch(&i, |
186 |
index, |
187 |
A->pattern->ptr[i + 1]-A->pattern->ptr[i], |
188 |
sizeof(index_t), |
189 |
Paso_comparIndex); |
190 |
if (where_p!=NULL) { |
191 |
a_ii=A->val[iPtr+(index_t)(where_p-index)]; |
192 |
} |
193 |
else { |
194 |
Paso_setError(VALUE_ERROR,"Paso_Solver_getAMG: Main diagonal element of system matrix is missing."); |
195 |
} |
196 |
alpha[k]=0; |
197 |
beta[k]=0; |
198 |
alpha_den[k]=0; |
199 |
beta_den[k]=0; |
200 |
/*compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i*/ |
201 |
for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) { |
202 |
j=A->pattern->index[iPtr]; |
203 |
if(j!=i) { |
204 |
if(A->val[iPtr]<0) { |
205 |
alpha[k]+=A->val[iPtr]; |
206 |
if(!mis_marker[j]) { |
207 |
alpha_den[k]+=A->val[iPtr]; |
208 |
} |
209 |
} |
210 |
else if(A->val[iPtr]>0) { |
211 |
beta[k]+=A->val[iPtr]; |
212 |
if(!mis_marker[j]) { |
213 |
beta_den[k]+=A->val[iPtr]; |
214 |
} |
215 |
} |
216 |
} |
217 |
} |
218 |
if(alpha_den[k]!=0) { |
219 |
alpha[k]=alpha[k]/(alpha_den[k]*a_ii); |
220 |
} |
221 |
if(beta_den[k]!=0) { |
222 |
beta[k]=beta[k]/(beta_den[k]*a_ii); |
223 |
} |
224 |
k++; |
225 |
/*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/ |
226 |
} |
227 |
} |
228 |
|
229 |
for (i = 0; i < W_FC->numRows; ++i) { |
230 |
for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) { |
231 |
j=W_FC->pattern->index[iPtr]; |
232 |
/*if(!mis_marker[j]) {*/ |
233 |
if(W_FC->val[iPtr]<0) { |
234 |
W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr]; |
235 |
} |
236 |
else if(W_FC->val[iPtr]>0) { |
237 |
W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr]; |
238 |
} |
239 |
/*}*/ |
240 |
} |
241 |
} |
242 |
|
243 |
TMPMEMFREE(alpha); |
244 |
TMPMEMFREE(beta); |
245 |
TMPMEMFREE(alpha_den); |
246 |
TMPMEMFREE(beta_den); |
247 |
} |
248 |
|
249 |
/*A_c=R*A*P taking into account coarseneing */ |
250 |
/* |
251 |
void Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A_c, Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) { |
252 |
|
253 |
index_t iptrA_c,iptrR,iptrP,iptrA; |
254 |
dim_t i,j,k,l,m; |
255 |
double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0; |
256 |
*/ |
257 |
/*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/ |
258 |
/* |
259 |
for(i = 0; i < A_c->numRows; i++) { |
260 |
for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) { |
261 |
j = A_c->pattern->index[iptrA_c]; |
262 |
second_sum=0; |
263 |
for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) { |
264 |
k = R->pattern->index[iptrR]; |
265 |
first_sum=0; |
266 |
for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) { |
267 |
l= A->pattern->index[iptrA]; |
268 |
*/ /*This loop has to be replaced by bsearch */ |
269 |
/* for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) { |
270 |
m=P->pattern->index[iptrP]; |
271 |
if(m==j) { |
272 |
p_lj=P->val[iptrP]; |
273 |
break; |
274 |
} |
275 |
} |
276 |
a_kl=A->val[iptrA]; |
277 |
first_sum+=a_kl*p_lj; |
278 |
} |
279 |
r_ik=R->val[iptrR]; |
280 |
second_sum+=r_ik*first_sum; |
281 |
if(i==0) { |
282 |
printf("row %d, col %d, k=%d, firstsum = %e, r_ik=%e \n",i,j,k,first_sum,r_ik); |
283 |
} |
284 |
} |
285 |
if(i==0) { |
286 |
printf("row %d, col %d, secondsum = %e \n",i,j,second_sum); |
287 |
} |
288 |
A_c->val[iptrA_c]=second_sum; |
289 |
} |
290 |
} |
291 |
|
292 |
} |
293 |
*/ |
294 |
|
295 |
/*A_c=R*A*P taking into account coarseneing */ |
296 |
/*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) { |
297 |
|
298 |
Paso_SparseMatrix *temp=NULL; |
299 |
Paso_SparseMatrix *out=NULL; |
300 |
|
301 |
temp=Paso_SparseMatrix_MatrixMatrix(A,P); |
302 |
out=Paso_SparseMatrix_MatrixMatrix(R,temp); |
303 |
|
304 |
Paso_SparseMatrix_free(temp); |
305 |
|
306 |
return out; |
307 |
} |
308 |
*/ |
309 |
|
310 |
Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) { |
311 |
index_t iptrA,iptrB,iptrC; |
312 |
dim_t i,j=0,k,l; |
313 |
Paso_Pattern* outpattern=NULL; |
314 |
Paso_SparseMatrix *out=NULL; |
315 |
double sum,b_lj=0; |
316 |
|
317 |
outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern); |
318 |
out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE); |
319 |
|
320 |
|
321 |
for(i = 0; i < out->numRows; i++) { |
322 |
for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) { |
323 |
j = out->pattern->index[iptrC]; |
324 |
sum=0; |
325 |
b_lj=0; |
326 |
for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) { |
327 |
k=A->pattern->index[iptrA]; |
328 |
/*can be replaced by bsearch */ |
329 |
for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) { |
330 |
l=B->pattern->index[iptrB]; |
331 |
if(l==j) { |
332 |
b_lj=B->val[iptrB]; |
333 |
break; |
334 |
} |
335 |
} |
336 |
sum+=A->val[iptrA]*b_lj; |
337 |
} |
338 |
out->val[iptrC]=sum; |
339 |
} |
340 |
} |
341 |
|
342 |
Paso_Pattern_free(outpattern); |
343 |
|
344 |
return out; |
345 |
} |
346 |
|
347 |
|
348 |
Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) { |
349 |
|
350 |
index_t iptrA_c,iptrR,iptrP,iptrA; |
351 |
dim_t i,j,k,l,m; |
352 |
double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0; |
353 |
Paso_Pattern* temp=NULL; |
354 |
Paso_Pattern* outpattern=NULL; |
355 |
Paso_SparseMatrix *A_c=NULL; |
356 |
|
357 |
temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern); |
358 |
outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp); |
359 |
A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE); |
360 |
|
361 |
|
362 |
/*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/ |
363 |
|
364 |
for(i = 0; i < A_c->numRows; i++) { |
365 |
for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) { |
366 |
j = A_c->pattern->index[iptrA_c]; |
367 |
second_sum=0; |
368 |
for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) { |
369 |
k = R->pattern->index[iptrR]; |
370 |
first_sum=0; |
371 |
for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) { |
372 |
l= A->pattern->index[iptrA]; |
373 |
p_lj=0; |
374 |
/*This loop has to be replaced by bsearch */ |
375 |
for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) { |
376 |
m=P->pattern->index[iptrP]; |
377 |
if(m==j) { |
378 |
p_lj=P->val[iptrP]; |
379 |
break; |
380 |
} |
381 |
} |
382 |
a_kl=A->val[iptrA]; |
383 |
first_sum+=a_kl*p_lj; |
384 |
} |
385 |
r_ik=R->val[iptrR]; |
386 |
second_sum+=r_ik*first_sum; |
387 |
} |
388 |
A_c->val[iptrA_c]=second_sum; |
389 |
} |
390 |
} |
391 |
|
392 |
Paso_Pattern_free(outpattern); |
393 |
Paso_Pattern_free(temp); |
394 |
|
395 |
return A_c; |
396 |
} |
397 |
|