/[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 2802 - (hide annotations)
Thu Dec 3 01:51:55 2009 UTC (9 years, 11 months ago) by artak
File MIME type: text/plain
File size: 11610 byte(s)
Aggressive coarsening is implemented. Removes more elements but adds more iteration steps.
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    
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 artak 2802 j=W->pattern->index[iptr];
62     Paso_IndexList_insertIndex(&(index_list[i]),j);
63     }
64     k++;
65 artak 2759 }
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 artak 2784
76 artak 2802 for (i=0;i<n;++i) {
77 artak 2759 if (mis_marker[i]) {
78     wptr=W->pattern->ptr[k];
79     for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
80     out->val[iptr]=W->val[wptr];
81     wptr++;
82     }
83     k++;
84     }
85     else {
86     for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr){
87     out->val[iptr]=1;
88     }
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 *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 artak 2784 #pragma omp parallel for private(i,iPtr,j) schedule(static)
229 artak 2759 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    
250     /*A_c=R*A*P taking into account coarseneing */
251     /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {
252    
253     Paso_SparseMatrix *temp=NULL;
254     Paso_SparseMatrix *out=NULL;
255    
256     temp=Paso_SparseMatrix_MatrixMatrix(A,P);
257     out=Paso_SparseMatrix_MatrixMatrix(R,temp);
258    
259     Paso_SparseMatrix_free(temp);
260    
261     return out;
262     }
263     */
264    
265     Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) {
266     index_t iptrA,iptrB,iptrC;
267     dim_t i,j=0,k,l;
268     Paso_Pattern* outpattern=NULL;
269     Paso_SparseMatrix *out=NULL;
270     double sum,b_lj=0;
271    
272 artak 2784 double time0=0;
273     bool_t verbose=0;
274    
275     time0=Paso_timer();
276    
277 artak 2759 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
278     out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
279    
280 artak 2784 time0=Paso_timer()-time0;
281     if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0);
282 artak 2759
283 artak 2784 time0=Paso_timer();
284    
285     #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static)
286 artak 2759 for(i = 0; i < out->numRows; i++) {
287     for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
288     j = out->pattern->index[iptrC];
289     sum=0;
290     for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
291     k=A->pattern->index[iptrA];
292     /*can be replaced by bsearch */
293 artak 2784 b_lj=0;
294 artak 2759 for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
295     l=B->pattern->index[iptrB];
296     if(l==j) {
297     b_lj=B->val[iptrB];
298     break;
299     }
300     }
301     sum+=A->val[iptrA]*b_lj;
302     }
303     out->val[iptrC]=sum;
304     }
305     }
306 artak 2784
307     time0=Paso_timer()-time0;
308     if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0);
309 artak 2759
310     Paso_Pattern_free(outpattern);
311    
312     return out;
313     }
314    
315    
316     Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {
317    
318     index_t iptrA_c,iptrR,iptrP,iptrA;
319     dim_t i,j,k,l,m;
320     double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
321     Paso_Pattern* temp=NULL;
322     Paso_Pattern* outpattern=NULL;
323     Paso_SparseMatrix *A_c=NULL;
324    
325 artak 2784 double time0=0;
326     bool_t verbose=0;
327    
328     time0=Paso_timer();
329    
330 artak 2759 temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
331     outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
332     A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
333    
334 artak 2784 time0=Paso_timer()-time0;
335     if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0);
336 artak 2759
337     /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
338    
339 artak 2784 time0=Paso_timer();
340    
341     #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)
342 artak 2759 for(i = 0; i < A_c->numRows; i++) {
343     for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
344     j = A_c->pattern->index[iptrA_c];
345     second_sum=0;
346     for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
347     k = R->pattern->index[iptrR];
348     first_sum=0;
349     for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
350     l= A->pattern->index[iptrA];
351     p_lj=0;
352 artak 2767 /*This loop can be replaced by bsearch */
353 artak 2759 for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
354     m=P->pattern->index[iptrP];
355     if(m==j) {
356     p_lj=P->val[iptrP];
357     break;
358     }
359     }
360     a_kl=A->val[iptrA];
361     first_sum+=a_kl*p_lj;
362     }
363     r_ik=R->val[iptrR];
364     second_sum+=r_ik*first_sum;
365     }
366     A_c->val[iptrA_c]=second_sum;
367     }
368 artak 2784 }
369     time0=Paso_timer()-time0;
370     if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0);
371 artak 2759
372 artak 2784 Paso_Pattern_free(outpattern);
373     Paso_Pattern_free(temp);
374 artak 2759
375     return A_c;
376     }
377    

  ViewVC Help
Powered by ViewVC 1.1.26