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

Contents of /trunk/paso/src/LocalAMG_Prolongation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4829 - (show annotations)
Thu Apr 3 04:02:53 2014 UTC (5 years, 5 months ago) by caltinay
File size: 26070 byte(s)
checkpointing some SparseMatrix cleanup.

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

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/LocalAMG_Prolongation.cpp:3531-3826 /branches/lapack2681/paso/src/LocalAMG_Prolongation.cpp:2682-2741 /branches/pasowrap/paso/src/LocalAMG_Prolongation.cpp:3661-3674 /branches/py3_attempt2/paso/src/LocalAMG_Prolongation.cpp:3871-3891 /branches/restext/paso/src/LocalAMG_Prolongation.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/LocalAMG_Prolongation.cpp:3669-3791 /branches/stage3.0/paso/src/LocalAMG_Prolongation.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/LocalAMG_Prolongation.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/LocalAMG_Prolongation.cpp:3517-3974 /release/3.0/paso/src/LocalAMG_Prolongation.cpp:2591-2601 /trunk/paso/src/LocalAMG_Prolongation.cpp:4257-4344 /trunk/ripley/test/python/paso/src/LocalAMG_Prolongation.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26