/[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 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 23740 byte(s)
Assorted spelling/comment fixes in paso.

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

  ViewVC Help
Powered by ViewVC 1.1.26