/[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 2828 - (show annotations)
Tue Dec 22 01:24:40 2009 UTC (9 years, 7 months ago) by artak
File MIME type: text/plain
File size: 11099 byte(s)
Smoother for AMG now can be selected from python. Now only Jacobi and Gauss-Seidel are available as smoothers.
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 iptr=out->pattern->ptr[i];
88 out->val[iptr]=1;
89 }
90 }
91
92 /* clean up */
93 if (index_list!=NULL) {
94 #pragma omp parallel for private(i) schedule(static)
95 for(i=0;i<n;++i) Paso_IndexList_free(index_list[i].extension);
96 }
97 TMPMEMFREE(index_list);
98 Paso_Pattern_free(outpattern);
99 return out;
100 }
101
102
103 /* Restriction matrix R=P^T */
104
105 Paso_SparseMatrix* Paso_SparseMatrix_getRestriction(Paso_SparseMatrix* P){
106
107 Paso_Pattern *outpattern=NULL;
108 Paso_SparseMatrix *out=NULL;
109
110 Paso_IndexList* index_list=NULL;
111 dim_t C=P->numCols;
112 dim_t F=P->numRows-C;
113 dim_t n=C+F;
114 dim_t i,j,k=0;
115 index_t iptr,jptr;
116
117 index_list=TMPMEMALLOC(C,Paso_IndexList);
118 if (! Paso_checkPtr(index_list)) {
119 #pragma omp parallel for private(i) schedule(static)
120 for(i=0;i<C;++i) {
121 index_list[i].extension=NULL;
122 index_list[i].n=0;
123 }
124 }
125
126
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 a_ii=0;
166 double sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos;
167
168 dim_t n=A->numRows;
169 dim_t F=W_FC->numRows;
170 dim_t i,j,k;
171
172 index_t iPtr,dptr=0;
173 /*index_t *index, *where_p;*/
174
175 alpha=TMPMEMALLOC(F,double);
176 beta=TMPMEMALLOC(F,double);
177
178
179 k=0;
180 for (i = 0; i < n; ++i) {
181 if(mis_marker[i]) {
182 alpha[k]=0;
183 beta[k]=0;
184 sum_all_neg=0;
185 sum_all_pos=0;
186 sum_strong_neg=0;
187 sum_strong_pos=0;
188 /*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*/
189 for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {
190 j=A->pattern->index[iPtr];
191 if(j!=i) {
192 if(A->val[iPtr]<0) {
193 sum_all_neg+=A->val[iPtr];
194 if(!mis_marker[j]) {
195 sum_strong_neg+=A->val[iPtr];
196 }
197 }
198 else if(A->val[iPtr]>0) {
199 sum_all_pos+=A->val[iPtr];
200 if(!mis_marker[j]) {
201 sum_strong_pos+=A->val[iPtr];
202 }
203 }
204 }
205 else {
206 a_ii=A->val[iPtr];
207 dptr=iPtr;
208 }
209 }
210
211 if(sum_strong_neg!=0) {
212 alpha[k]=sum_all_neg/(sum_strong_neg);
213 }
214 if(sum_strong_pos!=0) {
215 beta[k]=sum_all_pos/(sum_strong_pos);
216 }
217 else {
218 a_ii+=sum_all_pos;
219 beta[k]=0;
220 }
221 alpha[k]=alpha[k]/(a_ii);
222 beta[k]=beta[k]/(a_ii);
223 k++;
224 /*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]);*/
225 }
226 }
227 #pragma omp parallel for private(i,iPtr) schedule(static)
228 for (i = 0; i < W_FC->numRows; ++i) {
229 for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
230 if(W_FC->val[iPtr]<0) {
231 W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr];
232 }
233 else if(W_FC->val[iPtr]>0) {
234 W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr];
235 }
236 }
237 }
238
239 TMPMEMFREE(alpha);
240 TMPMEMFREE(beta);
241 }
242
243
244 /*A_c=R*A*P taking into account coarseneing */
245 /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {
246
247 Paso_SparseMatrix *temp=NULL;
248 Paso_SparseMatrix *out=NULL;
249
250 temp=Paso_SparseMatrix_MatrixMatrix(A,P);
251 out=Paso_SparseMatrix_MatrixMatrix(R,temp);
252
253 Paso_SparseMatrix_free(temp);
254
255 return out;
256 }
257 */
258
259 Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) {
260 index_t iptrA,iptrB,iptrC;
261 dim_t i,j=0,k,l;
262 Paso_Pattern* outpattern=NULL;
263 Paso_SparseMatrix *out=NULL;
264 double sum,b_lj=0;
265
266 double time0=0;
267 bool_t verbose=0;
268
269 time0=Paso_timer();
270
271 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
272 out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
273
274 time0=Paso_timer()-time0;
275 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0);
276
277 time0=Paso_timer();
278
279 #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static)
280 for(i = 0; i < out->numRows; i++) {
281 for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
282 j = out->pattern->index[iptrC];
283 sum=0;
284 for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
285 k=A->pattern->index[iptrA];
286 /*can be replaced by bsearch */
287 b_lj=0;
288 for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
289 l=B->pattern->index[iptrB];
290 if(l==j) {
291 b_lj=B->val[iptrB];
292 break;
293 }
294 }
295 sum+=A->val[iptrA]*b_lj;
296 }
297 out->val[iptrC]=sum;
298 }
299 }
300
301 time0=Paso_timer()-time0;
302 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0);
303
304 Paso_Pattern_free(outpattern);
305
306 return out;
307 }
308
309
310 Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {
311
312 index_t iptrA_c,iptrR,iptrP,iptrA;
313 dim_t i,j,k,l,m;
314 double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
315 Paso_Pattern* temp=NULL;
316 Paso_Pattern* outpattern=NULL;
317 Paso_SparseMatrix *A_c=NULL;
318
319 double time0=0;
320 bool_t verbose=0;
321
322 time0=Paso_timer();
323
324 temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
325 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
326 A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
327
328 time0=Paso_timer()-time0;
329 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0);
330
331 /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
332
333 time0=Paso_timer();
334
335 #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)
336 for(i = 0; i < A_c->numRows; i++) {
337 for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
338 j = A_c->pattern->index[iptrA_c];
339 second_sum=0;
340 for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
341 k = R->pattern->index[iptrR];
342 first_sum=0;
343 for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
344 l= A->pattern->index[iptrA];
345 p_lj=0;
346 /*This loop can be replaced by bsearch */
347 for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
348 m=P->pattern->index[iptrP];
349 if(m==j) {
350 p_lj=P->val[iptrP];
351 break;
352 }
353 }
354 a_kl=A->val[iptrA];
355 first_sum+=a_kl*p_lj;
356 }
357 r_ik=R->val[iptrR];
358 second_sum+=r_ik*first_sum;
359 }
360 A_c->val[iptrA_c]=second_sum;
361 }
362 }
363 time0=Paso_timer()-time0;
364 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0);
365
366 Paso_Pattern_free(outpattern);
367 Paso_Pattern_free(temp);
368
369 return A_c;
370 }
371

  ViewVC Help
Powered by ViewVC 1.1.26