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

Contents of /trunk/paso/src/LocalAMG_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: 23811 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 the n x n_C 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 strong connection.
45
46 Two settings for P are implemented (see below).
47 */
48
49
50 Paso_SparseMatrix* Paso_Preconditioner_LocalAMG_getProlongation(Paso_SparseMatrix* A_p,
51 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 {
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
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=offset_S[i];
74 for (p=0; p<degree_S[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=offset_S[i];
95 for (p=0; p<degree_S[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 ( (interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING) || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {
119 if (n_block == 1) {
120 Paso_Preconditioner_LocalAMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);
121 } else {
122 Paso_Preconditioner_LocalAMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
123 }
124 } else {
125 if (n_block == 1) {
126 Paso_Preconditioner_LocalAMG_setDirectProlongation(out, A_p, counter_C);
127 } else {
128 Paso_Preconditioner_LocalAMG_setDirectProlongation_Block(out, A_p, counter_C);
129 }
130 }
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 }
141
142 /*
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 void Paso_Preconditioner_LocalAMG_setDirectProlongation(Paso_SparseMatrix* P_p,
158 const Paso_SparseMatrix* A_p,
159 const index_t *counter_C) {
160 dim_t i;
161 const dim_t n =A_p->numRows;
162 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;
163 register index_t iPtr, j, offset;
164 index_t *where_p, *start_p;
165
166 #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 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 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
172 /* 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 strongly connected with j? We search 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 strongly connected 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
209 }
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 rtmp=(-1.)/A_ii;
224 alpha*=rtmp;
225 beta*=rtmp;
226 }
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 }
238
239 void Paso_Preconditioner_LocalAMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p,
240 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 register double A_ij, rtmp;
248 register index_t iPtr, j, offset, ib;
249 index_t *where_p, *start_p;
250
251 #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 {
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
261 #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 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
267 /* 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 strongly connected with j? We search 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 strongly connected 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 rtmp=(-1./A_ii[ib]);
326 alpha[ib]*=rtmp;
327 beta[ib]*=rtmp;
328 }
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 }
352
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 void Paso_Preconditioner_LocalAMG_setClassicProlongation(Paso_SparseMatrix* P_p,
369 Paso_SparseMatrix* A_p,
370 const index_t* offset_S, const dim_t* degree_S, const index_t* S,
371 const index_t *counter_C) {
372 dim_t i, q;
373 const dim_t n =A_p->numRows;
374 double *D_s=NULL;
375 index_t *D_s_offset=NULL, iPtr, iPtr_j;
376 const dim_t ll = Paso_Util_iMax(n, degree_S);
377 const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
378
379
380 #pragma omp parallel private(D_s, D_s_offset, iPtr, q, iPtr_j)
381 {
382 D_s=TMPMEMALLOC(ll,double);
383 D_s_offset=TMPMEMALLOC(ll,index_t);
384
385
386 #pragma omp for private(i) schedule(static)
387 for (i=0;i<n;++i) {
388 if (counter_C[i]>=0) {
389 P_p->val[P_p->pattern->ptr[i]]=1.; /* i is a C row */
390 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {
391 const index_t *start_s = &(S[offset_S[i]]);
392 const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);
393 const dim_t degree_P_i = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];
394 /* this loop sums up the weak connections in A and creates a
395 * list of the strong connected columns which are not in C
396 * (=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 accumulated */
407 a+=A_ij;
408 } else { /* yes i strongly connected 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 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_LocalAMG_setClassicProlongation: Interpolation point is missing.");
413 } 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 } /* end of row i loop */
460 TMPMEMFREE(D_s);
461 TMPMEMFREE(D_s_offset);
462 } /* end of parallel region */
463 }
464
465 void Paso_Preconditioner_LocalAMG_setClassicProlongation_Block(Paso_SparseMatrix* P_p,
466 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
495 * list of the strong connected columns which are not in C
496 * (=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_LocalAMG_setClassicProlongation_Block: 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 }
575

  ViewVC Help
Powered by ViewVC 1.1.26