/[escript]/branches/doubleplusgood/paso/src/LocalAMG_Prolongation.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/LocalAMG_Prolongation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 23869 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.26