/[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 4803 - (show annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 9 months ago) by caltinay
File size: 51878 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

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

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