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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3441 - (hide annotations)
Fri Jan 14 01:09:09 2011 UTC (9 years, 10 months ago) by gross
File MIME type: text/plain
File size: 23880 byte(s)
clarifiaction of function names as a first step towards a MPI version of AMG
1 artak 2759
2     /*******************************************************
3     *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 artak 2759 * 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 gross 3303 /* Paso: defines AMG prolongation */
18 artak 2759
19     /**************************************************************/
20    
21 gross 3303 /* Author: Artak Amirbekyan, artak@uq.edu.au, l.gross@uq.edu.au */
22 artak 2759
23     /**************************************************************/
24    
25     #include "Paso.h"
26     #include "SparseMatrix.h"
27     #include "PasoUtil.h"
28 gross 3303 #include "Preconditioner.h"
29 artak 2759
30     /**************************************************************
31    
32     Methods nessecary for AMG preconditioner
33    
34 gross 3303 construct n x n_C the prolongation matrix P from A_p.
35    
36 gross 3440 the columns in A_p to be considered are marked by counter_C[n] where
37     an unknown i is to be considered in P is marked by 0<= counter_C[i] < n_C
38     and counter_C[i] gives the new column number in P. S defines the strong connections.
39 gross 3303
40 gross 3440 The pattern of P is formed as follows:
41    
42     If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise
43     If row i is not C, then P[i,j] <> 0 if counter_C[k]==j (k in C) and (i,k) is strong connection.
44    
45     two settings for P are implemented (see below)
46 gross 3303
47 gross 3440 */
48 artak 2759
49    
50 gross 3441 Paso_SparseMatrix* Paso_Preconditioner_LocalAMG_getProlongation(Paso_SparseMatrix* A_p,
51 gross 3440 const index_t* offset_S, const dim_t* degree_S, const index_t* S,
52     const dim_t n_C, const index_t* counter_C, const index_t interpolation_method)
53 gross 3303 {
54     Paso_SparseMatrix* out=NULL;
55     Paso_Pattern *outpattern=NULL;
56     const dim_t n_block=A_p->row_block_size;
57     index_t *ptr=NULL, *index=NULL,j, iptr;
58     const dim_t n =A_p->numRows;
59     dim_t i,p,z, len_P;
60    
61     ptr=MEMALLOC(n+1,index_t);
62     if (! Esys_checkPtr(ptr)) {
63 gross 3323
64 gross 3303
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 gross 3440 iptr=offset_S[i];
74     for (p=0; p<degree_S[i]; ++p) {
75 gross 3303 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 artak 2759 }
81 gross 3303 len_P=Paso_Util_cumsum(n,ptr);
82     ptr[n]=len_P;
83 gross 3323
84 gross 3303 /* 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 gross 3440 iptr=offset_S[i];
95     for (p=0; p<degree_S[i]; ++p) {
96 gross 3303 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 artak 2759 }
105 gross 3303 }
106     if (Esys_noError()) {
107 gross 3312 outpattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n,n_C,ptr,index);
108 gross 3303 } 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 gross 3440 if ( (interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING) || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {
119     if (n_block == 1) {
120 gross 3441 Paso_Preconditioner_LocalAMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);
121 gross 3440 } else {
122 gross 3441 Paso_Preconditioner_LocalAMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
123 gross 3440 }
124 gross 3303 } else {
125 gross 3440 if (n_block == 1) {
126 gross 3441 Paso_Preconditioner_LocalAMG_setDirectProlongation(out, A_p, counter_C);
127 gross 3440 } else {
128 gross 3441 Paso_Preconditioner_LocalAMG_setDirectProlongation_Block(out, A_p, counter_C);
129 gross 3440 }
130 gross 3303 }
131     }
132     Paso_Pattern_free(outpattern);
133     if (Esys_noError()) {
134     return out;
135     } else {
136     Paso_SparseMatrix_free(out);
137     return NULL;
138     }
139    
140 artak 2759 }
141    
142 gross 3440 /*
143     Direct Prolongation:
144     -------------------
145    
146     If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise
147     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
148    
149     and a[i]=
150     alpha[i] = sum_s min(A[i,s],0)/(sum_{s in S and C} min(A[i,s],0)) A[i,k]<0
151     or if
152     beta[i] = sum_s max(A[i,s],0)/(sum_{s in S and C} max(A[i,s],0)) A[i,k]>0
153    
154    
155     */
156    
157 gross 3441 void Paso_Preconditioner_LocalAMG_setDirectProlongation(Paso_SparseMatrix* P_p,
158 gross 3303 const Paso_SparseMatrix* A_p,
159     const index_t *counter_C) {
160     dim_t i;
161     const dim_t n =A_p->numRows;
162 gross 3402 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;
163 gross 3303 register index_t iPtr, j, offset;
164     index_t *where_p, *start_p;
165 artak 3010
166 gross 3402 #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 , rtmp) schedule(static)
167 gross 3303 for (i=0;i<n;++i) {
168     if (counter_C[i]>=0) {
169     offset = P_p->pattern->ptr[i];
170     P_p->val[offset]=1.; /* i is a C row */
171 gross 3440 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
172 gross 3303 /* if i is an F row we first calculate alpha and beta: */
173     sum_all_neg=0; /* sum of all negative values in row i of A */
174     sum_all_pos=0; /* sum of all positive values in row i of A */
175     sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
176     sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
177     A_ii=0;
178     for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
179     j=A_p->pattern->index[iPtr];
180     A_ij=A_p->val[iPtr];
181     if(j==i) {
182     A_ii=A_ij;
183     } else {
184    
185     if(A_ij< 0) {
186     sum_all_neg+=A_ij;
187     } else {
188     sum_all_pos+=A_ij;
189     }
190    
191     if (counter_C[j]>=0) {
192     /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
193     start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
194     where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
195     P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
196     sizeof(index_t),
197     Paso_comparIndex);
198     if (! (where_p == NULL) ) { /* yes i stronly connect with j */
199     offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
200     P_p->val[offset]=A_ij; /* will be modified later */
201     if (A_ij< 0) {
202     sum_strong_neg+=A_ij;
203     } else {
204     sum_strong_pos+=A_ij;
205     }
206     }
207     }
208 artak 2759
209 gross 3303 }
210     }
211     if(sum_strong_neg<0) {
212     alpha= sum_all_neg/sum_strong_neg;
213     } else {
214     alpha=0;
215     }
216     if(sum_strong_pos>0) {
217     beta= sum_all_pos/sum_strong_pos;
218     } else {
219     beta=0;
220     A_ii+=sum_all_pos;
221     }
222     if ( A_ii > 0.) {
223 gross 3402 rtmp=(-1.)/A_ii;
224     alpha*=rtmp;
225     beta*=rtmp;
226 gross 3303 }
227     for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
228     A_ij=P_p->val[iPtr];
229     if (A_ij > 0 ) {
230     P_p->val[iPtr]*=beta;
231     } else {
232     P_p->val[iPtr]*=alpha;
233     }
234     }
235     }
236     }
237 artak 2759 }
238    
239 gross 3441 void Paso_Preconditioner_LocalAMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p,
240 gross 3303 const Paso_SparseMatrix* A_p,
241     const index_t *counter_C) {
242     dim_t i;
243     const dim_t n =A_p->numRows;
244     const dim_t row_block=A_p->row_block_size;
245     const dim_t A_block = A_p->block_size;
246     double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
247 gross 3402 register double A_ij, rtmp;
248 gross 3303 register index_t iPtr, j, offset, ib;
249     index_t *where_p, *start_p;
250 artak 3010
251 gross 3402 #pragma omp parallel private(ib, rtmp, 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 )
252 gross 3303 {
253     sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */
254     sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */
255     sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/
256     sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/
257     alpha=TMPMEMALLOC(row_block, double);
258     beta=TMPMEMALLOC(row_block, double);
259     A_ii=TMPMEMALLOC(row_block, double);
260 artak 2759
261 gross 3303 #pragma omp for schedule(static)
262     for (i=0;i<n;++i) {
263     if (counter_C[i]>=0) {
264     offset = P_p->pattern->ptr[i];
265     for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.; /* i is a C row */
266 gross 3440 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
267 gross 3303 /* if i is an F row we first calculate alpha and beta: */
268     for (ib =0; ib<row_block; ++ib) {
269     sum_all_neg[ib]=0;
270     sum_all_pos[ib]=0;
271     sum_strong_neg[ib]=0;
272     sum_strong_pos[ib]=0;
273     A_ii[ib]=0;
274     }
275     for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
276     j=A_p->pattern->index[iPtr];
277     if(j==i) {
278     for (ib =0; ib<row_block; ++ib) A_ii[ib]=A_p->val[A_block*iPtr+ib+row_block*ib];
279     } else {
280     for (ib =0; ib<row_block; ++ib) {
281     A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
282     if(A_ij< 0) {
283     sum_all_neg[ib]+=A_ij;
284     } else {
285     sum_all_pos[ib]+=A_ij;
286     }
287     }
288    
289     if (counter_C[j]>=0) {
290     /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
291     start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);
292     where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
293     P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],
294     sizeof(index_t),
295     Paso_comparIndex);
296     if (! (where_p == NULL) ) { /* yes i stronly connect with j */
297     offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
298     for (ib =0; ib<row_block; ++ib) {
299     A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];
300     P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */
301     if (A_ij< 0) {
302     sum_strong_neg[ib]+=A_ij;
303     } else {
304     sum_strong_pos[ib]+=A_ij;
305     }
306     }
307     }
308     }
309    
310     }
311     }
312     for (ib =0; ib<row_block; ++ib) {
313     if(sum_strong_neg[ib]<0) {
314     alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];
315     } else {
316     alpha[ib]=0;
317     }
318     if(sum_strong_pos[ib]>0) {
319     beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib];
320     } else {
321     beta[ib]=0;
322     A_ii[ib]+=sum_all_pos[ib];
323     }
324     if ( A_ii[ib] > 0.) {
325 gross 3402 rtmp=(-1./A_ii[ib]);
326     alpha[ib]*=rtmp;
327     beta[ib]*=rtmp;
328 gross 3303 }
329     }
330    
331     for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
332     for (ib =0; ib<row_block; ++ib) {
333     A_ij=P_p->val[row_block*iPtr+ib];
334     if (A_ij > 0 ) {
335     P_p->val[row_block*iPtr+ib]*=beta[ib];
336     } else {
337     P_p->val[row_block*iPtr+ib]*=alpha[ib];
338     }
339     }
340     }
341     }
342     }/* end i loop */
343     TMPMEMFREE(sum_all_neg);
344     TMPMEMFREE(sum_all_pos);
345     TMPMEMFREE(sum_strong_neg);
346     TMPMEMFREE(sum_strong_pos);
347     TMPMEMFREE(alpha);
348     TMPMEMFREE(beta);
349     TMPMEMFREE(A_ii);
350     } /* end parallel region */
351 gross 3312 }
352 gross 3440
353     /*
354     Classic Prolongation:
355     -------------------
356    
357     If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise
358     If row i is not C, then P[i,j] = - 1/a[i] * ( A[i,k] + sum_{l} A[i,l]*A+[l,k]/B[i,k])
359     where the summation over l is considering columns which are strongly connected
360     to i (l in S[i]) and not in C (counter_C[l]<0) and
361    
362     B[i,k]=sum_{m in S_i and in C} A+[k,m]
363     a[i]=A[i,i]+sum{l not strongly connected to i} A[i,l]
364    
365     A+[i,k]=A[i,k] if sign(A[i,k])==sign(A[i,i]) or 0 otherwise
366    
367    
368     */
369 gross 3441 void Paso_Preconditioner_LocalAMG_setClassicProlongation(Paso_SparseMatrix* P_p,
370 gross 3440 Paso_SparseMatrix* A_p,
371     const index_t* offset_S, const dim_t* degree_S, const index_t* S,
372     const index_t *counter_C) {
373     dim_t i, q;
374     const dim_t n =A_p->numRows;
375     double *D_s=NULL;
376     index_t *D_s_offset=NULL, iPtr, iPtr_j;
377     const dim_t ll = Paso_Util_iMax(n, degree_S);
378     const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
379    
380    
381     #pragma omp parallel private(D_s, D_s_offset, iPtr, q, iPtr_j)
382     {
383     D_s=TMPMEMALLOC(ll,double);
384     D_s_offset=TMPMEMALLOC(ll,index_t);
385    
386    
387     #pragma omp for private(i) schedule(static)
388     for (i=0;i<n;++i) {
389     if (counter_C[i]>=0) {
390     P_p->val[P_p->pattern->ptr[i]]=1.; /* i is a C row */
391     } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
392     const index_t *start_s = &(S[offset_S[i]]);
393     const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
394     const dim_t degree_P_i = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
395     /* this loop sums up the weak connections in a and creates a list of the strong connected columns
396     which are not in C (=no interpolation nodes) */
397     const double A_ii = A_p->val[ptr_main_A[i]];
398     double a=A_ii;
399    
400     for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
401     const index_t j=A_p->pattern->index[iPtr];
402     const double A_ij=A_p->val[iPtr];
403     if ( (i!=j) && (degree_S[j]>0) ) {
404     /* is (i,j) a strong connection ?*/
405     const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
406     if (where_s == NULL) { /* weak connections are accummulated */
407     a+=A_ij;
408     } else { /* yes i stronly connect with j */
409     if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */
410     const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);
411     if (where_p == NULL) {
412 gross 3441 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_LocalAMG_setBoomerProlongation: interpolation point is missing.");
413 gross 3440 } else {
414     const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
415     P_p->val[offset]+=A_ij;
416     }
417     } else { /* j is not an interpolation point */
418     /* find all interpolation points m of k */
419     const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
420     const dim_t degree_P_j = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
421     double s=0.;
422     dim_t len_D_s=0;
423     for (iPtr_j=A_p->pattern->ptr[j];iPtr_j<A_p->pattern->ptr[j + 1]; ++iPtr_j) {
424     const double A_jm=A_p->val[iPtr_j];
425     const index_t m=A_p->pattern->index[iPtr_j];
426     /* is m an interpolation point ? */
427     const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex);
428     if (! (where_p_m==NULL)) {
429     const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j);
430     if (! SAMESIGN(A_ii,A_jm)) {
431     D_s[len_D_s]=A_jm;
432     } else {
433     D_s[len_D_s]=0.;
434     }
435     D_s_offset[len_D_s]=offset_m;
436     len_D_s++;
437     }
438     }
439     for (q=0;q<len_D_s;++q) s+=D_s[q];
440     if (ABS(s)>0) {
441     s=A_ij/s;
442     for (q=0;q<len_D_s;++q) {
443     P_p->val[D_s_offset[q]]+=s*D_s[q];
444     }
445     } else {
446     a+=A_ij;
447     }
448     }
449     }
450     }
451     } /* i has been processed, now we need to do some rescaling */
452     if (ABS(a)>0.) {
453     a=-1./a;
454     for (iPtr=P_p->pattern->ptr[i]; iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
455     P_p->val[iPtr]*=a;
456     }
457     }
458     }
459     } /* endo of row i loop */
460     TMPMEMFREE(D_s);
461     TMPMEMFREE(D_s_offset);
462     } /* end of parallel region */
463     }
464    
465 gross 3441 void Paso_Preconditioner_LocalAMG_setClassicProlongation_Block(Paso_SparseMatrix* P_p,
466 gross 3440 Paso_SparseMatrix* A_p,
467     const index_t* offset_S, const dim_t* degree_S, const index_t* S,
468     const index_t *counter_C) {
469     dim_t i, q, ib;
470     const dim_t row_block=A_p->row_block_size;
471     const dim_t A_block = A_p->block_size;
472     const dim_t n =A_p->numRows;
473     double *D_s=NULL;
474     index_t *D_s_offset=NULL, iPtr, iPtr_j;
475     const dim_t ll = Paso_Util_iMax(n, degree_S);
476     const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
477    
478    
479     #pragma omp parallel private(D_s, D_s_offset, iPtr, q, iPtr_j,ib)
480     {
481     double *a=TMPMEMALLOC(row_block, double);
482     D_s=TMPMEMALLOC(row_block*ll,double);
483     D_s_offset=TMPMEMALLOC(row_block*ll,index_t);
484    
485     #pragma omp for private(i) schedule(static)
486     for (i=0;i<n;++i) {
487     if (counter_C[i]>=0) {
488     const index_t offset = P_p->pattern->ptr[i];
489     for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.; /* i is a C row */
490     } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
491     const index_t *start_s = &(S[offset_S[i]]);
492     const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
493     const dim_t degree_P_i = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
494     /* this loop sums up the weak connections in a and creates a list of the strong connected columns
495     which are not in C (=no interpolation nodes) */
496     const double *A_ii = &(A_p->val[ptr_main_A[i]*A_block]);
497     for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];
498    
499    
500     for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {
501     const index_t j=A_p->pattern->index[iPtr];
502     const double* A_ij=&(A_p->val[iPtr*A_block]);
503    
504     if ( (i!=j) && (degree_S[j]>0) ) {
505     /* is (i,j) a strong connection ?*/
506     const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
507     if (where_s == NULL) { /* weak connections are accummulated */
508     for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
509     } else { /* yes i stronly connect with j */
510     if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */
511     const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);
512     if (where_p == NULL) {
513 gross 3441 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_LocalAMG_setClassicProlongation: interpolation point is missing.");
514 gross 3440 } else {
515     const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);
516     for (ib=0; ib<row_block; ib++) P_p->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
517     }
518     } else { /* j is not an interpolation point */
519     /* find all interpolation points m of k */
520     const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
521     const dim_t degree_P_j = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
522     dim_t len_D_s=0;
523     for (iPtr_j=A_p->pattern->ptr[j];iPtr_j<A_p->pattern->ptr[j + 1]; ++iPtr_j) {
524     const double* A_jm=&(A_p->val[iPtr_j*A_block]);
525     const index_t m=A_p->pattern->index[iPtr_j];
526     /* is m an interpolation point ? */
527     const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex);
528     if (! (where_p_m==NULL)) {
529     const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j);
530     for (ib=0; ib<row_block; ib++) {
531     if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
532     D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
533     } else {
534     D_s[len_D_s*row_block+ib]=0.;
535     }
536     }
537     D_s_offset[len_D_s]=offset_m;
538     len_D_s++;
539     }
540     }
541     for (ib=0; ib<row_block; ib++) {
542     double s=0;
543     for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
544    
545     if (ABS(s)>0) {
546     s=A_ij[(row_block+1)*ib]/s;
547     for (q=0;q<len_D_s;++q) {
548     P_p->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
549     }
550     } else {
551     a[ib]+=A_ij[(row_block+1)*ib];
552     }
553     }
554     }
555     }
556     }
557     } /* i has been processed, now we need to do some rescaling */
558     for (ib=0; ib<row_block; ib++) {
559     register double a2=a[ib];
560     if (ABS(a2)>0.) {
561     a2=-1./a2;
562     for (iPtr=P_p->pattern->ptr[i]; iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {
563     P_p->val[iPtr*row_block+ib]*=a2;
564     }
565     }
566     }
567     }
568     } /* endo of row i loop */
569     TMPMEMFREE(D_s);
570     TMPMEMFREE(D_s_offset);
571     TMPMEMFREE(a);
572     } /* end of parallel region */
573     }

  ViewVC Help
Powered by ViewVC 1.1.26