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

Contents of /trunk/paso/src/AMG_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: 51336 byte(s)
checkpointing some SparseMatrix cleanup.

1 /*****************************************************************************
2 *
3 * Copyright (c) 2003-2014 by University of Queensland
4 * http://www.uq.edu.au
5 *
6 * Primary Business: Queensland, Australia
7 * Licensed under the Open Software License version 3.0
8 * http://www.opensource.org/licenses/osl-3.0.php
9 *
10 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 * Development 2012-2013 by School of Earth Sciences
12 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
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 n x n_C the 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 a strong connection.
47
48 Two settings for P are implemented (see below)
49
50 */
51
52 Paso_SystemMatrix* Paso_Preconditioner_AMG_getProlongation(
53 Paso_SystemMatrix* A_p, const index_t* offset_S, const dim_t* degree_S,
54 const index_t* S, const dim_t n_C, index_t* counter_C,
55 const index_t interpolation_method)
56 {
57 Esys_MPIInfo *mpi_info=Esys_MPIInfo_getReference(A_p->mpi_info);
58 Paso_SystemMatrix *out=NULL;
59 paso::SystemMatrixPattern_ptr pattern;
60 paso::Distribution_ptr input_dist, output_dist;
61 paso::SharedComponents_ptr send, recv;
62 paso::Connector_ptr col_connector;
63 paso::Pattern_ptr main_pattern, couple_pattern;
64 const dim_t row_block_size=A_p->row_block_size;
65 const dim_t col_block_size=A_p->col_block_size;
66 const dim_t my_n=A_p->mainBlock->numCols;
67 const dim_t overlap_n=A_p->col_coupleBlock->numCols;
68 const dim_t num_threads=omp_get_max_threads();
69 index_t size=mpi_info->size, *dist=NULL;
70 index_t *main_p=NULL, *couple_p=NULL, *main_idx=NULL, *couple_idx=NULL;
71 index_t *shared=NULL, *offsetInShared=NULL;
72 index_t *recv_shared=NULL, *send_shared=NULL;
73 index_t sum, i, j, k, l, p, q, iptr;
74 index_t my_n_C, global_label, num_neighbors;
75 #ifdef ESYS_MPI
76 index_t rank=mpi_info->rank;
77 #endif
78 Esys_MPI_rank *neighbor=NULL;
79 #ifdef ESYS_MPI
80 MPI_Request* mpi_requests=NULL;
81 MPI_Status* mpi_stati=NULL;
82 #else
83 int *mpi_requests=NULL, *mpi_stati=NULL;
84 #endif
85
86 /* number of C points in current distribution */
87 my_n_C = 0;
88 sum=0;
89 if (num_threads>1) {
90 #pragma omp parallel private(i,sum)
91 {
92 sum=0;
93 #pragma omp for schedule(static)
94 for (i=0;i<my_n;++i) {
95 if (counter_C[i] != -1) {
96 sum++;
97 }
98 }
99 #pragma omp critical
100 {
101 my_n_C += sum;
102 }
103 }
104 } else { /* num_threads=1 */
105 for (i=0;i<my_n;++i) {
106 if (counter_C[i] != -1) {
107 my_n_C++;
108 }
109 }
110 }
111
112 /* create row distribution (output_distribution) and col distribution
113 (input_distribution) */
114 /* ??? should I alloc an new Esys_MPIInfo object or reuse the one in
115 system matrix A. for now, I'm reuse A->mpi_info ??? */
116 dist = A_p->pattern->output_distribution->first_component;
117 output_dist.reset(new paso::Distribution(mpi_info, dist, 1, 0));
118 dist = new index_t[size+1]; /* now prepare for col distribution */
119 #ifdef ESYS_MPI
120 MPI_Allgather(&my_n_C, 1, MPI_INT, dist, 1, MPI_INT, mpi_info->comm);
121 #endif
122 global_label=0;
123 for (i=0; i<size; i++) {
124 k = dist[i];
125 dist[i] = global_label;
126 global_label += k;
127 }
128 dist[size] = global_label;
129
130 input_dist.reset(new paso::Distribution(mpi_info, dist, 1, 0));
131 delete[] dist;
132
133 /* create pattern for mainBlock and coupleBlock */
134 main_p = new index_t[my_n+1];
135 couple_p = new index_t[my_n+1];
136 /* count the number of entries per row in the Prolongation matrix :*/
137 #pragma omp parallel for private(i,l,k,iptr,j,p) schedule(static)
138 for (i=0; i<my_n; i++) {
139 l = 0;
140 if (counter_C[i]>=0) {
141 k = 1; /* i is a C unknown */
142 } else {
143 k = 0;
144 iptr = offset_S[i];
145 for (p=0; p<degree_S[i]; p++) {
146 j = S[iptr+p]; /* this is a strong connection */
147 if (counter_C[j]>=0) { /* and is in C */
148 if (j <my_n) k++;
149 else {
150 l++;
151 }
152 }
153 }
154 }
155 main_p[i] = k;
156 couple_p[i] = l;
157 }
158
159 /* number of unknowns in the col-coupleBlock of the interpolation matrix */
160 sum = 0;
161 for (i=0;i<overlap_n;++i) {
162 if (counter_C[i+my_n] > -1) {
163 counter_C[i+my_n] -= my_n_C;
164 sum++;
165 }
166 }
167
168 /* allocate and create index vector for prolongation: */
169 p = Paso_Util_cumsum(my_n, main_p);
170 main_p[my_n] = p;
171 main_idx = new index_t[p];
172 p = Paso_Util_cumsum(my_n, couple_p);
173 couple_p[my_n] = p;
174 couple_idx = new index_t[p];
175 #pragma omp parallel for private(i,k,l,iptr,j,p) schedule(static)
176 for (i=0; i<my_n; i++) {
177 if (counter_C[i]>=0) {
178 main_idx[main_p[i]]=counter_C[i]; /* i is a C unknown */
179 } else {
180 k = 0;
181 l = 0;
182 iptr = offset_S[i];
183 for (p=0; p<degree_S[i]; p++) {
184 j = S[iptr+p]; /* this is a strong connection */
185 if (counter_C[j] >=0) { /* and is in C */
186 if (j < my_n) {
187 main_idx[main_p[i]+k] = counter_C[j];
188 k++;
189 } else {
190 couple_idx[couple_p[i]+l] = counter_C[j];
191 l++;
192 }
193 }
194 }
195 }
196 }
197
198 if (Esys_noError()) {
199 main_pattern.reset(new paso::Pattern(MATRIX_FORMAT_DEFAULT, my_n,
200 my_n_C, main_p, main_idx));
201 couple_pattern.reset(new paso::Pattern(MATRIX_FORMAT_DEFAULT, my_n,
202 sum, couple_p, couple_idx));
203 } else {
204 delete[] main_p;
205 delete[] main_idx;
206 delete[] couple_p;
207 delete[] couple_idx;
208 }
209
210 /* prepare the receiver for the col_connector.
211 Note that the allocation for "shared" assumes the send and receive buffer
212 of the interpolation matrix P is no larger than that of matrix A_p. */
213 neighbor = new Esys_MPI_rank[size];
214 offsetInShared = new index_t[size+1];
215 recv = A_p->col_coupler->connector->recv;
216 send = A_p->col_coupler->connector->send;
217 i = recv->numSharedComponents;
218 recv_shared = new index_t[i];
219 memset(recv_shared, 0, sizeof(index_t)*i);
220 k = send->numSharedComponents;
221 send_shared = new index_t[k];
222 if (k > i) i = k;
223 shared = new index_t[i];
224
225 #ifdef ESYS_MPI
226 mpi_requests=new MPI_Request[size*2];
227 mpi_stati=new MPI_Status[size*2];
228 #else
229 mpi_requests=new int[size*2];
230 mpi_stati=new int[size*2];
231 #endif
232
233 for (p=0; p<send->numNeighbors; p++) {
234 i = send->offsetInShared[p];
235 #ifdef ESYS_MPI
236 MPI_Irecv (&(send_shared[i]), send->offsetInShared[p+1]-i, MPI_INT,
237 send->neighbor[p], mpi_info->msg_tag_counter+send->neighbor[p],
238 mpi_info->comm, &mpi_requests[p]);
239 #endif
240 }
241
242 num_neighbors = 0;
243 q = 0;
244 p = recv->numNeighbors;
245 offsetInShared[0]=0;
246 for (i=0; i<p; i++) {
247 l = 0;
248 k = recv->offsetInShared[i+1];
249 for (j=recv->offsetInShared[i]; j<k; j++) {
250 if (counter_C[recv->shared[j]] > -1) {
251 shared[q] = my_n_C + q;
252 recv_shared[recv->shared[j]-my_n] = 1;
253 q++;
254 l = 1;
255 }
256 }
257 if (l == 1) {
258 iptr = recv->neighbor[i];
259 neighbor[num_neighbors] = iptr;
260 num_neighbors++;
261 offsetInShared[num_neighbors] = q;
262 }
263 #ifdef ESYS_MPI
264 MPI_Issend(&(recv_shared[recv->offsetInShared[i]]),
265 k-recv->offsetInShared[i], MPI_INT, recv->neighbor[i],
266 mpi_info->msg_tag_counter+rank, mpi_info->comm,
267 &mpi_requests[i+send->numNeighbors]);
268 #endif
269 }
270 recv.reset(new paso::SharedComponents(my_n_C, num_neighbors, neighbor,
271 shared, offsetInShared, 1, 0, mpi_info));
272
273 /* now we can build the sender */
274 #ifdef ESYS_MPI
275 MPI_Waitall(recv->numNeighbors+send->numNeighbors, mpi_requests, mpi_stati);
276 #endif
277 ESYS_MPI_INC_COUNTER(*mpi_info, size)
278 delete[] mpi_requests;
279 delete[] mpi_stati;
280
281 num_neighbors = 0;
282 q = 0;
283 p = send->numNeighbors;
284 offsetInShared[0]=0;
285 for (i=0; i<p; i++) {
286 l = 0;
287 k = send->offsetInShared[i+1];
288 for (j=send->offsetInShared[i]; j<k; j++) {
289 if (send_shared[j] == 1) {
290 shared[q] = counter_C[send->shared[j]];
291 q++;
292 l = 1;
293 }
294 }
295 if (l == 1) {
296 iptr = send->neighbor[i];
297 neighbor[num_neighbors] = iptr;
298 num_neighbors++;
299 offsetInShared[num_neighbors] = q;
300 }
301 }
302
303 send.reset(new paso::SharedComponents(my_n_C, num_neighbors, neighbor,
304 shared, offsetInShared, 1, 0, mpi_info));
305 col_connector.reset(new paso::Connector(send, recv));
306 delete[] recv_shared;
307 delete[] send_shared;
308 delete[] neighbor;
309 delete[] offsetInShared;
310 delete[] shared;
311
312 /* now we need to create the System Matrix
313 TO BE FIXED: at this stage, we only construction col_couple_pattern
314 and col_connector for interpolation matrix P. To be completed,
315 row_couple_pattern and row_connector need to be constructed as well */
316 if (Esys_noError()) {
317 pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
318 output_dist, input_dist, main_pattern, couple_pattern,
319 couple_pattern, col_connector, col_connector));
320 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK, pattern,
321 row_block_size, col_block_size, FALSE);
322 }
323
324 /* now fill in the matrix */
325 if (Esys_noError()) {
326 if ((interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
327 || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {
328 if (row_block_size == 1) {
329 Paso_Preconditioner_AMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);
330 } else {
331 Paso_Preconditioner_AMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
332 }
333 } else {
334 if (row_block_size == 1) {
335 Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, offset_S, degree_S, S, counter_C);
336 } else {
337 Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
338 }
339 }
340 }
341
342 /* clean up */
343 if (Esys_noError()) {
344 return out;
345 } else {
346 Paso_SystemMatrix_free(out);
347 return NULL;
348 }
349
350 }
351
352 /*
353 Direct Prolongation:
354 -------------------
355
356 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
357 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
358
359 and a[i]=
360 alpha[i] = sum_s min(A[i,s],0)/(sum_{s in S and C} min(A[i,s],0)) A[i,k]<0
361 or if
362 beta[i] = sum_s max(A[i,s],0)/(sum_{s in S and C} max(A[i,s],0)) A[i,k]>0
363
364
365 */
366
367 void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SystemMatrix* P,
368 Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
369 const index_t* S, const index_t *counter_C)
370 {
371 paso::SparseMatrix_ptr main_block(P->mainBlock);
372 paso::SparseMatrix_ptr couple_block(P->col_coupleBlock);
373 paso::Pattern_ptr main_pattern(main_block->pattern);
374 paso::Pattern_ptr couple_pattern(couple_block->pattern);
375 const dim_t my_n=A->mainBlock->numRows;
376 index_t range;
377
378 dim_t i;
379 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;
380 register index_t iPtr, j, offset;
381 index_t *where_p, *start_p;
382
383 #pragma omp parallel for private(i,offset,sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos,A_ii,range,iPtr,j,A_ij,start_p,where_p,alpha,beta,rtmp) schedule(static)
384 for (i=0; i<my_n; i++) {
385 if (counter_C[i]>=0) {
386 offset = main_pattern->ptr[i];
387 main_block->val[offset]=1.; /* i is a C row */
388 } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
389 (couple_pattern->ptr[i +1] > couple_pattern->ptr[i])) {
390 /* if i is an F row we first calculate alpha and beta: */
391
392 sum_all_neg=0; /* sum of all negative values in row i of A */
393 sum_all_pos=0; /* sum of all positive values in row i of A */
394 sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
395 sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
396 A_ii=0;
397
398 /* first check the mainBlock */
399 range = A->mainBlock->pattern->ptr[i + 1];
400 for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
401 j=A->mainBlock->pattern->index[iPtr];
402 A_ij=A->mainBlock->val[iPtr];
403 if(j==i) {
404 A_ii=A_ij;
405 } else {
406
407 if(A_ij< 0) {
408 sum_all_neg+=A_ij;
409 } else {
410 sum_all_pos+=A_ij;
411 }
412
413 if (counter_C[j]>=0) {
414 /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
415 start_p=&(main_pattern->index[main_pattern->ptr[i]]);
416 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
417 main_pattern->ptr[i + 1] - main_pattern->ptr[i],
418 sizeof(index_t),
419 paso::comparIndex);
420 if (! (where_p == NULL) ) { /* yes i strongly connected with j */
421 offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p);
422 main_block->val[offset]=A_ij; /* will be modified later */
423 if (A_ij< 0) {
424 sum_strong_neg+=A_ij;
425 } else {
426 sum_strong_pos+=A_ij;
427 }
428 }
429 }
430
431 }
432 }
433
434 /* now we deal with the col_coupleBlock */
435 range = A->col_coupleBlock->pattern->ptr[i + 1];
436 for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
437 j=A->col_coupleBlock->pattern->index[iPtr];
438 A_ij=A->col_coupleBlock->val[iPtr];
439 if(A_ij< 0) {
440 sum_all_neg+=A_ij;
441 } else {
442 sum_all_pos+=A_ij;
443 }
444
445 if (counter_C[j+my_n]>=0) {
446 /* is i strongly connect with j? We serach for counter_C[j] in P[i,:] */
447 start_p=&(couple_pattern->index[couple_pattern->ptr[i]]);
448 where_p=(index_t*)bsearch(&(counter_C[j+my_n]), start_p,
449 couple_pattern->ptr[i + 1] - couple_pattern->ptr[i],
450 sizeof(index_t),
451 paso::comparIndex);
452 if (! (where_p == NULL) ) { /* yes i strongly connect with j */
453 offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p);
454 couple_block->val[offset]=A_ij; /* will be modified later */
455 if (A_ij< 0) {
456 sum_strong_neg+=A_ij;
457 } else {
458 sum_strong_pos+=A_ij;
459 }
460 }
461 }
462 }
463
464 if(sum_strong_neg<0) {
465 alpha= sum_all_neg/sum_strong_neg;
466 } else {
467 alpha=0;
468 }
469 if(sum_strong_pos>0) {
470 beta= sum_all_pos/sum_strong_pos;
471 } else {
472 beta=0;
473 A_ii+=sum_all_pos;
474 }
475 if ( A_ii > 0.) {
476 rtmp=(-1.)/A_ii;
477 alpha*=rtmp;
478 beta*=rtmp;
479 }
480
481 range = main_pattern->ptr[i + 1];
482 for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
483 A_ij=main_block->val[iPtr];
484 if (A_ij > 0 ) {
485 main_block->val[iPtr] = A_ij * beta;
486 } else {
487 main_block->val[iPtr] = A_ij * alpha;
488 }
489 }
490
491 range = couple_pattern->ptr[i + 1];
492 for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
493 A_ij=couple_block->val[iPtr];
494 if (A_ij > 0 ) {
495 couple_block->val[iPtr] = A_ij * beta;
496 } else {
497 couple_block->val[iPtr] = A_ij * alpha;
498 }
499 }
500 }
501 }
502 }
503
504 void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SystemMatrix* P,
505 Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
506 const index_t* S, const index_t *counter_C)
507 {
508 paso::SparseMatrix_ptr main_block(P->mainBlock);
509 paso::SparseMatrix_ptr couple_block(P->col_coupleBlock);
510 paso::Pattern_ptr main_pattern(main_block->pattern);
511 paso::Pattern_ptr couple_pattern(couple_block->pattern);
512 const dim_t row_block_size=A->row_block_size;
513 const dim_t A_block = A->block_size;
514 const dim_t my_n=A->mainBlock->numRows;
515 index_t range;
516
517 dim_t i;
518 double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
519 register double A_ij, rtmp;
520 register index_t iPtr, j, offset, ib;
521 index_t *where_p, *start_p;
522
523 #pragma omp parallel private(i,offset,ib,sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos,A_ii,range,iPtr,j,A_ij,start_p,where_p,alpha,beta,rtmp)
524 {
525 sum_all_neg=new double[row_block_size]; /* sum of all negative values in row i of A */
526 sum_all_pos=new double[row_block_size]; /* sum of all positive values in row i of A */
527 sum_strong_neg=new double[row_block_size]; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
528 sum_strong_pos=new double[row_block_size]; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
529 alpha=new double[row_block_size];
530 beta=new double[row_block_size];
531 A_ii=new double[row_block_size];
532
533 #pragma omp for schedule(static)
534 for (i=0;i<my_n;++i) {
535 if (counter_C[i]>=0) { /* i is a C row */
536 offset = main_pattern->ptr[i];
537 for (ib =0; ib<row_block_size; ++ib)
538 main_block->val[row_block_size*offset+ib]=1.;
539 } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
540 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
541 /* if i is an F row we first calculate alpha and beta: */
542 for (ib =0; ib<row_block_size; ++ib) {
543 sum_all_neg[ib]=0;
544 sum_all_pos[ib]=0;
545 sum_strong_neg[ib]=0;
546 sum_strong_pos[ib]=0;
547 A_ii[ib]=0;
548 }
549
550 /* first check the mainBlock */
551 range = A->mainBlock->pattern->ptr[i + 1];
552 for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
553 j=A->mainBlock->pattern->index[iPtr];
554 if(j==i) {
555 for (ib =0; ib<row_block_size; ++ib)
556 A_ii[ib]=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
557 } else {
558 for (ib =0; ib<row_block_size; ++ib) {
559 A_ij=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
560 if(A_ij< 0) {
561 sum_all_neg[ib]+=A_ij;
562 } else {
563 sum_all_pos[ib]+=A_ij;
564 }
565 }
566
567 if (counter_C[j]>=0) {
568 /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
569 start_p=&(main_pattern->index[main_pattern->ptr[i]]);
570 where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
571 main_pattern->ptr[i + 1]-main_pattern->ptr[i],
572 sizeof(index_t),
573 paso::comparIndex);
574 if (! (where_p == NULL) ) { /* yes i strongly connected with j */
575 offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p);
576 for (ib =0; ib<row_block_size; ++ib) {
577 A_ij=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
578 main_block->val[row_block_size*offset+ib]=A_ij; /* will be modified later */
579 if (A_ij< 0) {
580 sum_strong_neg[ib]+=A_ij;
581 } else {
582 sum_strong_pos[ib]+=A_ij;
583 }
584 }
585 }
586 }
587
588 }
589 }
590
591 /* now we deal with the col_coupleBlock */
592 range = A->col_coupleBlock->pattern->ptr[i + 1];
593 for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
594 j=A->col_coupleBlock->pattern->index[iPtr];
595 for (ib =0; ib<row_block_size; ++ib) {
596 A_ij=A->col_coupleBlock->val[A_block*iPtr+ib+row_block_size*ib];
597 if(A_ij< 0) {
598 sum_all_neg[ib]+=A_ij;
599 } else {
600 sum_all_pos[ib]+=A_ij;
601 }
602 }
603
604 if (counter_C[j+my_n]>=0) {
605 /* is i strongly connect with j? We serach for counter_C[j] in P[i,:] */
606 start_p=&(couple_pattern->index[couple_pattern->ptr[i]]);
607 where_p=(index_t*)bsearch(&(counter_C[j+my_n]), start_p,
608 couple_pattern->ptr[i + 1]-couple_pattern->ptr[i],
609 sizeof(index_t),
610 paso::comparIndex);
611 if (! (where_p == NULL) ) { /* yes i strongly connect with j */
612 offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p);
613 for (ib =0; ib<row_block_size; ++ib) {
614 A_ij=A->col_coupleBlock->val[A_block*iPtr+ib+row_block_size*ib];
615 couple_block->val[row_block_size*offset+ib]=A_ij; /* will be modified later */
616
617 if (A_ij< 0) {
618 sum_strong_neg[ib]+=A_ij;
619 } else {
620 sum_strong_pos[ib]+=A_ij;
621 }
622 }
623 }
624 }
625 }
626
627 for (ib =0; ib<row_block_size; ++ib) {
628 if(sum_strong_neg[ib]<0) {
629 alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];
630 } else {
631 alpha[ib]=0;
632 }
633 if(sum_strong_pos[ib]>0) {
634 beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib];
635 } else {
636 beta[ib]=0;
637 A_ii[ib]+=sum_all_pos[ib];
638 }
639 if ( A_ii[ib] > 0.) {
640 rtmp=(-1./A_ii[ib]);
641 alpha[ib]*=rtmp;
642 beta[ib]*=rtmp;
643 }
644 }
645
646 range = main_pattern->ptr[i + 1];
647 for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
648 for (ib =0; ib<row_block_size; ++ib) {
649 A_ij=main_block->val[row_block_size*iPtr+ib];
650 if (A_ij > 0 ) {
651 main_block->val[row_block_size*iPtr+ib] = A_ij * beta[ib];
652 } else {
653 main_block->val[row_block_size*iPtr+ib] = A_ij * alpha[ib];
654 }
655 }
656 }
657
658 range = couple_pattern->ptr[i + 1];
659 for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
660 for (ib =0; ib<row_block_size; ++ib) {
661 A_ij=couple_block->val[row_block_size*iPtr+ib];
662 if (A_ij > 0 ) {
663 couple_block->val[row_block_size*iPtr+ib] = A_ij * beta[ib];
664 } else {
665 couple_block->val[row_block_size*iPtr+ib] = A_ij * alpha[ib];
666 }
667 }
668 }
669
670
671 }
672 }/* end i loop */
673 delete[] sum_all_neg;
674 delete[] sum_all_pos;
675 delete[] sum_strong_neg;
676 delete[] sum_strong_pos;
677 delete[] alpha;
678 delete[] beta;
679 delete[] A_ii;
680 } /* end parallel region */
681 }
682
683 /*
684 Classic Prolongation:
685 -------------------
686
687 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
688 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])
689 where the summation over l is considering columns which are strongly connected
690 to i (l in S[i]) and not in C (counter_C[l]<0) and
691
692 B[i,k]=sum_{m in S_i and in C} A+[k,m]
693 a[i]=A[i,i]+sum{l not strongly connected to i} A[i,l]
694
695 A+[i,k]=A[i,k] if sign(A[i,k])==sign(A[i,i]) or 0 otherwise
696
697
698 */
699 void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SystemMatrix* P,
700 Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
701 const index_t* S, const index_t *counter_C)
702 {
703 paso::SparseMatrix_ptr main_block(P->mainBlock);
704 paso::SparseMatrix_ptr couple_block(P->col_coupleBlock);
705 paso::Pattern_ptr main_pattern(main_block->pattern);
706 paso::Pattern_ptr couple_pattern(couple_block->pattern);
707 const dim_t my_n=A->mainBlock->numRows;
708 index_t range, range_j;
709
710 dim_t i, q;
711 double *D_s=NULL;
712 index_t *D_s_offset=NULL, iPtr, iPtr_j;
713 const index_t *ptr_main_A = A->mainBlock->borrowMainDiagonalPointer();
714 index_t *start_p_main_i, *start_p_couple_i;
715 dim_t degree_p_main_i, degree_p_couple_i;
716 dim_t len, ll, main_len, len_D_s;
717
718 len = A->col_coupleBlock->numCols;
719 len = MAX(A->remote_coupleBlock->numCols, len);
720 ll = len + my_n;
721 main_len = main_pattern->len;
722
723 #pragma omp parallel private(D_s,D_s_offset,i,start_p_main_i,start_p_couple_i,degree_p_main_i,degree_p_couple_i,range,iPtr,q,range_j,iPtr_j,len_D_s)
724 {
725 D_s=new double[ll];
726 D_s_offset=new index_t[ll];
727
728 #pragma omp for schedule(static)
729 for (i=0;i<my_n;++i) {
730 if (counter_C[i]>=0) {
731 main_block->val[main_pattern->ptr[i]]=1.; /* i is a C row */
732 } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
733 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
734 /* this loop sums up the weak connections in a and creates
735 a list of the strong connected columns which are not in
736 C (=no interpolation nodes) */
737 const index_t *start_s = &(S[offset_S[i]]);
738 const double A_ii = A->mainBlock->val[ptr_main_A[i]];
739 double a=A_ii;
740
741 start_p_main_i = &(main_pattern->index[main_pattern->ptr[i]]);
742 start_p_couple_i = &(couple_pattern->index[couple_pattern->ptr[i]]);
743 degree_p_main_i = main_pattern->ptr[i+1] - main_pattern->ptr[i];
744 degree_p_couple_i = couple_pattern->ptr[i+1] - couple_pattern->ptr[i];
745
746 /* first, check the mainBlock */
747 range = A->mainBlock->pattern->ptr[i + 1];
748 for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
749 const index_t j=A->mainBlock->pattern->index[iPtr];
750 const double A_ij=A->mainBlock->val[iPtr];
751 if ( (i!=j) && (degree_S[j]>0) ) {
752 /* is (i,j) a strong connection ?*/
753 const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), paso::comparIndex);
754 if (where_s == NULL) { /* weak connections are accumulated */
755 a+=A_ij;
756 } else { /* yes i strongly connected with j */
757 if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */
758 const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
759 if (where_p == NULL) {
760 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_setClassicProlongation: interpolation point is missing.");
761 } else {
762 const index_t offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p_main_i);
763 main_block->val[offset]+=A_ij;
764 }
765 } else { /* j is not an interpolation point */
766 /* find all interpolation points m of k */
767 double s=0.;
768 len_D_s=0;
769
770 /* first, the mainBlock part */
771 range_j = A->mainBlock->pattern->ptr[j + 1];
772 for (iPtr_j=A->mainBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
773 const double A_jm=A->mainBlock->val[iPtr_j];
774 const index_t m=A->mainBlock->pattern->index[iPtr_j];
775 /* is m an interpolation point ? */
776 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
777 if (! (where_p_m==NULL)) {
778 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
779 if (! SAMESIGN(A_ii,A_jm)) {
780 D_s[len_D_s]=A_jm;
781 } else {
782 D_s[len_D_s]=0.;
783 }
784 D_s_offset[len_D_s]=offset_m;
785 len_D_s++;
786 }
787 }
788
789 /* then the coupleBlock part */
790 if (degree_p_couple_i) {
791 range_j = A->col_coupleBlock->pattern->ptr[j + 1];
792 for (iPtr_j=A->col_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
793 const double A_jm=A->col_coupleBlock->val[iPtr_j];
794 const index_t m=A->col_coupleBlock->pattern->index[iPtr_j];
795 /* is m an interpolation point ? */
796 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), paso::comparIndex);
797 if (! (where_p_m==NULL)) {
798 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
799 if (! SAMESIGN(A_ii,A_jm)) {
800 D_s[len_D_s]=A_jm;
801 } else {
802 D_s[len_D_s]=0.;
803 }
804 D_s_offset[len_D_s]=offset_m + main_len;
805 len_D_s++;
806 }
807 }
808 }
809
810 for (q=0;q<len_D_s;++q) s+=D_s[q];
811 if (ABS(s)>0) {
812 s=A_ij/s;
813 for (q=0;q<len_D_s;++q) {
814 if (D_s_offset[q] < main_len)
815 main_block->val[D_s_offset[q]]+=s*D_s[q];
816 else
817 couple_block->val[D_s_offset[q]-main_len]+=s*D_s[q];
818 }
819 } else {
820 a+=A_ij;
821 }
822 }
823 }
824 }
825 }
826
827 if (A->mpi_info->size > 1) {
828 /* now, deal with the coupleBlock */
829 range = A->col_coupleBlock->pattern->ptr[i + 1];
830 for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
831 const index_t j=A->col_coupleBlock->pattern->index[iPtr];
832 const double A_ij=A->col_coupleBlock->val[iPtr];
833 if ( (i!=j) && (degree_S[j]>0) ) {
834 /* is (i,j) a strong connection ?*/
835 index_t t=j+my_n;
836 const index_t *where_s=(index_t*)bsearch(&t, start_s,degree_S[i],sizeof(index_t), paso::comparIndex);
837 if (where_s == NULL) { /* weak connections are accumulated */
838 a+=A_ij;
839 } else { /* yes i strongly connect with j */
840 if (counter_C[t]>=0) { /* j is an interpolation point : add A_ij into P */
841 const index_t *where_p=(index_t*)bsearch(&counter_C[t], start_p_couple_i,degree_p_couple_i, sizeof(index_t), paso::comparIndex);
842 if (where_p == NULL) {
843 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");
844 } else {
845 const index_t offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p_couple_i);
846 couple_block->val[offset]+=A_ij;
847 }
848 } else { /* j is not an interpolation point */
849 /* find all interpolation points m of k */
850 double s=0.;
851 len_D_s=0;
852
853 /* first, the row_coupleBlock part */
854 range_j = A->row_coupleBlock->pattern->ptr[j + 1];
855 for (iPtr_j=A->row_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
856 const double A_jm=A->row_coupleBlock->val[iPtr_j];
857 const index_t m=A->row_coupleBlock->pattern->index[iPtr_j];
858 /* is m an interpolation point ? */
859 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
860 if (! (where_p_m==NULL)) {
861 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
862 if (! SAMESIGN(A_ii,A_jm)) {
863 D_s[len_D_s]=A_jm;
864 } else {
865 D_s[len_D_s]=0.;
866 }
867 D_s_offset[len_D_s]=offset_m;
868 len_D_s++;
869 }
870 }
871
872 /* then the remote_coupleBlock part */
873 range_j = A->remote_coupleBlock->pattern->ptr[j + 1];
874 for (iPtr_j=A->remote_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
875 const double A_jm=A->remote_coupleBlock->val[iPtr_j];
876 const index_t m=A->remote_coupleBlock->pattern->index[iPtr_j];
877 /* is m an interpolation point ? */
878 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i, degree_p_couple_i, sizeof(index_t), paso::comparIndex);
879 if (! (where_p_m==NULL)) {
880 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
881 if (! SAMESIGN(A_ii,A_jm)) {
882 D_s[len_D_s]=A_jm;
883 } else {
884 D_s[len_D_s]=0.;
885 }
886 D_s_offset[len_D_s]=offset_m + main_len;
887 len_D_s++;
888 }
889 }
890
891 for (q=0;q<len_D_s;++q) s+=D_s[q];
892 if (ABS(s)>0) {
893 s=A_ij/s;
894 for (q=0;q<len_D_s;++q) {
895 if (D_s_offset[q] < main_len)
896 main_block->val[D_s_offset[q]]+=s*D_s[q];
897 else
898 couple_block->val[D_s_offset[q]-main_len]+=s*D_s[q];
899 }
900 } else {
901 a+=A_ij;
902 }
903 }
904 }
905 }
906 }
907 }
908
909 /* i has been processed, now we need to do some rescaling */
910 if (ABS(a)>0.) {
911 a=-1./a;
912 range = main_pattern->ptr[i + 1];
913 for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
914 main_block->val[iPtr]*=a;
915 }
916
917 range = couple_pattern->ptr[i + 1];
918 for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
919 couple_block->val[iPtr]*=a;
920 }
921 }
922 }
923 } /* end of row i loop */
924 delete[] D_s;
925 delete[] D_s_offset;
926 } /* end of parallel region */
927 }
928
929 void Paso_Preconditioner_AMG_setClassicProlongation_Block(
930 Paso_SystemMatrix* P, Paso_SystemMatrix* A, const index_t* offset_S,
931 const dim_t* degree_S, const index_t* S, const index_t *counter_C)
932 {
933 paso::SparseMatrix_ptr main_block(P->mainBlock);
934 paso::SparseMatrix_ptr couple_block(P->col_coupleBlock);
935 paso::Pattern_ptr main_pattern(main_block->pattern);
936 paso::Pattern_ptr couple_pattern(couple_block->pattern);
937 const dim_t row_block=A->row_block_size;
938 const dim_t my_n=A->mainBlock->numRows;
939 const dim_t A_block = A->block_size;
940 index_t range, range_j;
941
942 dim_t i, q, ib;
943 double *D_s=NULL;
944 index_t *D_s_offset=NULL, iPtr, iPtr_j;
945 index_t *start_p_main_i, *start_p_couple_i;
946 dim_t degree_p_main_i, degree_p_couple_i;
947 dim_t len, ll, main_len, len_D_s;
948 const index_t *ptr_main_A = A->mainBlock->borrowMainDiagonalPointer();
949
950 len = A->col_coupleBlock->numCols;
951 len = MAX(A->remote_coupleBlock->numCols, len);
952 ll = len + my_n;
953 main_len = main_pattern->len;
954 #pragma omp parallel private(D_s,D_s_offset,i,ib,start_p_main_i,start_p_couple_i,degree_p_main_i,degree_p_couple_i,range,iPtr,q,range_j,iPtr_j,len_D_s)
955 {
956 double *a=new double[row_block];
957 D_s=new double[row_block*ll];
958 D_s_offset=new index_t[row_block*ll];
959
960 #pragma omp for private(i) schedule(static)
961 for (i=0;i<my_n;++i) {
962 if (counter_C[i]>=0) {
963 const index_t offset = main_pattern->ptr[i];
964 for (ib =0; ib<row_block; ++ib) main_block->val[row_block*offset+ib]=1.; /* i is a C row */
965 } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
966 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
967 /* this loop sums up the weak connections in a and creates
968 a list of the strong connected columns which are not in
969 C (=no interpolation nodes) */
970 const index_t *start_s = &(S[offset_S[i]]);
971 const double *A_ii = &(A->mainBlock->val[ptr_main_A[i]*A_block]);
972 start_p_main_i = &(main_pattern->index[main_pattern->ptr[i]]);
973 start_p_couple_i = &(couple_pattern->index[couple_pattern->ptr[i]]);
974 degree_p_main_i = main_pattern->ptr[i+1] - main_pattern->ptr[i];
975 degree_p_couple_i = couple_pattern->ptr[i+1] - couple_pattern->ptr[i];
976 for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];
977
978 /* first, check the mainBlock */
979 range = A->mainBlock->pattern->ptr[i + 1];
980 for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
981 const index_t j=A->mainBlock->pattern->index[iPtr];
982 const double* A_ij=&(A->mainBlock->val[iPtr*A_block]);
983
984 if ( (i!=j) && (degree_S[j]>0) ) {
985 /* is (i,j) a strong connection ?*/
986 const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), paso::comparIndex);
987 if (where_s == NULL) { /* weak connections are accumulated */
988 for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
989 } else { /* yes i strongly connected with j */
990 if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */
991 const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
992 if (where_p == NULL) {
993 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
994 } else {
995 const index_t offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p_main_i);
996 for (ib=0; ib<row_block; ib++) main_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
997 }
998 } else { /* j is not an interpolation point */
999 /* find all interpolation points m of k */
1000 len_D_s=0;
1001
1002 /* first, the mainBlock part */
1003 range_j = A->mainBlock->pattern->ptr[j + 1];
1004 for (iPtr_j=A->mainBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1005 const double* A_jm=&(A->mainBlock->val[iPtr_j*A_block]);
1006 const index_t m=A->mainBlock->pattern->index[iPtr_j];
1007 /* is m an interpolation point ? */
1008 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
1009 if (! (where_p_m==NULL)) {
1010 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1011 for (ib=0; ib<row_block; ib++) {
1012 if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1013 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1014 } else {
1015 D_s[len_D_s*row_block+ib]=0.;
1016 }
1017 }
1018 D_s_offset[len_D_s] = offset_m;
1019 len_D_s++;
1020 }
1021 }
1022
1023 /* then the coupleBlock part */
1024 range_j = A->col_coupleBlock->pattern->ptr[j+1];
1025 for (iPtr_j=A->col_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1026 const double* A_jm=&(A->col_coupleBlock->val[iPtr_j*A_block]);
1027 const index_t m=A->col_coupleBlock->pattern->index[iPtr_j];
1028 /* is m an interpolation point ? */
1029 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), paso::comparIndex);
1030 if (! (where_p_m==NULL)) {
1031 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1032 for (ib=0; ib<row_block; ib++) {
1033 if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1034 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1035 } else {
1036 D_s[len_D_s*row_block+ib]=0.;
1037 }
1038 }
1039 D_s_offset[len_D_s]=offset_m + main_len;
1040 len_D_s++;
1041 }
1042 }
1043
1044 for (ib=0; ib<row_block; ib++) {
1045 double s=0;
1046 for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1047
1048 if (ABS(s)>0) {
1049 s=A_ij[(row_block+1)*ib]/s;
1050 for (q=0; q<len_D_s; q++) {
1051 if (D_s_offset[q] < main_len)
1052 main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1053 else{
1054 couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1055 }
1056 }
1057 } else {
1058 a[ib]+=A_ij[(row_block+1)*ib];
1059 }
1060 }
1061 }
1062 }
1063 }
1064 }
1065
1066 if (A->mpi_info->size > 1) {
1067 /* now, deal with the coupleBlock */
1068 range = A->col_coupleBlock->pattern->ptr[i + 1];
1069 for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
1070 const index_t j=A->col_coupleBlock->pattern->index[iPtr];
1071 const double* A_ij=&(A->col_coupleBlock->val[iPtr*A_block]);
1072
1073 if ( (i!=j) && (degree_S[j]>0) ) {
1074 /* is (i,j) a strong connection ?*/
1075 index_t t=j+my_n;
1076 const index_t *where_s=(index_t*)bsearch(&t, start_s,degree_S[i],sizeof(index_t), paso::comparIndex);
1077 if (where_s == NULL) { /* weak connections are accumulated */
1078 for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
1079 } else { /* yes i strongly connected with j */
1080 if (counter_C[t]>=0) { /* j is an interpolation point : add A_ij into P */
1081 const index_t *where_p=(index_t*)bsearch(&counter_C[t], start_p_couple_i,degree_p_couple_i, sizeof(index_t), paso::comparIndex);
1082 if (where_p == NULL) {
1083 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
1084
1085 } else {
1086 const index_t offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p_couple_i);
1087 for (ib=0; ib<row_block; ib++) couple_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
1088 }
1089 } else { /* j is not an interpolation point */
1090 /* find all interpolation points m of k */
1091 len_D_s=0;
1092
1093 /* first, the row_coupleBlock part */
1094 range_j = A->row_coupleBlock->pattern->ptr[j + 1];
1095 for (iPtr_j=A->row_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1096 /* is m an interpolation point ? */
1097 index_t m, *where_p_m;
1098 double *A_jm;
1099 A_jm=&(A->row_coupleBlock->val[iPtr_j*A_block]);
1100 m=A->row_coupleBlock->pattern->index[iPtr_j];
1101
1102 where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), paso::comparIndex);
1103 if (! (where_p_m==NULL)) {
1104 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1105 for (ib=0; ib<row_block; ib++) {
1106 if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1107 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1108 } else {
1109 D_s[len_D_s*row_block+ib]=0.;
1110 }
1111 }
1112 D_s_offset[len_D_s]=offset_m;
1113 len_D_s++;
1114 }
1115 }
1116
1117 /* then the remote_coupleBlock part */
1118 if (degree_p_couple_i) {
1119 range_j = A->remote_coupleBlock->pattern->ptr[j + 1];
1120 for (iPtr_j=A->remote_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1121 const double* A_jm=&(A->remote_coupleBlock->val[iPtr_j*A_block]);
1122 const index_t m=A->remote_coupleBlock->pattern->index[iPtr_j];
1123 /* is m an interpolation point ? */
1124 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), paso::comparIndex);
1125 if (! (where_p_m==NULL)) {
1126 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1127 for (ib=0; ib<row_block; ib++) {
1128 if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1129 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1130 } else {
1131 D_s[len_D_s*row_block+ib]=0.;
1132 }
1133 }
1134 D_s_offset[len_D_s]=offset_m + main_len;
1135 len_D_s++;
1136 }
1137 }
1138 }
1139
1140 for (ib=0; ib<row_block; ib++) {
1141 double s=0;
1142 for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1143
1144 if (ABS(s)>0) {
1145 s=A_ij[(row_block+1)*ib]/s;
1146 for (q=0;q<len_D_s;++q) {
1147 if (D_s_offset[q] < main_len)
1148 main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1149 else
1150 couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1151 }
1152 } else {
1153 a[ib]+=A_ij[(row_block+1)*ib];
1154 }
1155 }
1156 }
1157 }
1158 }
1159 }
1160 }
1161
1162 /* i has been processed, now we need to do some rescaling */
1163 for (ib=0; ib<row_block; ib++) {
1164 register double a2=a[ib];
1165 if (ABS(a2)>0.) {
1166 a2=-1./a2;
1167 range = main_pattern->ptr[i + 1];
1168 for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
1169 main_block->val[iPtr*row_block+ib]*=a2;
1170 }
1171
1172 range = couple_pattern->ptr[i + 1];
1173 for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
1174 couple_block->val[iPtr*row_block+ib]*=a2;
1175 }
1176 }
1177 }
1178 }
1179 } /* end of row i loop */
1180 delete[] D_s;
1181 delete[] D_s_offset;
1182 delete[] a;
1183 } /* end of parallel region */
1184 }
1185

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26