/[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 2807 - (show annotations)
Mon Dec 7 00:02:55 2009 UTC (9 years, 10 months ago) by artak
File MIME type: text/plain
File size: 11237 byte(s)
Post and pre smoothing parameters are added for AMLI. Little bit reengineering according to Axelson's original AMLI paper.
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
59 for (i=0;i<n;++i) {
60 if (mis_marker[i]) {
61 for (iptr=W->pattern->ptr[k];iptr<W->pattern->ptr[k+1]; ++iptr) {
62 j=W->pattern->index[iptr];
63 Paso_IndexList_insertIndex(&(index_list[i]),j);
64 }
65 k++;
66 }
67 else {
68 Paso_IndexList_insertIndex(&(index_list[i]),i-k);
69 }
70 }
71
72 outpattern=Paso_IndexList_createPattern(0, n,index_list,0,W->numCols,0);
73 out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE);
74
75 k=0;
76
77 for (i=0;i<n;++i) {
78 if (mis_marker[i]) {
79 wptr=W->pattern->ptr[k];
80 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
81 out->val[iptr]=W->val[wptr];
82 wptr++;
83 }
84 k++;
85 }
86 else {
87 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr){
88 out->val[iptr]=1;
89 }
90 }
91 }
92
93 /* clean up */
94 if (index_list!=NULL) {
95 #pragma omp parallel for private(i) schedule(static)
96 for(i=0;i<n;++i) Paso_IndexList_free(index_list[i].extension);
97 }
98 TMPMEMFREE(index_list);
99 Paso_Pattern_free(outpattern);
100 return out;
101 }
102
103
104 /* Restriction matrix R=P^T */
105
106 Paso_SparseMatrix* Paso_SparseMatrix_getRestriction(Paso_SparseMatrix* P){
107
108 Paso_Pattern *outpattern=NULL;
109 Paso_SparseMatrix *out=NULL;
110
111 Paso_IndexList* index_list=NULL;
112 dim_t C=P->numCols;
113 dim_t F=P->numRows-C;
114 dim_t n=C+F;
115 dim_t i,j,k=0;
116 index_t iptr,jptr;
117
118 index_list=TMPMEMALLOC(C,Paso_IndexList);
119 if (! Paso_checkPtr(index_list)) {
120 #pragma omp parallel for private(i) schedule(static)
121 for(i=0;i<C;++i) {
122 index_list[i].extension=NULL;
123 index_list[i].n=0;
124 }
125 }
126
127
128 for (i=0;i<n;++i) {
129 for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {
130 j=P->pattern->index[iptr];
131 Paso_IndexList_insertIndex(&(index_list[j]),i);
132 }
133 }
134
135 outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0);
136 out=Paso_SparseMatrix_alloc(P->type,outpattern,1,1,TRUE);
137
138 for (i=0;i<out->numRows;++i) {
139 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
140 j=out->pattern->index[iptr];
141 /*This for later will be replaced by bsearch!!*/
142 for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
143 k=P->pattern->index[jptr];
144 if(k==i) {
145 out->val[iptr]=P->val[jptr];
146 }
147 }
148 }
149 }
150
151 /* clean up */
152 if (index_list!=NULL) {
153 #pragma omp parallel for private(i) schedule(static)
154 for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
155 }
156 TMPMEMFREE(index_list);
157 Paso_Pattern_free(outpattern);
158 return out;
159 }
160
161
162 void Paso_SparseMatrix_updateWeights(Paso_SparseMatrix* A,Paso_SparseMatrix* W_FC, index_t* mis_marker){
163
164 double *alpha;
165 double *beta;
166 double a_ii=0;
167 double sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos;
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,dptr=0;
174 /*index_t *index, *where_p;*/
175
176 alpha=TMPMEMALLOC(F,double);
177 beta=TMPMEMALLOC(F,double);
178
179
180 k=0;
181 for (i = 0; i < n; ++i) {
182 if(mis_marker[i]) {
183 alpha[k]=0;
184 beta[k]=0;
185 sum_all_neg=0;
186 sum_all_pos=0;
187 sum_strong_neg=0;
188 sum_strong_pos=0;
189 /*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*/
190 for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {
191 j=A->pattern->index[iPtr];
192 if(j!=i) {
193 if(A->val[iPtr]<0) {
194 sum_all_neg+=A->val[iPtr];
195 if(!mis_marker[j]) {
196 sum_strong_neg+=A->val[iPtr];
197 }
198 }
199 else if(A->val[iPtr]>0) {
200 sum_all_pos+=A->val[iPtr];
201 if(!mis_marker[j]) {
202 sum_strong_pos+=A->val[iPtr];
203 }
204 }
205 }
206 else {
207 a_ii=A->val[iPtr];
208 dptr=iPtr;
209 }
210 }
211
212 if(sum_strong_neg!=0) {
213 alpha[k]=sum_all_neg/(sum_strong_neg);
214 }
215 if(sum_strong_pos!=0) {
216 beta[k]=sum_all_pos/(sum_strong_pos);
217 }
218 /* else {
219 a_ii+=beta[k];
220 A->val[dptr]=a_ii;
221 beta[k]=0;
222 }
223 */ alpha[k]=alpha[k]/(a_ii);
224 beta[k]=beta[k]/(a_ii);
225 k++;
226 /*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]);*/
227 }
228 }
229 #pragma omp parallel for private(i,iPtr,j) schedule(static)
230 for (i = 0; i < W_FC->numRows; ++i) {
231 for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
232 j=W_FC->pattern->index[iPtr];
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 TMPMEMFREE(alpha);
243 TMPMEMFREE(beta);
244 }
245
246
247 /*A_c=R*A*P taking into account coarseneing */
248 /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {
249
250 Paso_SparseMatrix *temp=NULL;
251 Paso_SparseMatrix *out=NULL;
252
253 temp=Paso_SparseMatrix_MatrixMatrix(A,P);
254 out=Paso_SparseMatrix_MatrixMatrix(R,temp);
255
256 Paso_SparseMatrix_free(temp);
257
258 return out;
259 }
260 */
261
262 Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) {
263 index_t iptrA,iptrB,iptrC;
264 dim_t i,j=0,k,l;
265 Paso_Pattern* outpattern=NULL;
266 Paso_SparseMatrix *out=NULL;
267 double sum,b_lj=0;
268
269 double time0=0;
270 bool_t verbose=0;
271
272 time0=Paso_timer();
273
274 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
275 out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
276
277 time0=Paso_timer()-time0;
278 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0);
279
280 time0=Paso_timer();
281
282 #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static)
283 for(i = 0; i < out->numRows; i++) {
284 for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
285 j = out->pattern->index[iptrC];
286 sum=0;
287 for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
288 k=A->pattern->index[iptrA];
289 /*can be replaced by bsearch */
290 b_lj=0;
291 for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
292 l=B->pattern->index[iptrB];
293 if(l==j) {
294 b_lj=B->val[iptrB];
295 break;
296 }
297 }
298 sum+=A->val[iptrA]*b_lj;
299 }
300 out->val[iptrC]=sum;
301 }
302 }
303
304 time0=Paso_timer()-time0;
305 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0);
306
307 Paso_Pattern_free(outpattern);
308
309 return out;
310 }
311
312
313 Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {
314
315 index_t iptrA_c,iptrR,iptrP,iptrA;
316 dim_t i,j,k,l,m;
317 double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
318 Paso_Pattern* temp=NULL;
319 Paso_Pattern* outpattern=NULL;
320 Paso_SparseMatrix *A_c=NULL;
321
322 double time0=0;
323 bool_t verbose=0;
324
325 time0=Paso_timer();
326
327 temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
328 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
329 A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
330
331 time0=Paso_timer()-time0;
332 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0);
333
334 /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
335
336 time0=Paso_timer();
337
338 #pragma omp parallel for private(i,iptrA_c,j,second_sum,iptrR,k,first_sum,p_lj,iptrP,m,a_kl,r_ik) schedule(static)
339 for(i = 0; i < A_c->numRows; i++) {
340 for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
341 j = A_c->pattern->index[iptrA_c];
342 second_sum=0;
343 for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
344 k = R->pattern->index[iptrR];
345 first_sum=0;
346 for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
347 l= A->pattern->index[iptrA];
348 p_lj=0;
349 /*This loop can be replaced by bsearch */
350 for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
351 m=P->pattern->index[iptrP];
352 if(m==j) {
353 p_lj=P->val[iptrP];
354 break;
355 }
356 }
357 a_kl=A->val[iptrA];
358 first_sum+=a_kl*p_lj;
359 }
360 r_ik=R->val[iptrR];
361 second_sum+=r_ik*first_sum;
362 }
363 A_c->val[iptrA_c]=second_sum;
364 }
365 }
366 time0=Paso_timer()-time0;
367 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0);
368
369 Paso_Pattern_free(outpattern);
370 Paso_Pattern_free(temp);
371
372 return A_c;
373 }
374

  ViewVC Help
Powered by ViewVC 1.1.26