/[escript]/trunk/paso/src/SparseMatrix_AMGcomponents.c
ViewVC logotype

Contents of /trunk/paso/src/SparseMatrix_AMGcomponents.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2759 - (show annotations)
Thu Nov 19 05:19:47 2009 UTC (9 years, 10 months ago) by artak
File MIME type: text/plain
File size: 12167 byte(s)
methods needed for the new AMG preconditioner
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

  ViewVC Help
Powered by ViewVC 1.1.26