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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3312 - (show annotations)
Tue Oct 26 07:54:58 2010 UTC (10 years, 5 months ago) by gross
File MIME type: text/plain
File size: 10589 byte(s)
last step for a clean up version of the AMG
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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: defines AMG prolongation */
18
19 /**************************************************************/
20
21 /* Author: Artak Amirbekyan, artak@uq.edu.au, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "SparseMatrix.h"
27 #include "PasoUtil.h"
28 #include "Preconditioner.h"
29
30 /**************************************************************
31
32 Methods nessecary for AMG preconditioner
33
34 construct n x n_C the prolongation matrix P from A_p.
35
36 the columns in A_p to be considered ar marked by counter_C[n] where
37 an unknown i to be considered in P is marked by 0<= counter_C[i] < n_C
38 gives the new column number in P. S defines the strong connections.
39
40 If row i is in C, then P[i,j]=counter_C[i] if counter_C[i]>=0 and 0 otherwise
41 If row i is not C, then P[i,j] = - a[i] * A[i,k]/A[i,i] with j=counter_C[k]>=0 and k in S
42
43 and a[i]=
44 alpha[i] = sum_s min(A[i,s],0)/(sum_{s in S and C} min(A[i,s],0)) A[i,k]<0
45 or if
46 beta[i] = sum_s max(A[i,s],0)/(sum_{s in S and C} max(A[i,s],0)) A[i,k]>0
47
48
49 */
50
51 Paso_SparseMatrix* Paso_Preconditioner_AMG_getDirectProlongation(const Paso_SparseMatrix* A_p,
52 const dim_t* degree, const index_t* S,
53 const dim_t n_C, const index_t* counter_C)
54 {
55 Paso_SparseMatrix* out=NULL;
56 Paso_Pattern *outpattern=NULL;
57 const dim_t n_block=A_p->row_block_size;
58 index_t *ptr=NULL, *index=NULL,j, iptr;
59 const dim_t n =A_p->numRows;
60 dim_t i,p,z, len_P;
61
62 ptr=MEMALLOC(n+1,index_t);
63 if (! Esys_checkPtr(ptr)) {
64
65 /* count the number of entries per row in the Prolongation matrix :*/
66
67 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static)
68 for (i=0;i<n;++i) {
69 if (counter_C[i]>=0) {
70 z=1; /* i is a C unknown */
71 } else {
72 z=0;
73 iptr=A_p->pattern->ptr[i];
74 for (p=0; p<degree[i]; ++p) {
75 j=S[iptr+p]; /* this is a strong connection */
76 if (counter_C[j]>=0) z++; /* and is in C */
77 }
78 }
79 ptr[i]=z;
80 }
81 len_P=Paso_Util_cumsum(n,ptr);
82 ptr[n]=len_P;
83
84 /* allocate and create index vector for prolongation: */
85 index=MEMALLOC(len_P,index_t);
86
87 if (! Esys_checkPtr(index)) {
88 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static)
89 for (i=0;i<n;++i) {
90 if (counter_C[i]>=0) {
91 index[ptr[i]]=counter_C[i]; /* i is a C unknown */
92 } else {
93 z=0;
94 iptr=A_p->pattern->ptr[i];
95 for (p=0; p<degree[i]; ++p) {
96 j=S[iptr+p]; /* this is a strong connection */
97 if (counter_C[j]>=0) { /* and is in C */
98 index[ptr[i]+z]=counter_C[j];
99 z++; /* and is in C */
100 }
101 }
102 }
103 }
104 }
105 }
106 if (Esys_noError()) {
107 outpattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n,n_C,ptr,index);
108 } else {
109 MEMFREE(ptr);
110 MEMFREE(index);
111 }
112 /* now we need to create a matrix and fill it */
113 if (Esys_noError()) {
114 out=Paso_SparseMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK,outpattern,n_block,n_block,FALSE);
115 }
116
117 if (Esys_noError()) {
118 if (n_block == 1) {
119 Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, counter_C);
120 } else {
121 Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, counter_C);
122 }
123 }
124 Paso_Pattern_free(outpattern);
125 if (Esys_noError()) {
126 return out;
127 } else {
128 Paso_SparseMatrix_free(out);
129 return NULL;
130 }
131
132 }
133
134 void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SparseMatrix* P_p,
135 const Paso_SparseMatrix* A_p,
136 const index_t *counter_C) {
137 dim_t i;
138 const dim_t n =A_p->numRows;
139 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii;
140 register index_t iPtr, j, offset;
141 index_t *where_p, *start_p;
142
143 #pragma omp parallel for private(A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij ) schedule(static)
144 for (i=0;i<n;++i) {
145 if (counter_C[i]>=0) {
146 offset = P_p->pattern->ptr[i];
147 P_p->val[offset]=1.; /* i is a C row */
148 } else {
149 /* if i is an F row we first calculate alpha and beta: */
150 sum_all_neg=0; /* sum of all negative values in row i of A */
151 sum_all_pos=0; /* sum of all positive values in row i of A */
152 sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
153 sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
154 A_ii=0;
155 for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
156 j=A_p->pattern->index[iPtr];
157 A_ij=A_p->val[iPtr];
158 if(j==i) {
159 A_ii=A_ij;
160 } else {
161
162 if(A_ij< 0) {
163 sum_all_neg+=A_ij;
164 } else {
165 sum_all_pos+=A_ij;
166 }
167
168 if (counter_C[j]>=0) {
169 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
170 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
171 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
172 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
173 sizeof(index_t),
174 Paso_comparIndex);
175 if (! (where_p == NULL) ) { /* yes i stronly connect with j */
176 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
177 P_p->val[offset]=A_ij; /* will be modified later */
178 if (A_ij< 0) {
179 sum_strong_neg+=A_ij;
180 } else {
181 sum_strong_pos+=A_ij;
182 }
183 }
184 }
185
186 }
187 }
188 if(sum_strong_neg<0) {
189 alpha= sum_all_neg/sum_strong_neg;
190 } else {
191 alpha=0;
192 }
193 if(sum_strong_pos>0) {
194 beta= sum_all_pos/sum_strong_pos;
195 } else {
196 beta=0;
197 A_ii+=sum_all_pos;
198 }
199 if ( A_ii > 0.) {
200 alpha*=(-1./A_ii);
201 beta*=(-1./A_ii);
202 }
203
204 for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
205 A_ij=P_p->val[iPtr];
206 if (A_ij > 0 ) {
207 P_p->val[iPtr]*=beta;
208 } else {
209 P_p->val[iPtr]*=alpha;
210 }
211 }
212 }
213 }
214 }
215
216 void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p,
217 const Paso_SparseMatrix* A_p,
218 const index_t *counter_C) {
219 dim_t i;
220 const dim_t n =A_p->numRows;
221 const dim_t row_block=A_p->row_block_size;
222 const dim_t A_block = A_p->block_size;
223 double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
224 register double A_ij;
225 register index_t iPtr, j, offset, ib;
226 index_t *where_p, *start_p;
227
228 #pragma omp parallel private(ib, A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij )
229 {
230 sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */
231 sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */
232 sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/
233 sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/
234 alpha=TMPMEMALLOC(row_block, double);
235 beta=TMPMEMALLOC(row_block, double);
236 A_ii=TMPMEMALLOC(row_block, double);
237
238 #pragma omp for schedule(static)
239 for (i=0;i<n;++i) {
240 if (counter_C[i]>=0) {
241 offset = P_p->pattern->ptr[i];
242 for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.; /* i is a C row */
243 } else {
244 /* if i is an F row we first calculate alpha and beta: */
245 for (ib =0; ib<row_block; ++ib) {
246 sum_all_neg[ib]=0;
247 sum_all_pos[ib]=0;
248 sum_strong_neg[ib]=0;
249 sum_strong_pos[ib]=0;
250 A_ii[ib]=0;
251 }
252 for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
253 j=A_p->pattern->index[iPtr];
254 if(j==i) {
255 for (ib =0; ib<row_block; ++ib) A_ii[ib]=A_p->val[A_block*iPtr+ib+row_block*ib];
256 } else {
257 for (ib =0; ib<row_block; ++ib) {
258 A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
259 if(A_ij< 0) {
260 sum_all_neg[ib]+=A_ij;
261 } else {
262 sum_all_pos[ib]+=A_ij;
263 }
264 }
265
266 if (counter_C[j]>=0) {
267 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
268 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
269 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
270 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
271 sizeof(index_t),
272 Paso_comparIndex);
273 if (! (where_p == NULL) ) { /* yes i stronly connect with j */
274 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
275 for (ib =0; ib<row_block; ++ib) {
276 A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
277 P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */
278 if (A_ij< 0) {
279 sum_strong_neg[ib]+=A_ij;
280 } else {
281 sum_strong_pos[ib]+=A_ij;
282 }
283 }
284 }
285 }
286
287 }
288 }
289 for (ib =0; ib<row_block; ++ib) {
290 if(sum_strong_neg[ib]<0) {
291 alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];
292 } else {
293 alpha[ib]=0;
294 }
295 if(sum_strong_pos[ib]>0) {
296 beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib];
297 } else {
298 beta[ib]=0;
299 A_ii[ib]+=sum_all_pos[ib];
300 }
301 if ( A_ii[ib] > 0.) {
302 alpha[ib]*=(-1./A_ii[ib]);
303 beta[ib]*=(-1./A_ii[ib]);
304 }
305 }
306
307 for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
308 for (ib =0; ib<row_block; ++ib) {
309 A_ij=P_p->val[row_block*iPtr+ib];
310 if (A_ij > 0 ) {
311 P_p->val[row_block*iPtr+ib]*=beta[ib];
312 } else {
313 P_p->val[row_block*iPtr+ib]*=alpha[ib];
314 }
315 }
316 }
317 }
318 }/* end i loop */
319 TMPMEMFREE(sum_all_neg);
320 TMPMEMFREE(sum_all_pos);
321 TMPMEMFREE(sum_strong_neg);
322 TMPMEMFREE(sum_strong_pos);
323 TMPMEMFREE(alpha);
324 TMPMEMFREE(beta);
325 TMPMEMFREE(A_ii);
326 } /* end parallel region */
327 }

  ViewVC Help
Powered by ViewVC 1.1.26