/[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 3323 - (show annotations)
Thu Oct 28 09:53:46 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 10593 byte(s)
some openmp fizes.
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
66 /* count the number of entries per row in the Prolongation matrix :*/
67
68 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static)
69 for (i=0;i<n;++i) {
70 if (counter_C[i]>=0) {
71 z=1; /* i is a C unknown */
72 } else {
73 z=0;
74 iptr=A_p->pattern->ptr[i];
75 for (p=0; p<degree[i]; ++p) {
76 j=S[iptr+p]; /* this is a strong connection */
77 if (counter_C[j]>=0) z++; /* and is in C */
78 }
79 }
80 ptr[i]=z;
81 }
82 len_P=Paso_Util_cumsum(n,ptr);
83 ptr[n]=len_P;
84
85 /* allocate and create index vector for prolongation: */
86 index=MEMALLOC(len_P,index_t);
87
88 if (! Esys_checkPtr(index)) {
89 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static)
90 for (i=0;i<n;++i) {
91 if (counter_C[i]>=0) {
92 index[ptr[i]]=counter_C[i]; /* i is a C unknown */
93 } else {
94 z=0;
95 iptr=A_p->pattern->ptr[i];
96 for (p=0; p<degree[i]; ++p) {
97 j=S[iptr+p]; /* this is a strong connection */
98 if (counter_C[j]>=0) { /* and is in C */
99 index[ptr[i]+z]=counter_C[j];
100 z++; /* and is in C */
101 }
102 }
103 }
104 }
105 }
106 }
107 if (Esys_noError()) {
108 outpattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n,n_C,ptr,index);
109 } else {
110 MEMFREE(ptr);
111 MEMFREE(index);
112 }
113 /* now we need to create a matrix and fill it */
114 if (Esys_noError()) {
115 out=Paso_SparseMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK,outpattern,n_block,n_block,FALSE);
116 }
117
118 if (Esys_noError()) {
119 if (n_block == 1) {
120 Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, counter_C);
121 } else {
122 Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, counter_C);
123 }
124 }
125 Paso_Pattern_free(outpattern);
126 if (Esys_noError()) {
127 return out;
128 } else {
129 Paso_SparseMatrix_free(out);
130 return NULL;
131 }
132
133 }
134
135 void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SparseMatrix* P_p,
136 const Paso_SparseMatrix* A_p,
137 const index_t *counter_C) {
138 dim_t i;
139 const dim_t n =A_p->numRows;
140 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii;
141 register index_t iPtr, j, offset;
142 index_t *where_p, *start_p;
143
144 #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)
145 for (i=0;i<n;++i) {
146 if (counter_C[i]>=0) {
147 offset = P_p->pattern->ptr[i];
148 P_p->val[offset]=1.; /* i is a C row */
149 } else {
150 /* if i is an F row we first calculate alpha and beta: */
151 sum_all_neg=0; /* sum of all negative values in row i of A */
152 sum_all_pos=0; /* sum of all positive values in row i of A */
153 sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
154 sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
155 A_ii=0;
156 for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
157 j=A_p->pattern->index[iPtr];
158 A_ij=A_p->val[iPtr];
159 if(j==i) {
160 A_ii=A_ij;
161 } else {
162
163 if(A_ij< 0) {
164 sum_all_neg+=A_ij;
165 } else {
166 sum_all_pos+=A_ij;
167 }
168
169 if (counter_C[j]>=0) {
170 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
171 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
172 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
173 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
174 sizeof(index_t),
175 Paso_comparIndex);
176 if (! (where_p == NULL) ) { /* yes i stronly connect with j */
177 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
178 P_p->val[offset]=A_ij; /* will be modified later */
179 if (A_ij< 0) {
180 sum_strong_neg+=A_ij;
181 } else {
182 sum_strong_pos+=A_ij;
183 }
184 }
185 }
186
187 }
188 }
189 if(sum_strong_neg<0) {
190 alpha= sum_all_neg/sum_strong_neg;
191 } else {
192 alpha=0;
193 }
194 if(sum_strong_pos>0) {
195 beta= sum_all_pos/sum_strong_pos;
196 } else {
197 beta=0;
198 A_ii+=sum_all_pos;
199 }
200 if ( A_ii > 0.) {
201 alpha*=(-1./A_ii);
202 beta*=(-1./A_ii);
203 }
204
205 for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
206 A_ij=P_p->val[iPtr];
207 if (A_ij > 0 ) {
208 P_p->val[iPtr]*=beta;
209 } else {
210 P_p->val[iPtr]*=alpha;
211 }
212 }
213 }
214 }
215 }
216
217 void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p,
218 const Paso_SparseMatrix* A_p,
219 const index_t *counter_C) {
220 dim_t i;
221 const dim_t n =A_p->numRows;
222 const dim_t row_block=A_p->row_block_size;
223 const dim_t A_block = A_p->block_size;
224 double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
225 register double A_ij;
226 register index_t iPtr, j, offset, ib;
227 index_t *where_p, *start_p;
228
229 #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 )
230 {
231 sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */
232 sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */
233 sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/
234 sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/
235 alpha=TMPMEMALLOC(row_block, double);
236 beta=TMPMEMALLOC(row_block, double);
237 A_ii=TMPMEMALLOC(row_block, double);
238
239 #pragma omp for schedule(static)
240 for (i=0;i<n;++i) {
241 if (counter_C[i]>=0) {
242 offset = P_p->pattern->ptr[i];
243 for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.; /* i is a C row */
244 } else {
245 /* if i is an F row we first calculate alpha and beta: */
246 for (ib =0; ib<row_block; ++ib) {
247 sum_all_neg[ib]=0;
248 sum_all_pos[ib]=0;
249 sum_strong_neg[ib]=0;
250 sum_strong_pos[ib]=0;
251 A_ii[ib]=0;
252 }
253 for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
254 j=A_p->pattern->index[iPtr];
255 if(j==i) {
256 for (ib =0; ib<row_block; ++ib) A_ii[ib]=A_p->val[A_block*iPtr+ib+row_block*ib];
257 } else {
258 for (ib =0; ib<row_block; ++ib) {
259 A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
260 if(A_ij< 0) {
261 sum_all_neg[ib]+=A_ij;
262 } else {
263 sum_all_pos[ib]+=A_ij;
264 }
265 }
266
267 if (counter_C[j]>=0) {
268 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
269 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
270 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
271 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
272 sizeof(index_t),
273 Paso_comparIndex);
274 if (! (where_p == NULL) ) { /* yes i stronly connect with j */
275 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
276 for (ib =0; ib<row_block; ++ib) {
277 A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
278 P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */
279 if (A_ij< 0) {
280 sum_strong_neg[ib]+=A_ij;
281 } else {
282 sum_strong_pos[ib]+=A_ij;
283 }
284 }
285 }
286 }
287
288 }
289 }
290 for (ib =0; ib<row_block; ++ib) {
291 if(sum_strong_neg[ib]<0) {
292 alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];
293 } else {
294 alpha[ib]=0;
295 }
296 if(sum_strong_pos[ib]>0) {
297 beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib];
298 } else {
299 beta[ib]=0;
300 A_ii[ib]+=sum_all_pos[ib];
301 }
302 if ( A_ii[ib] > 0.) {
303 alpha[ib]*=(-1./A_ii[ib]);
304 beta[ib]*=(-1./A_ii[ib]);
305 }
306 }
307
308 for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
309 for (ib =0; ib<row_block; ++ib) {
310 A_ij=P_p->val[row_block*iPtr+ib];
311 if (A_ij > 0 ) {
312 P_p->val[row_block*iPtr+ib]*=beta[ib];
313 } else {
314 P_p->val[row_block*iPtr+ib]*=alpha[ib];
315 }
316 }
317 }
318 }
319 }/* end i loop */
320 TMPMEMFREE(sum_all_neg);
321 TMPMEMFREE(sum_all_pos);
322 TMPMEMFREE(sum_strong_neg);
323 TMPMEMFREE(sum_strong_pos);
324 TMPMEMFREE(alpha);
325 TMPMEMFREE(beta);
326 TMPMEMFREE(A_ii);
327 } /* end parallel region */
328 }

  ViewVC Help
Powered by ViewVC 1.1.26