/[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 4873 - (show annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 1 month ago) by caltinay
File size: 55375 byte(s)
whitespace only changes.

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 #include <cstring> // memset
33
34 namespace paso {
35
36 /****************************************************************************
37
38 Methods necessary for AMG preconditioner
39
40 construct n x n_C the prolongation matrix P from A_p.
41
42 The columns in A_p to be considered are marked by counter_C[n] where
43 an unknown i to be considered in P is marked by 0<= counter_C[i] < n_C
44 and counter_C[i] gives the new column number in P.
45 S defines the strong connections.
46
47 The pattern of P is formed as follows:
48
49 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise.
50 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.
51
52 Two settings for P are implemented (see below)
53
54 */
55
56 SystemMatrix_ptr Preconditioner_AMG_getProlongation(
57 SystemMatrix_ptr A_p, const index_t* offset_S, const dim_t* degree_S,
58 const index_t* S, const dim_t n_C, index_t* counter_C,
59 const index_t interpolation_method)
60 {
61 Esys_MPIInfo *mpi_info=Esys_MPIInfo_getReference(A_p->mpi_info);
62 Distribution_ptr input_dist, output_dist;
63 SharedComponents_ptr send, recv;
64 Connector_ptr col_connector;
65 Pattern_ptr main_pattern, couple_pattern;
66 const dim_t row_block_size=A_p->row_block_size;
67 const dim_t col_block_size=A_p->col_block_size;
68 const dim_t my_n=A_p->mainBlock->numCols;
69 const dim_t overlap_n=A_p->col_coupleBlock->numCols;
70 const dim_t num_threads=omp_get_max_threads();
71 index_t size=mpi_info->size, *dist=NULL;
72 index_t *main_p=NULL, *couple_p=NULL, *main_idx=NULL, *couple_idx=NULL;
73 index_t *shared=NULL, *offsetInShared=NULL;
74 index_t *recv_shared=NULL, *send_shared=NULL;
75 index_t sum, i, j, k, l, p, q, iptr;
76 index_t my_n_C, global_label, num_neighbors;
77 #ifdef ESYS_MPI
78 index_t rank=mpi_info->rank;
79 #endif
80 Esys_MPI_rank *neighbor=NULL;
81 #ifdef ESYS_MPI
82 MPI_Request* mpi_requests=NULL;
83 MPI_Status* mpi_stati=NULL;
84 #else
85 int *mpi_requests=NULL, *mpi_stati=NULL;
86 #endif
87
88 /* number of C points in current distribution */
89 my_n_C = 0;
90 sum=0;
91 if (num_threads>1) {
92 #pragma omp parallel private(i,sum)
93 {
94 sum=0;
95 #pragma omp for schedule(static)
96 for (i=0;i<my_n;++i) {
97 if (counter_C[i] != -1) {
98 sum++;
99 }
100 }
101 #pragma omp critical
102 {
103 my_n_C += sum;
104 }
105 }
106 } else { /* num_threads=1 */
107 for (i=0;i<my_n;++i) {
108 if (counter_C[i] != -1) {
109 my_n_C++;
110 }
111 }
112 }
113
114 /* create row distribution (output_distribution) and col distribution
115 (input_distribution) */
116 /* ??? should I alloc an new Esys_MPIInfo object or reuse the one in
117 system matrix A. for now, I'm reuse A->mpi_info ??? */
118 dist = A_p->pattern->output_distribution->first_component;
119 output_dist.reset(new Distribution(mpi_info, dist, 1, 0));
120 dist = new index_t[size+1]; /* now prepare for col distribution */
121 #ifdef ESYS_MPI
122 MPI_Allgather(&my_n_C, 1, MPI_INT, dist, 1, MPI_INT, mpi_info->comm);
123 #endif
124 global_label=0;
125 for (i=0; i<size; i++) {
126 k = dist[i];
127 dist[i] = global_label;
128 global_label += k;
129 }
130 dist[size] = global_label;
131
132 input_dist.reset(new Distribution(mpi_info, dist, 1, 0));
133 delete[] dist;
134
135 /* create pattern for mainBlock and coupleBlock */
136 main_p = new index_t[my_n+1];
137 couple_p = new index_t[my_n+1];
138 /* count the number of entries per row in the Prolongation matrix :*/
139 #pragma omp parallel for private(i,l,k,iptr,j,p) schedule(static)
140 for (i=0; i<my_n; i++) {
141 l = 0;
142 if (counter_C[i]>=0) {
143 k = 1; /* i is a C unknown */
144 } else {
145 k = 0;
146 iptr = offset_S[i];
147 for (p=0; p<degree_S[i]; p++) {
148 j = S[iptr+p]; /* this is a strong connection */
149 if (counter_C[j]>=0) { /* and is in C */
150 if (j <my_n) k++;
151 else {
152 l++;
153 }
154 }
155 }
156 }
157 main_p[i] = k;
158 couple_p[i] = l;
159 }
160
161 /* number of unknowns in the col-coupleBlock of the interpolation matrix */
162 sum = 0;
163 for (i=0;i<overlap_n;++i) {
164 if (counter_C[i+my_n] > -1) {
165 counter_C[i+my_n] -= my_n_C;
166 sum++;
167 }
168 }
169
170 /* allocate and create index vector for prolongation: */
171 p = util::cumsum(my_n, main_p);
172 main_p[my_n] = p;
173 main_idx = new index_t[p];
174 p = util::cumsum(my_n, couple_p);
175 couple_p[my_n] = p;
176 couple_idx = new index_t[p];
177 #pragma omp parallel for private(i,k,l,iptr,j,p) schedule(static)
178 for (i=0; i<my_n; i++) {
179 if (counter_C[i]>=0) {
180 main_idx[main_p[i]]=counter_C[i]; /* i is a C unknown */
181 } else {
182 k = 0;
183 l = 0;
184 iptr = offset_S[i];
185 for (p=0; p<degree_S[i]; p++) {
186 j = S[iptr+p]; /* this is a strong connection */
187 if (counter_C[j] >=0) { /* and is in C */
188 if (j < my_n) {
189 main_idx[main_p[i]+k] = counter_C[j];
190 k++;
191 } else {
192 couple_idx[couple_p[i]+l] = counter_C[j];
193 l++;
194 }
195 }
196 }
197 }
198 }
199
200 if (Esys_noError()) {
201 main_pattern.reset(new Pattern(MATRIX_FORMAT_DEFAULT, my_n,
202 my_n_C, main_p, main_idx));
203 couple_pattern.reset(new Pattern(MATRIX_FORMAT_DEFAULT, my_n,
204 sum, couple_p, couple_idx));
205 } else {
206 delete[] main_p;
207 delete[] main_idx;
208 delete[] couple_p;
209 delete[] couple_idx;
210 }
211
212 /* prepare the receiver for the col_connector.
213 Note that the allocation for "shared" assumes the send and receive buffer
214 of the interpolation matrix P is no larger than that of matrix A_p. */
215 neighbor = new Esys_MPI_rank[size];
216 offsetInShared = new index_t[size+1];
217 recv = A_p->col_coupler->connector->recv;
218 send = A_p->col_coupler->connector->send;
219 i = recv->numSharedComponents;
220 recv_shared = new index_t[i];
221 memset(recv_shared, 0, sizeof(index_t)*i);
222 k = send->numSharedComponents;
223 send_shared = new index_t[k];
224 if (k > i) i = k;
225 shared = new index_t[i];
226
227 #ifdef ESYS_MPI
228 mpi_requests=new MPI_Request[size*2];
229 mpi_stati=new MPI_Status[size*2];
230 #else
231 mpi_requests=new int[size*2];
232 mpi_stati=new int[size*2];
233 #endif
234
235 for (p=0; p<send->numNeighbors; p++) {
236 i = send->offsetInShared[p];
237 #ifdef ESYS_MPI
238 MPI_Irecv (&(send_shared[i]), send->offsetInShared[p+1]-i, MPI_INT,
239 send->neighbor[p], mpi_info->msg_tag_counter+send->neighbor[p],
240 mpi_info->comm, &mpi_requests[p]);
241 #endif
242 }
243
244 num_neighbors = 0;
245 q = 0;
246 p = recv->numNeighbors;
247 offsetInShared[0]=0;
248 for (i=0; i<p; i++) {
249 l = 0;
250 k = recv->offsetInShared[i+1];
251 for (j=recv->offsetInShared[i]; j<k; j++) {
252 if (counter_C[recv->shared[j]] > -1) {
253 shared[q] = my_n_C + q;
254 recv_shared[recv->shared[j]-my_n] = 1;
255 q++;
256 l = 1;
257 }
258 }
259 if (l == 1) {
260 iptr = recv->neighbor[i];
261 neighbor[num_neighbors] = iptr;
262 num_neighbors++;
263 offsetInShared[num_neighbors] = q;
264 }
265 #ifdef ESYS_MPI
266 MPI_Issend(&(recv_shared[recv->offsetInShared[i]]),
267 k-recv->offsetInShared[i], MPI_INT, recv->neighbor[i],
268 mpi_info->msg_tag_counter+rank, mpi_info->comm,
269 &mpi_requests[i+send->numNeighbors]);
270 #endif
271 }
272 recv.reset(new SharedComponents(my_n_C, num_neighbors, neighbor,
273 shared, offsetInShared, 1, 0, mpi_info));
274
275 /* now we can build the sender */
276 #ifdef ESYS_MPI
277 MPI_Waitall(recv->numNeighbors+send->numNeighbors, mpi_requests, mpi_stati);
278 #endif
279 ESYS_MPI_INC_COUNTER(*mpi_info, size)
280 delete[] mpi_requests;
281 delete[] mpi_stati;
282
283 num_neighbors = 0;
284 q = 0;
285 p = send->numNeighbors;
286 offsetInShared[0]=0;
287 for (i=0; i<p; i++) {
288 l = 0;
289 k = send->offsetInShared[i+1];
290 for (j=send->offsetInShared[i]; j<k; j++) {
291 if (send_shared[j] == 1) {
292 shared[q] = counter_C[send->shared[j]];
293 q++;
294 l = 1;
295 }
296 }
297 if (l == 1) {
298 iptr = send->neighbor[i];
299 neighbor[num_neighbors] = iptr;
300 num_neighbors++;
301 offsetInShared[num_neighbors] = q;
302 }
303 }
304
305 send.reset(new SharedComponents(my_n_C, num_neighbors, neighbor,
306 shared, offsetInShared, 1, 0, mpi_info));
307 col_connector.reset(new Connector(send, recv));
308 delete[] recv_shared;
309 delete[] send_shared;
310 delete[] neighbor;
311 delete[] offsetInShared;
312 delete[] shared;
313
314 /* now we need to create the System Matrix
315 TO BE FIXED: at this stage, we only construct col_couple_pattern
316 and col_connector for interpolation matrix P. To be completed,
317 row_couple_pattern and row_connector need to be constructed as well */
318 SystemMatrix_ptr out;
319 SystemMatrixPattern_ptr pattern;
320 if (Esys_noError()) {
321 pattern.reset(new SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
322 output_dist, input_dist, main_pattern, couple_pattern,
323 couple_pattern, col_connector, col_connector));
324 out.reset(new SystemMatrix(MATRIX_FORMAT_DIAGONAL_BLOCK, pattern,
325 row_block_size, col_block_size, false));
326 }
327
328 /* now fill in the matrix */
329 if (Esys_noError()) {
330 if ((interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
331 || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {
332 if (row_block_size == 1) {
333 Preconditioner_AMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);
334 } else {
335 Preconditioner_AMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
336 }
337 } else {
338 if (row_block_size == 1) {
339 Preconditioner_AMG_setDirectProlongation(out, A_p, offset_S, degree_S, S, counter_C);
340 } else {
341 Preconditioner_AMG_setDirectProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
342 }
343 }
344 }
345
346 if (!Esys_noError()) {
347 out.reset();
348 }
349 return out;
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 Preconditioner_AMG_setDirectProlongation(SystemMatrix_ptr P,
368 SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S,
369 const index_t* S, const index_t *counter_C)
370 {
371 SparseMatrix_ptr main_block(P->mainBlock);
372 SparseMatrix_ptr couple_block(P->col_coupleBlock);
373 Pattern_ptr main_pattern(main_block->pattern);
374 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 util::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 util::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 Preconditioner_AMG_setDirectProlongation_Block(SystemMatrix_ptr P,
505 SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S,
506 const index_t* S, const index_t *counter_C)
507 {
508 SparseMatrix_ptr main_block(P->mainBlock);
509 SparseMatrix_ptr couple_block(P->col_coupleBlock);
510 Pattern_ptr main_pattern(main_block->pattern);
511 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 util::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 util::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 Preconditioner_AMG_setClassicProlongation(SystemMatrix_ptr P,
700 SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S,
701 const index_t* S, const index_t *counter_C)
702 {
703 SparseMatrix_ptr main_block(P->mainBlock);
704 SparseMatrix_ptr couple_block(P->col_coupleBlock);
705 Pattern_ptr main_pattern(main_block->pattern);
706 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 = std::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), util::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), util::comparIndex);
759 if (where_p == NULL) {
760 Esys_setError(SYSTEM_ERROR, "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), util::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 (!util::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), util::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 (!util::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 (std::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), util::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), util::comparIndex);
842 if (where_p == NULL) {
843 Esys_setError(SYSTEM_ERROR, "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), util::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 (!util::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), util::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 (!util::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 (std::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 (std::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 Preconditioner_AMG_setClassicProlongation_Block(
930 SystemMatrix_ptr P, SystemMatrix_ptr A, const index_t* offset_S,
931 const dim_t* degree_S, const index_t* S, const index_t *counter_C)
932 {
933 SparseMatrix_ptr main_block(P->mainBlock);
934 SparseMatrix_ptr couple_block(P->col_coupleBlock);
935 Pattern_ptr main_pattern(main_block->pattern);
936 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 dim_t i, q, ib;
942 double *D_s=NULL;
943 index_t *D_s_offset=NULL, iPtr, iPtr_j;
944 index_t *start_p_main_i, *start_p_couple_i;
945 dim_t degree_p_main_i, degree_p_couple_i;
946 dim_t len, ll, main_len, len_D_s;
947 const index_t *ptr_main_A = A->mainBlock->borrowMainDiagonalPointer();
948
949 len = A->col_coupleBlock->numCols;
950 len = std::max(A->remote_coupleBlock->numCols, len);
951 ll = len + my_n;
952 main_len = main_pattern->len;
953 #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)
954 {
955 double *a=new double[row_block];
956 D_s=new double[row_block*ll];
957 D_s_offset=new index_t[row_block*ll];
958
959 #pragma omp for private(i) schedule(static)
960 for (i=0;i<my_n;++i) {
961 if (counter_C[i]>=0) {
962 const index_t offset = main_pattern->ptr[i];
963 for (ib =0; ib<row_block; ++ib) main_block->val[row_block*offset+ib]=1.; /* i is a C row */
964 } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
965 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
966 /* this loop sums up the weak connections in a and creates
967 a list of the strong connected columns which are not in
968 C (=no interpolation nodes) */
969 const index_t *start_s = &(S[offset_S[i]]);
970 const double *A_ii = &(A->mainBlock->val[ptr_main_A[i]*A_block]);
971 start_p_main_i = &(main_pattern->index[main_pattern->ptr[i]]);
972 start_p_couple_i = &(couple_pattern->index[couple_pattern->ptr[i]]);
973 degree_p_main_i = main_pattern->ptr[i+1] - main_pattern->ptr[i];
974 degree_p_couple_i = couple_pattern->ptr[i+1] - couple_pattern->ptr[i];
975 for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];
976
977 /* first, check the mainBlock */
978 range = A->mainBlock->pattern->ptr[i + 1];
979 for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
980 const index_t j=A->mainBlock->pattern->index[iPtr];
981 const double* A_ij=&(A->mainBlock->val[iPtr*A_block]);
982
983 if ( (i!=j) && (degree_S[j]>0) ) {
984 /* is (i,j) a strong connection ?*/
985 const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), util::comparIndex);
986 if (where_s == NULL) { /* weak connections are accumulated */
987 for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
988 } else { /* yes i strongly connected with j */
989 if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */
990 const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p_main_i,degree_p_main_i, sizeof(index_t), util::comparIndex);
991 if (where_p == NULL) {
992 Esys_setError(SYSTEM_ERROR, "Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
993 } else {
994 const index_t offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p_main_i);
995 for (ib=0; ib<row_block; ib++) main_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
996 }
997 } else { /* j is not an interpolation point */
998 /* find all interpolation points m of k */
999 len_D_s=0;
1000
1001 /* first, the mainBlock part */
1002 range_j = A->mainBlock->pattern->ptr[j + 1];
1003 for (iPtr_j=A->mainBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1004 const double* A_jm=&(A->mainBlock->val[iPtr_j*A_block]);
1005 const index_t m=A->mainBlock->pattern->index[iPtr_j];
1006 /* is m an interpolation point ? */
1007 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), util::comparIndex);
1008 if (! (where_p_m==NULL)) {
1009 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1010 for (ib=0; ib<row_block; ib++) {
1011 if (!util::samesign(A_ii[(row_block+1)*ib], A_jm[(row_block+1)*ib])) {
1012 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1013 } else {
1014 D_s[len_D_s*row_block+ib]=0.;
1015 }
1016 }
1017 D_s_offset[len_D_s] = offset_m;
1018 len_D_s++;
1019 }
1020 }
1021
1022 /* then the coupleBlock part */
1023 range_j = A->col_coupleBlock->pattern->ptr[j+1];
1024 for (iPtr_j=A->col_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1025 const double* A_jm=&(A->col_coupleBlock->val[iPtr_j*A_block]);
1026 const index_t m=A->col_coupleBlock->pattern->index[iPtr_j];
1027 /* is m an interpolation point ? */
1028 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), util::comparIndex);
1029 if (! (where_p_m==NULL)) {
1030 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1031 for (ib=0; ib<row_block; ib++) {
1032 if (!util::samesign(A_ii[(row_block+1)*ib], A_jm[(row_block+1)*ib]) ) {
1033 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1034 } else {
1035 D_s[len_D_s*row_block+ib]=0.;
1036 }
1037 }
1038 D_s_offset[len_D_s]=offset_m + main_len;
1039 len_D_s++;
1040 }
1041 }
1042
1043 for (ib=0; ib<row_block; ib++) {
1044 double s=0;
1045 for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1046
1047 if (std::abs(s)>0) {
1048 s=A_ij[(row_block+1)*ib]/s;
1049 for (q=0; q<len_D_s; q++) {
1050 if (D_s_offset[q] < main_len)
1051 main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1052 else{
1053 couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1054 }
1055 }
1056 } else {
1057 a[ib]+=A_ij[(row_block+1)*ib];
1058 }
1059 }
1060 }
1061 }
1062 }
1063 }
1064
1065 if (A->mpi_info->size > 1) {
1066 /* now, deal with the coupleBlock */
1067 range = A->col_coupleBlock->pattern->ptr[i + 1];
1068 for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
1069 const index_t j=A->col_coupleBlock->pattern->index[iPtr];
1070 const double* A_ij=&(A->col_coupleBlock->val[iPtr*A_block]);
1071
1072 if ( (i!=j) && (degree_S[j]>0) ) {
1073 /* is (i,j) a strong connection ?*/
1074 index_t t=j+my_n;
1075 const index_t *where_s=(index_t*)bsearch(&t, start_s,degree_S[i],sizeof(index_t), util::comparIndex);
1076 if (where_s == NULL) { /* weak connections are accumulated */
1077 for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
1078 } else { /* yes i strongly connected with j */
1079 if (counter_C[t]>=0) { /* j is an interpolation point : add A_ij into P */
1080 const index_t *where_p=(index_t*)bsearch(&counter_C[t], start_p_couple_i,degree_p_couple_i, sizeof(index_t), util::comparIndex);
1081 if (where_p == NULL) {
1082 Esys_setError(SYSTEM_ERROR, "Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
1083
1084 } else {
1085 const index_t offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p_couple_i);
1086 for (ib=0; ib<row_block; ib++) couple_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
1087 }
1088 } else { /* j is not an interpolation point */
1089 /* find all interpolation points m of k */
1090 len_D_s=0;
1091
1092 /* first, the row_coupleBlock part */
1093 range_j = A->row_coupleBlock->pattern->ptr[j + 1];
1094 for (iPtr_j=A->row_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1095 /* is m an interpolation point ? */
1096 index_t m, *where_p_m;
1097 double *A_jm;
1098 A_jm=&(A->row_coupleBlock->val[iPtr_j*A_block]);
1099 m=A->row_coupleBlock->pattern->index[iPtr_j];
1100
1101 where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), util::comparIndex);
1102 if (! (where_p_m==NULL)) {
1103 const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1104 for (ib=0; ib<row_block; ib++) {
1105 if (!util::samesign(A_ii[(row_block+1)*ib], A_jm[(row_block+1)*ib])) {
1106 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1107 } else {
1108 D_s[len_D_s*row_block+ib]=0.;
1109 }
1110 }
1111 D_s_offset[len_D_s]=offset_m;
1112 len_D_s++;
1113 }
1114 }
1115
1116 /* then the remote_coupleBlock part */
1117 if (degree_p_couple_i) {
1118 range_j = A->remote_coupleBlock->pattern->ptr[j + 1];
1119 for (iPtr_j=A->remote_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1120 const double* A_jm=&(A->remote_coupleBlock->val[iPtr_j*A_block]);
1121 const index_t m=A->remote_coupleBlock->pattern->index[iPtr_j];
1122 /* is m an interpolation point ? */
1123 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), util::comparIndex);
1124 if (! (where_p_m==NULL)) {
1125 const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1126 for (ib=0; ib<row_block; ib++) {
1127 if (!util::samesign(A_ii[(row_block+1)*ib], A_jm[(row_block+1)*ib]) ) {
1128 D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1129 } else {
1130 D_s[len_D_s*row_block+ib]=0.;
1131 }
1132 }
1133 D_s_offset[len_D_s]=offset_m + main_len;
1134 len_D_s++;
1135 }
1136 }
1137 }
1138
1139 for (ib=0; ib<row_block; ib++) {
1140 double s=0;
1141 for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1142
1143 if (std::abs(s)>0) {
1144 s=A_ij[(row_block+1)*ib]/s;
1145 for (q=0;q<len_D_s;++q) {
1146 if (D_s_offset[q] < main_len)
1147 main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1148 else
1149 couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1150 }
1151 } else {
1152 a[ib]+=A_ij[(row_block+1)*ib];
1153 }
1154 }
1155 }
1156 }
1157 }
1158 }
1159 }
1160
1161 /* i has been processed, now we need to do some rescaling */
1162 for (ib=0; ib<row_block; ib++) {
1163 register double a2=a[ib];
1164 if (std::abs(a2)>0.) {
1165 a2=-1./a2;
1166 range = main_pattern->ptr[i + 1];
1167 for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
1168 main_block->val[iPtr*row_block+ib]*=a2;
1169 }
1170
1171 range = couple_pattern->ptr[i + 1];
1172 for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
1173 couple_block->val[iPtr*row_block+ib]*=a2;
1174 }
1175 }
1176 }
1177 }
1178 } /* end of row i loop */
1179 delete[] D_s;
1180 delete[] D_s_offset;
1181 delete[] a;
1182 } /* end of parallel region */
1183 }
1184
1185 } // namespace paso
1186

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