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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2828 - (hide annotations)
Tue Dec 22 01:24:40 2009 UTC (9 years, 9 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 artak 2759
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 artak 2807
58    
59 artak 2759 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 artak 2802 j=W->pattern->index[iptr];
63     Paso_IndexList_insertIndex(&(index_list[i]),j);
64     }
65     k++;
66 artak 2759 }
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 artak 2784
77 artak 2802 for (i=0;i<n;++i) {
78 artak 2759 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 artak 2828 iptr=out->pattern->ptr[i];
88     out->val[iptr]=1;
89 artak 2759 }
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 artak 2804 double sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos;
167 artak 2759
168     dim_t n=A->numRows;
169     dim_t F=W_FC->numRows;
170     dim_t i,j,k;
171    
172 artak 2804 index_t iPtr,dptr=0;
173     /*index_t *index, *where_p;*/
174 artak 2759
175     alpha=TMPMEMALLOC(F,double);
176     beta=TMPMEMALLOC(F,double);
177    
178 artak 2804
179     k=0;
180 artak 2759 for (i = 0; i < n; ++i) {
181     if(mis_marker[i]) {
182     alpha[k]=0;
183     beta[k]=0;
184 artak 2804 sum_all_neg=0;
185     sum_all_pos=0;
186     sum_strong_neg=0;
187     sum_strong_pos=0;
188 artak 2759 /*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 artak 2804 sum_all_neg+=A->val[iPtr];
194 artak 2759 if(!mis_marker[j]) {
195 artak 2804 sum_strong_neg+=A->val[iPtr];
196 artak 2759 }
197     }
198     else if(A->val[iPtr]>0) {
199 artak 2804 sum_all_pos+=A->val[iPtr];
200 artak 2759 if(!mis_marker[j]) {
201 artak 2804 sum_strong_pos+=A->val[iPtr];
202 artak 2759 }
203     }
204     }
205 artak 2804 else {
206     a_ii=A->val[iPtr];
207     dptr=iPtr;
208     }
209 artak 2759 }
210 artak 2804
211     if(sum_strong_neg!=0) {
212     alpha[k]=sum_all_neg/(sum_strong_neg);
213 artak 2759 }
214 artak 2804 if(sum_strong_pos!=0) {
215     beta[k]=sum_all_pos/(sum_strong_pos);
216 artak 2759 }
217 artak 2828 else {
218     a_ii+=sum_all_pos;
219 artak 2804 beta[k]=0;
220     }
221 artak 2828 alpha[k]=alpha[k]/(a_ii);
222 artak 2804 beta[k]=beta[k]/(a_ii);
223     k++;
224 artak 2759 /*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 artak 2828 #pragma omp parallel for private(i,iPtr) schedule(static)
228 artak 2759 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 artak 2784 double time0=0;
267     bool_t verbose=0;
268    
269     time0=Paso_timer();
270    
271 artak 2759 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
272     out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
273    
274 artak 2784 time0=Paso_timer()-time0;
275     if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0);
276 artak 2759
277 artak 2784 time0=Paso_timer();
278    
279     #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static)
280 artak 2759 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 artak 2784 b_lj=0;
288 artak 2759 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 artak 2784
301     time0=Paso_timer()-time0;
302     if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0);
303 artak 2759
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 artak 2784 double time0=0;
320     bool_t verbose=0;
321    
322     time0=Paso_timer();
323    
324 artak 2759 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 artak 2784 time0=Paso_timer()-time0;
329     if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0);
330 artak 2759
331     /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
332    
333 artak 2784 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 artak 2759 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 artak 2767 /*This loop can be replaced by bsearch */
347 artak 2759 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 artak 2784 }
363     time0=Paso_timer()-time0;
364     if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0);
365 artak 2759
366 artak 2784 Paso_Pattern_free(outpattern);
367     Paso_Pattern_free(temp);
368 artak 2759
369     return A_c;
370     }
371    

  ViewVC Help
Powered by ViewVC 1.1.26