source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_lr_interp.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 163.5 KB
RevLine 
[2aa6644]1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14#include "headers.h"
15#include "aux_interp.h"
16
17#define MAX_C_CONNECTIONS 100
18#define HAVE_COMMON_C 1
19
20/*---------------------------------------------------------------------------
21 * hypre_BoomerAMGBuildStdInterp
22 * Comment: The interpolatory weighting can be changed with the sep_weight
23 * variable. This can enable not separating negative and positive
24 * off diagonals in the weight formula.
25 *--------------------------------------------------------------------------*/
26
27int
28hypre_BoomerAMGBuildStdInterp(hypre_ParCSRMatrix *A, int *CF_marker,
29 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
30 int num_functions, int *dof_func, int debug_flag,
31 double trunc_factor, int max_elmts,
32 int sep_weight, int *col_offd_S_to_A,
33 hypre_ParCSRMatrix **P_ptr)
34{
35 /* Communication Variables */
36 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
37 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
38 int my_id, num_procs;
39
40 /* Variables to store input variables */
41 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
42 double *A_diag_data = hypre_CSRMatrixData(A_diag);
43 int *A_diag_i = hypre_CSRMatrixI(A_diag);
44 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
45
46 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
47 double *A_offd_data = hypre_CSRMatrixData(A_offd);
48 int *A_offd_i = hypre_CSRMatrixI(A_offd);
49 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
50
51 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
52 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
53 int n_fine = hypre_CSRMatrixNumRows(A_diag);
54 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
55 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
56 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt)local_numrows;
57 HYPRE_BigInt total_global_cpts, my_first_cpt;
58
59 /* Variables to store strong connection matrix info */
60 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
61 int *S_diag_i = hypre_CSRMatrixI(S_diag);
62 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
63
64 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
65 int *S_offd_i = hypre_CSRMatrixI(S_offd);
66 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
67
68 /* Interpolation matrix P */
69 hypre_ParCSRMatrix *P;
70 hypre_CSRMatrix *P_diag;
71 hypre_CSRMatrix *P_offd;
72
73 double *P_diag_data = NULL;
74 int *P_diag_i, *P_diag_j = NULL;
75 double *P_offd_data = NULL;
76 int *P_offd_i, *P_offd_j = NULL;
77
78 HYPRE_BigInt *col_map_offd_P;
79 int *tmp_map_offd;
80 int P_diag_size;
81 int P_offd_size;
82 int *P_marker = NULL;
83 int *P_marker_offd = NULL;
84 int *CF_marker_offd = NULL;
85 int *tmp_CF_marker_offd = NULL;
86 int *dof_func_offd = NULL;
87
88 /* Full row information for columns of A that are off diag*/
89 hypre_BigCSRMatrix *A_ext;
90 double *A_ext_data;
91 int *A_ext_i;
92 HYPRE_BigInt *big_A_ext_j;
93 int *A_ext_j;
94
95 int *fine_to_coarse = NULL;
96 HYPRE_BigInt *fine_to_coarse_offd = NULL;
97 HYPRE_BigInt *found;
98
99 int num_cols_P_offd;
100 int newoff, loc_col;
101 int A_ext_rows, full_off_procNodes;
102
103 hypre_BigCSRMatrix *Sop;
104 int *Sop_i;
105 HYPRE_BigInt *big_Sop_j;
106 int *Sop_j;
107
108 int Soprows;
109
110 /* Variables to keep count of interpolatory points */
111 int jj_counter, jj_counter_offd;
112 int jj_begin_row, jj_end_row;
113 int jj_begin_row_offd = 0;
114 int jj_end_row_offd = 0;
115 int coarse_counter, coarse_counter_offd;
116 int *ihat = NULL;
117 int *ihat_offd = NULL;
118 int *ipnt = NULL;
119 int *ipnt_offd = NULL;
120 int strong_f_marker = -2;
121
122 /* Interpolation weight variables */
123 double *ahat = NULL;
124 double *ahat_offd = NULL;
125 double sum_pos, sum_pos_C, sum_neg, sum_neg_C, sum, sum_C;
126 double diagonal, distribute;
127 double alfa, beta;
128
129 /* Loop variables */
130 int index;
131 int start_indexing = 0;
132 int i, i1, j1, jj, kk, k1;
133 int cnt_c, cnt_f, cnt_c_offd, cnt_f_offd, indx;
134
135 /* Definitions */
136 double zero = 0.0;
137 double one = 1.0;
138 double wall_time;
139 double wall_1 = 0;
140 double wall_2 = 0;
141 double wall_3 = 0;
142
143
144 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
145
146 if (debug_flag== 4) wall_time = time_getWallclockSeconds();
147
148 /* BEGIN */
149 MPI_Comm_size(comm, &num_procs);
150 MPI_Comm_rank(comm,&my_id);
151
152#ifdef HYPRE_NO_GLOBAL_PARTITION
153 my_first_cpt = num_cpts_global[0];
154 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
155 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
156#else
157 my_first_cpt = num_cpts_global[my_id];
158 total_global_cpts = num_cpts_global[num_procs];
159#endif
160
161 if (!comm_pkg)
162 {
163#ifdef HYPRE_NO_GLOBAL_PARTITION
164 hypre_NewCommPkgCreate(A);
165#else
166 hypre_MatvecCommPkgCreate(A);
167#endif
168 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
169 }
170
171 /* Set up off processor information (specifically for neighbors of
172 * neighbors */
173 newoff = 0;
174 full_off_procNodes = 0;
175 if (num_procs > 1)
176 {
177 /*----------------------------------------------------------------------
178 * Get the off processors rows for A and S, associated with columns in
179 * A_offd and S_offd.
180 *---------------------------------------------------------------------*/
181 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
182 A_ext_i = hypre_BigCSRMatrixI(A_ext);
183 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
184 A_ext_data = hypre_BigCSRMatrixData(A_ext);
185 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
186
187 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
188 Sop_i = hypre_BigCSRMatrixI(Sop);
189 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
190 Soprows = hypre_BigCSRMatrixNumRows(Sop);
191
192 /* Find nodes that are neighbors of neighbors, not found in offd */
193 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
194 &A_ext_j, Soprows, col_map_offd, col_1,
195 col_n, Sop_i, big_Sop_j, &Sop_j, CF_marker,
196 comm_pkg);
197
198 hypre_BigCSRMatrixJ(A_ext) = NULL;
199 hypre_BigCSRMatrixJ(Sop) = NULL;
200
201 if(newoff >= 0)
202 full_off_procNodes = newoff + num_cols_A_offd;
203 else
204 return(1);
205
206 /* Possibly add new points and new processors to the comm_pkg, all
207 * processors need new_comm_pkg */
208
209 /* AHB - create a new comm package just for extended info -
210 this will work better with the assumed partition*/
211 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
212 &extend_comm_pkg);
213
214 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
215
216 if (num_functions > 1 && full_off_procNodes > 0)
217 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
218
219 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
220 full_off_procNodes, CF_marker_offd);
221
222 if(num_functions > 1)
223 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
224 full_off_procNodes, dof_func_offd);
225 }
226
227
228 /*-----------------------------------------------------------------------
229 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
230 *-----------------------------------------------------------------------*/
231
232 /*-----------------------------------------------------------------------
233 * Intialize counters and allocate mapping vector.
234 *-----------------------------------------------------------------------*/
235 P_diag_i = hypre_CTAlloc(int, n_fine+1);
236 P_offd_i = hypre_CTAlloc(int, n_fine+1);
237
238 if (n_fine)
239 {
240 fine_to_coarse = hypre_CTAlloc(int, n_fine);
241 P_marker = hypre_CTAlloc(int, n_fine);
242 }
243
244 if (full_off_procNodes)
245 {
246 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
247 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
248 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
249 }
250
251 for (i=0; i < full_off_procNodes; i++)
252 {
253 P_marker_offd[i] = -1;
254 tmp_CF_marker_offd[i] = -1;
255 fine_to_coarse_offd[i] = -1;
256 }
257
258 for (i=0; i < n_fine; i++)
259 {
260 P_marker[i] = -1;
261 fine_to_coarse[i] = -1;
262 }
263
264 jj_counter = start_indexing;
265 jj_counter_offd = start_indexing;
266 coarse_counter = 0;
267 coarse_counter_offd = 0;
268
269 /*-----------------------------------------------------------------------
270 * Loop over fine grid.
271 *-----------------------------------------------------------------------*/
272 for (i = 0; i < n_fine; i++)
273 {
274 P_diag_i[i] = jj_counter;
275 if (num_procs > 1)
276 P_offd_i[i] = jj_counter_offd;
277
278 if (CF_marker[i] >= 0)
279 {
280 jj_counter++;
281 fine_to_coarse[i] = coarse_counter;
282 coarse_counter++;
283 }
284
285 /*--------------------------------------------------------------------
286 * If i is an F-point, interpolation is from the C-points that
287 * strongly influence i, or C-points that stronly influence F-points
288 * that strongly influence i.
289 *--------------------------------------------------------------------*/
290 else if (CF_marker[i] != -3)
291 {
292 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
293 {
294 i1 = S_diag_j[jj];
295 if (CF_marker[i1] >= 0)
296 { /* i1 is a C point */
297 if (P_marker[i1] < P_diag_i[i])
298 {
299 P_marker[i1] = jj_counter;
300 jj_counter++;
301 }
302 }
303 else if (CF_marker[i1] != -3)
304 { /* i1 is a F point, loop through it's strong neighbors */
305 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
306 {
307 k1 = S_diag_j[kk];
308 if (CF_marker[k1] >= 0)
309 {
310 if(P_marker[k1] < P_diag_i[i])
311 {
312 P_marker[k1] = jj_counter;
313 jj_counter++;
314 }
315 }
316 }
317 if(num_procs > 1)
318 {
319 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
320 {
321 if(col_offd_S_to_A)
322 k1 = col_offd_S_to_A[S_offd_j[kk]];
323 else
324 k1 = S_offd_j[kk];
325 if (CF_marker_offd[k1] > 0)
326 {
327 if(P_marker_offd[k1] < P_offd_i[i])
328 {
329 tmp_CF_marker_offd[k1] = 1;
330 P_marker_offd[k1] = jj_counter_offd;
331 jj_counter_offd++;
332 }
333 }
334 }
335 }
336 }
337 }
338 /* Look at off diag strong connections of i */
339 if (num_procs > 1)
340 {
341 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
342 {
343 i1 = S_offd_j[jj];
344 if(col_offd_S_to_A)
345 i1 = col_offd_S_to_A[i1];
346 if (CF_marker_offd[i1] > 0)
347 {
348 if(P_marker_offd[i1] < P_offd_i[i])
349 {
350 tmp_CF_marker_offd[i1] = 1;
351 P_marker_offd[i1] = jj_counter_offd;
352 jj_counter_offd++;
353 }
354 }
355 else if (CF_marker_offd[i1] != -3)
356 { /* F point; look at neighbors of i1. Sop contains global col
357 * numbers and entries that could be in S_diag or S_offd or
358 * neither. */
359 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
360 {
361 loc_col = Sop_j[kk];
362 if(loc_col > -1)
363 { /* In S_diag */
364 if(CF_marker[loc_col] > 0)
365 {
366 if(P_marker[loc_col] < P_diag_i[i])
367 {
368 P_marker[loc_col] = jj_counter;
369 jj_counter++;
370 }
371 }
372 }
373 else
374 {
375 loc_col = - loc_col - 1;
376 if(CF_marker_offd[loc_col] > 0)
377 {
378 if(P_marker_offd[loc_col] < P_offd_i[i])
379 {
380 P_marker_offd[loc_col] = jj_counter_offd;
381 tmp_CF_marker_offd[loc_col] = 1;
382 jj_counter_offd++;
383 }
384 }
385 }
386 }
387 }
388 }
389 }
390 }
391 }
392
393 if (debug_flag==4)
394 {
395 wall_time = time_getWallclockSeconds() - wall_time;
396 printf("Proc = %d determine structure %f\n",
397 my_id, wall_time);
398 fflush(NULL);
399 }
400 /*-----------------------------------------------------------------------
401 * Allocate arrays.
402 *-----------------------------------------------------------------------*/
403
404
405 P_diag_size = jj_counter;
406 P_offd_size = jj_counter_offd;
407
408 if (P_diag_size)
409 {
410 P_diag_j = hypre_CTAlloc(int, P_diag_size);
411 P_diag_data = hypre_CTAlloc(double, P_diag_size);
412 }
413 if (P_offd_size)
414 {
415 P_offd_j = hypre_CTAlloc(int, P_offd_size);
416 P_offd_data = hypre_CTAlloc(double, P_offd_size);
417 }
418
419 P_diag_i[n_fine] = jj_counter;
420 P_offd_i[n_fine] = jj_counter_offd;
421
422 jj_counter = start_indexing;
423 jj_counter_offd = start_indexing;
424
425 /* Fine to coarse mapping */
426 if(num_procs > 1)
427 {
428 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
429 full_off_procNodes, my_first_cpt,
430 fine_to_coarse_offd);
431 }
432
433 /* Initialize ahat, which is a modification to a, used in the standard
434 * interpolation routine. */
435 if (n_fine)
436 {
437 ahat = hypre_CTAlloc(double, n_fine);
438 ihat = hypre_CTAlloc(int, n_fine);
439 ipnt = hypre_CTAlloc(int, n_fine);
440 }
441 if (full_off_procNodes)
442 {
443 ahat_offd = hypre_CTAlloc(double, full_off_procNodes);
444 ihat_offd = hypre_CTAlloc(int, full_off_procNodes);
445 ipnt_offd = hypre_CTAlloc(int, full_off_procNodes);
446 }
447
448 for (i = 0; i < n_fine; i++)
449 {
450 P_marker[i] = -1;
451 ahat[i] = 0;
452 ihat[i] = -1;
453 }
454 for (i = 0; i < full_off_procNodes; i++)
455 {
456 P_marker_offd[i] = -1;
457 ahat_offd[i] = 0;
458 ihat_offd[i] = -1;
459 }
460
461 /*-----------------------------------------------------------------------
462 * Loop over fine grid points.
463 *-----------------------------------------------------------------------*/
464 for (i = 0; i < n_fine; i++)
465 {
466 jj_begin_row = jj_counter;
467 if(num_procs > 1)
468 jj_begin_row_offd = jj_counter_offd;
469
470 /*--------------------------------------------------------------------
471 * If i is a c-point, interpolation is the identity.
472 *--------------------------------------------------------------------*/
473
474 if (CF_marker[i] >= 0)
475 {
476 P_diag_j[jj_counter] = fine_to_coarse[i];
477 P_diag_data[jj_counter] = one;
478 jj_counter++;
479 }
480
481 /*--------------------------------------------------------------------
482 * If i is an F-point, build interpolation.
483 *--------------------------------------------------------------------*/
484
485 else if (CF_marker[i] != -3)
486 {
487 if (debug_flag==4) wall_time = time_getWallclockSeconds();
488 strong_f_marker--;
489 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
490 {
491 i1 = S_diag_j[jj];
492
493 /*--------------------------------------------------------------
494 * If neighbor i1 is a C-point, set column number in P_diag_j
495 * and initialize interpolation weight to zero.
496 *--------------------------------------------------------------*/
497
498 if (CF_marker[i1] >= 0)
499 {
500 if (P_marker[i1] < jj_begin_row)
501 {
502 P_marker[i1] = jj_counter;
503 P_diag_j[jj_counter] = i1;
504 P_diag_data[jj_counter] = zero;
505 jj_counter++;
506 }
507 }
508 else if (CF_marker[i1] != -3)
509 {
510 P_marker[i1] = strong_f_marker;
511 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
512 {
513 k1 = S_diag_j[kk];
514 if (CF_marker[k1] >= 0)
515 {
516 if(P_marker[k1] < jj_begin_row)
517 {
518 P_marker[k1] = jj_counter;
519 P_diag_j[jj_counter] = k1;
520 P_diag_data[jj_counter] = zero;
521 jj_counter++;
522 }
523 }
524 }
525 if(num_procs > 1)
526 {
527 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
528 {
529 if(col_offd_S_to_A)
530 k1 = col_offd_S_to_A[S_offd_j[kk]];
531 else
532 k1 = S_offd_j[kk];
533 if(CF_marker_offd[k1] >= 0)
534 {
535 if(P_marker_offd[k1] < jj_begin_row_offd)
536 {
537 P_marker_offd[k1] = jj_counter_offd;
538 P_offd_j[jj_counter_offd] = k1;
539 P_offd_data[jj_counter_offd] = zero;
540 jj_counter_offd++;
541 }
542 }
543 }
544 }
545 }
546 }
547
548 if ( num_procs > 1)
549 {
550 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
551 {
552 i1 = S_offd_j[jj];
553 if(col_offd_S_to_A)
554 i1 = col_offd_S_to_A[i1];
555 if ( CF_marker_offd[i1] >= 0)
556 {
557 if(P_marker_offd[i1] < jj_begin_row_offd)
558 {
559 P_marker_offd[i1] = jj_counter_offd;
560 P_offd_j[jj_counter_offd]=i1;
561 P_offd_data[jj_counter_offd] = zero;
562 jj_counter_offd++;
563 }
564 }
565 else if (CF_marker_offd[i1] != -3)
566 {
567 P_marker_offd[i1] = strong_f_marker;
568 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
569 {
570 loc_col = Sop_j[kk];
571 if(loc_col > -1)
572 {
573 if(CF_marker[loc_col] > 0)
574 {
575 if(P_marker[loc_col] < jj_begin_row)
576 {
577 P_marker[loc_col] = jj_counter;
578 P_diag_j[jj_counter] = loc_col;
579 P_diag_data[jj_counter] = zero;
580 jj_counter++;
581 }
582 }
583 }
584 else
585 {
586 loc_col = -loc_col - 1;
587 if(CF_marker_offd[loc_col] > 0)
588 {
589 if(P_marker_offd[loc_col] < jj_begin_row_offd)
590 {
591 P_marker_offd[loc_col] = jj_counter_offd;
592 P_offd_j[jj_counter_offd]=loc_col;
593 P_offd_data[jj_counter_offd] = zero;
594 jj_counter_offd++;
595 }
596 }
597 }
598 }
599 }
600 }
601 }
602
603 jj_end_row = jj_counter;
604 jj_end_row_offd = jj_counter_offd;
605
606 if (debug_flag==4)
607 {
608 wall_time = time_getWallclockSeconds() - wall_time;
609 wall_1 += wall_time;
610 fflush(NULL);
611 }
612 if (debug_flag==4) wall_time = time_getWallclockSeconds();
613 cnt_c = 0;
614 cnt_f = jj_end_row-jj_begin_row;
615 cnt_c_offd = 0;
616 cnt_f_offd = jj_end_row_offd-jj_begin_row_offd;
617 ihat[i] = cnt_f;
618 ipnt[cnt_f] = i;
619 ahat[cnt_f++] = A_diag_data[A_diag_i[i]];
620 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
621 { /* i1 is direct neighbor */
622 i1 = A_diag_j[jj];
623 if (P_marker[i1] != strong_f_marker)
624 {
625 indx = ihat[i1];
626 if (indx > -1)
627 ahat[indx] += A_diag_data[jj];
628 else if (P_marker[i1] >= jj_begin_row)
629 {
630 ihat[i1] = cnt_c;
631 ipnt[cnt_c] = i1;
632 ahat[cnt_c++] += A_diag_data[jj];
633 }
634 else if (CF_marker[i1] != -3)
635 {
636 ihat[i1] = cnt_f;
637 ipnt[cnt_f] = i1;
638 ahat[cnt_f++] += A_diag_data[jj];
639 }
640 }
641 else
642 {
643 if(num_functions == 1 || dof_func[i] == dof_func[i1])
644 {
645 distribute = A_diag_data[jj]/A_diag_data[A_diag_i[i1]];
646 for (kk = A_diag_i[i1]+1; kk < A_diag_i[i1+1]; kk++)
647 {
648 k1 = A_diag_j[kk];
649 indx = ihat[k1];
650 if (indx > -1)
651 ahat[indx] -= A_diag_data[kk]*distribute;
652 else if (P_marker[k1] >= jj_begin_row)
653 {
654 ihat[k1] = cnt_c;
655 ipnt[cnt_c] = k1;
656 ahat[cnt_c++] -= A_diag_data[kk]*distribute;
657 }
658 else
659 {
660 ihat[k1] = cnt_f;
661 ipnt[cnt_f] = k1;
662 ahat[cnt_f++] -= A_diag_data[kk]*distribute;
663 }
664 }
665 if(num_procs > 1)
666 {
667 for (kk = A_offd_i[i1]; kk < A_offd_i[i1+1]; kk++)
668 {
669 k1 = A_offd_j[kk];
670 indx = ihat_offd[k1];
671 if(num_functions == 1 || dof_func[i1] == dof_func_offd[k1])
672 {
673 if (indx > -1)
674 ahat_offd[indx] -= A_offd_data[kk]*distribute;
675 else if (P_marker_offd[k1] >= jj_begin_row_offd)
676 {
677 ihat_offd[k1] = cnt_c_offd;
678 ipnt_offd[cnt_c_offd] = k1;
679 ahat_offd[cnt_c_offd++] -= A_offd_data[kk]*distribute;
680 }
681 else
682 {
683 ihat_offd[k1] = cnt_f_offd;
684 ipnt_offd[cnt_f_offd] = k1;
685 ahat_offd[cnt_f_offd++] -= A_offd_data[kk]*distribute;
686 }
687 }
688 }
689 }
690 }
691 }
692 }
693 if(num_procs > 1)
694 {
695 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
696 {
697 i1 = A_offd_j[jj];
698 if(P_marker_offd[i1] != strong_f_marker)
699 {
700 indx = ihat_offd[i1];
701 if (indx > -1)
702 ahat_offd[indx] += A_offd_data[jj];
703 else if (P_marker_offd[i1] >= jj_begin_row_offd)
704 {
705 ihat_offd[i1] = cnt_c_offd;
706 ipnt_offd[cnt_c_offd] = i1;
707 ahat_offd[cnt_c_offd++] += A_offd_data[jj];
708 }
709 else if (CF_marker_offd[i1] != -3)
710 {
711 ihat_offd[i1] = cnt_f_offd;
712 ipnt_offd[cnt_f_offd] = i1;
713 ahat_offd[cnt_f_offd++] += A_offd_data[jj];
714 }
715 }
716 else
717 {
718 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
719 {
720 distribute = A_offd_data[jj]/A_ext_data[A_ext_i[i1]];
721 for (kk = A_ext_i[i1]+1; kk < A_ext_i[i1+1]; kk++)
722 {
723 loc_col = A_ext_j[kk];
724 if(loc_col > -1)
725 { /*diag*/
726 indx = ihat[loc_col];
727 if (indx > -1)
728 ahat[indx] -= A_ext_data[kk]*distribute;
729 else if (P_marker[loc_col] >= jj_begin_row)
730 {
731 ihat[loc_col] = cnt_c;
732 ipnt[cnt_c] = loc_col;
733 ahat[cnt_c++] -= A_ext_data[kk]*distribute;
734 }
735 else
736 {
737 ihat[loc_col] = cnt_f;
738 ipnt[cnt_f] = loc_col;
739 ahat[cnt_f++] -= A_ext_data[kk]*distribute;
740 }
741 }
742 else
743 {
744 loc_col = - loc_col - 1;
745 if(num_functions == 1 ||
746 dof_func_offd[loc_col] == dof_func_offd[i1])
747 {
748 indx = ihat_offd[loc_col];
749 if (indx > -1)
750 ahat_offd[indx] -= A_ext_data[kk]*distribute;
751 else if(P_marker_offd[loc_col] >= jj_begin_row_offd)
752 {
753 ihat_offd[loc_col] = cnt_c_offd;
754 ipnt_offd[cnt_c_offd] = loc_col;
755 ahat_offd[cnt_c_offd++] -= A_ext_data[kk]*distribute;
756 }
757 else
758 {
759 ihat_offd[loc_col] = cnt_f_offd;
760 ipnt_offd[cnt_f_offd] = loc_col;
761 ahat_offd[cnt_f_offd++] -= A_ext_data[kk]*distribute;
762 }
763 }
764 }
765 }
766 }
767 }
768 }
769 }
770 if (debug_flag==4)
771 {
772 wall_time = time_getWallclockSeconds() - wall_time;
773 wall_2 += wall_time;
774 fflush(NULL);
775 }
776
777 if (debug_flag==4) wall_time = time_getWallclockSeconds();
778 diagonal = ahat[cnt_c];
779 ahat[cnt_c] = 0;
780 sum_pos = 0;
781 sum_pos_C = 0;
782 sum_neg = 0;
783 sum_neg_C = 0;
784 sum = 0;
785 sum_C = 0;
786 if(sep_weight == 1)
787 {
788 for (jj=0; jj < cnt_c; jj++)
789 {
790 if (ahat[jj] > 0)
791 {
792 sum_pos_C += ahat[jj];
793 }
794 else
795 {
796 sum_neg_C += ahat[jj];
797 }
798 }
799 if(num_procs > 1)
800 {
801 for (jj=0; jj < cnt_c_offd; jj++)
802 {
803 if (ahat_offd[jj] > 0)
804 {
805 sum_pos_C += ahat_offd[jj];
806 }
807 else
808 {
809 sum_neg_C += ahat_offd[jj];
810 }
811 }
812 }
813 sum_pos = sum_pos_C;
814 sum_neg = sum_neg_C;
815 for (jj=cnt_c+1; jj < cnt_f; jj++)
816 {
817 if (ahat[jj] > 0)
818 {
819 sum_pos += ahat[jj];
820 }
821 else
822 {
823 sum_neg += ahat[jj];
824 }
825 ahat[jj] = 0;
826 }
827 if(num_procs > 1)
828 {
829 for (jj=cnt_c_offd; jj < cnt_f_offd; jj++)
830 {
831 if (ahat_offd[jj] > 0)
832 {
833 sum_pos += ahat_offd[jj];
834 }
835 else
836 {
837 sum_neg += ahat_offd[jj];
838 }
839 ahat_offd[jj] = 0;
840 }
841 }
842 if (sum_neg_C) alfa = sum_neg/sum_neg_C/diagonal;
843 if (sum_pos_C) beta = sum_pos/sum_pos_C/diagonal;
844
845 /*-----------------------------------------------------------------
846 * Set interpolation weight by dividing by the diagonal.
847 *-----------------------------------------------------------------*/
848
849 for (jj = jj_begin_row; jj < jj_end_row; jj++)
850 {
851 j1 = ihat[P_diag_j[jj]];
852 if (ahat[j1] > 0)
853 P_diag_data[jj] = -beta*ahat[j1];
854 else
855 P_diag_data[jj] = -alfa*ahat[j1];
856
857 P_diag_j[jj] = fine_to_coarse[P_diag_j[jj]];
858 ahat[j1] = 0;
859 }
860 for (jj=0; jj < cnt_f; jj++)
861 ihat[ipnt[jj]] = -1;
862 if(num_procs > 1)
863 {
864 for (jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
865 {
866 j1 = ihat_offd[P_offd_j[jj]];
867 if (ahat_offd[j1] > 0)
868 P_offd_data[jj] = -beta*ahat_offd[j1];
869 else
870 P_offd_data[jj] = -alfa*ahat_offd[j1];
871
872 ahat_offd[j1] = 0;
873 }
874 for (jj=0; jj < cnt_f_offd; jj++)
875 ihat_offd[ipnt_offd[jj]] = -1;
876 }
877 }
878 else
879 {
880 for (jj=0; jj < cnt_c; jj++)
881 {
882 sum_C += ahat[jj];
883 }
884 if(num_procs > 1)
885 {
886 for (jj=0; jj < cnt_c_offd; jj++)
887 {
888 sum_C += ahat_offd[jj];
889 }
890 }
891 sum = sum_C;
892 for (jj=cnt_c+1; jj < cnt_f; jj++)
893 {
894 sum += ahat[jj];
895 ahat[jj] = 0;
896 }
897 if(num_procs > 1)
898 {
899 for (jj=cnt_c_offd; jj < cnt_f_offd; jj++)
900 {
901 sum += ahat_offd[jj];
902 ahat_offd[jj] = 0;
903 }
904 }
905 if (sum_C) alfa = sum/sum_C/diagonal;
906
907 /*-----------------------------------------------------------------
908 * Set interpolation weight by dividing by the diagonal.
909 *-----------------------------------------------------------------*/
910
911 for (jj = jj_begin_row; jj < jj_end_row; jj++)
912 {
913 j1 = ihat[P_diag_j[jj]];
914 P_diag_data[jj] = -alfa*ahat[j1];
915 P_diag_j[jj] = fine_to_coarse[P_diag_j[jj]];
916 ahat[j1] = 0;
917 }
918 for (jj=0; jj < cnt_f; jj++)
919 ihat[ipnt[jj]] = -1;
920 if(num_procs > 1)
921 {
922 for (jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
923 {
924 j1 = ihat_offd[P_offd_j[jj]];
925 P_offd_data[jj] = -alfa*ahat_offd[j1];
926 ahat_offd[j1] = 0;
927 }
928 for (jj=0; jj < cnt_f_offd; jj++)
929 ihat_offd[ipnt_offd[jj]] = -1;
930 }
931 }
932 if (debug_flag==4)
933 {
934 wall_time = time_getWallclockSeconds() - wall_time;
935 wall_3 += wall_time;
936 fflush(NULL);
937 }
938 }
939 }
940
941 if (debug_flag==4)
942 {
943 printf("Proc = %d fill part 1 %f part 2 %f part 3 %f\n",
944 my_id, wall_1, wall_2, wall_3);
945 fflush(NULL);
946 }
947 P = hypre_ParCSRMatrixCreate(comm,
948 hypre_ParCSRMatrixGlobalNumRows(A),
949 total_global_cpts,
950 hypre_ParCSRMatrixColStarts(A),
951 num_cpts_global,
952 0,
953 P_diag_i[n_fine],
954 P_offd_i[n_fine]);
955
956 P_diag = hypre_ParCSRMatrixDiag(P);
957 hypre_CSRMatrixData(P_diag) = P_diag_data;
958 hypre_CSRMatrixI(P_diag) = P_diag_i;
959 hypre_CSRMatrixJ(P_diag) = P_diag_j;
960 P_offd = hypre_ParCSRMatrixOffd(P);
961 hypre_CSRMatrixData(P_offd) = P_offd_data;
962 hypre_CSRMatrixI(P_offd) = P_offd_i;
963 hypre_CSRMatrixJ(P_offd) = P_offd_j;
964 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
965
966 /* Compress P, removing coefficients smaller than trunc_factor * Max */
967 if (trunc_factor != 0.0 || max_elmts > 0)
968 {
969 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
970 P_diag_data = hypre_CSRMatrixData(P_diag);
971 P_diag_i = hypre_CSRMatrixI(P_diag);
972 P_diag_j = hypre_CSRMatrixJ(P_diag);
973 P_offd_data = hypre_CSRMatrixData(P_offd);
974 P_offd_i = hypre_CSRMatrixI(P_offd);
975 P_offd_j = hypre_CSRMatrixJ(P_offd);
976 P_diag_size = P_diag_i[n_fine];
977 P_offd_size = P_offd_i[n_fine];
978 }
979
980 /* This builds col_map, col_map should be monotone increasing and contain
981 * global numbers. */
982 num_cols_P_offd = 0;
983 if(P_offd_size)
984 {
985 num_cols_P_offd = 0;
986 for (i=0; i < P_offd_size; i++)
987 {
988 index = P_offd_j[i];
989 if(tmp_CF_marker_offd[index] == 1)
990 {
991 num_cols_P_offd++;
992 tmp_CF_marker_offd[index] = 2;
993 }
994 }
995
996 if (num_cols_P_offd)
997 {
998 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
999 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
1000 }
1001
1002 index = 0;
1003 for(i = 0; i < num_cols_P_offd; i++)
1004 {
1005 while( tmp_CF_marker_offd[index] == -1) index++;
1006 col_map_offd_P[i] = fine_to_coarse_offd[index];
1007 tmp_map_offd[i] = index++;
1008 }
1009
1010 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
1011
1012 for (i = 0; i < num_cols_P_offd; i++)
1013 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
1014
1015 for(i = 0; i < P_offd_size; i++)
1016 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
1017 /*index = 0;
1018 for(i = 0; i < num_cols_P_offd; i++)
1019 {
1020 while (P_marker[index] == 0) index++;
1021
1022 col_map_offd_P[i] = fine_to_coarse_offd[index];
1023 index++;
1024 }*/
1025
1026 /* Sort the col_map_offd_P and P_offd_j correctly */
1027 /* Check if sort actually changed anything */
1028 /*for(i = 0; i < num_cols_P_offd; i++)
1029 P_marker[i] = col_map_offd_P[i];
1030
1031 if(hypre_ssort(col_map_offd_P,num_cols_P_offd))
1032 {
1033 for(i = 0; i < P_offd_size; i++)
1034 for(j = 0; j < num_cols_P_offd; j++)
1035 if(P_marker[P_offd_j[i]] == col_map_offd_P[j])
1036 {
1037 P_offd_j[i] = j;
1038 j = num_cols_P_offd;
1039 }
1040 }
1041 hypre_TFree(P_marker); */
1042 }
1043
1044 if (num_cols_P_offd)
1045 {
1046 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
1047 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
1048 hypre_TFree(tmp_map_offd);
1049 }
1050
1051#ifdef HYPRE_NO_GLOBAL_PARTITION
1052 hypre_NewCommPkgCreate(P);
1053#else
1054 hypre_MatvecCommPkgCreate(P);
1055#endif
1056
1057 for (i=0; i < n_fine; i++)
1058 if (CF_marker[i] == -3) CF_marker[i] = -1;
1059
1060 *P_ptr = P;
1061
1062 /* Deallocate memory */
1063 hypre_TFree(fine_to_coarse);
1064 hypre_TFree(P_marker);
1065 hypre_TFree(ahat);
1066 hypre_TFree(ihat);
1067 hypre_TFree(ipnt);
1068
1069 if (full_off_procNodes)
1070 {
1071 hypre_TFree(ahat_offd);
1072 hypre_TFree(ihat_offd);
1073 hypre_TFree(ipnt_offd);
1074 }
1075 if (num_procs > 1)
1076 {
1077 hypre_BigCSRMatrixDestroy(Sop);
1078 hypre_BigCSRMatrixDestroy(A_ext);
1079 hypre_TFree(fine_to_coarse_offd);
1080 hypre_TFree(P_marker_offd);
1081 hypre_TFree(CF_marker_offd);
1082 hypre_TFree(tmp_CF_marker_offd);
1083
1084 if(num_functions > 1)
1085 hypre_TFree(dof_func_offd);
1086 hypre_TFree(found);
1087
1088 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
1089
1090 }
1091
1092
1093 return hypre_error_flag;
1094}
1095
1096/*---------------------------------------------------------------------------
1097 * hypre_BoomerAMGBuildExtPIInterp
1098 * Comment:
1099 *--------------------------------------------------------------------------*/
1100int
1101hypre_BoomerAMGBuildExtPIInterp(hypre_ParCSRMatrix *A, int *CF_marker,
1102 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
1103 int num_functions, int *dof_func, int debug_flag,
1104 double trunc_factor, int max_elmts,
1105 int *col_offd_S_to_A,
1106 hypre_ParCSRMatrix **P_ptr)
1107{
1108 /* Communication Variables */
1109 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
1110 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1111
1112
1113 int my_id, num_procs;
1114
1115 /* Variables to store input variables */
1116 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1117 double *A_diag_data = hypre_CSRMatrixData(A_diag);
1118 int *A_diag_i = hypre_CSRMatrixI(A_diag);
1119 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
1120
1121 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1122 double *A_offd_data = hypre_CSRMatrixData(A_offd);
1123 int *A_offd_i = hypre_CSRMatrixI(A_offd);
1124 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
1125
1126 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
1127 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
1128 int n_fine = hypre_CSRMatrixNumRows(A_diag);
1129 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
1130 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
1131 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt) local_numrows;
1132 HYPRE_BigInt total_global_cpts, my_first_cpt;
1133
1134 /* Variables to store strong connection matrix info */
1135 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
1136 int *S_diag_i = hypre_CSRMatrixI(S_diag);
1137 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
1138
1139 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
1140 int *S_offd_i = hypre_CSRMatrixI(S_offd);
1141 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
1142
1143 /* Interpolation matrix P */
1144 hypre_ParCSRMatrix *P;
1145 hypre_CSRMatrix *P_diag;
1146 hypre_CSRMatrix *P_offd;
1147
1148 double *P_diag_data = NULL;
1149 int *P_diag_i, *P_diag_j = NULL;
1150 double *P_offd_data = NULL;
1151 int *P_offd_i, *P_offd_j = NULL;
1152
1153 HYPRE_BigInt *col_map_offd_P;
1154 int *tmp_map_offd;
1155 int P_diag_size;
1156 int P_offd_size;
1157 int *P_marker = NULL;
1158 int *P_marker_offd = NULL;
1159 int *CF_marker_offd = NULL;
1160 int *tmp_CF_marker_offd = NULL;
1161 int *dof_func_offd = NULL;
1162
1163 /* Full row information for columns of A that are off diag*/
1164 hypre_BigCSRMatrix *A_ext;
1165 double *A_ext_data;
1166 int *A_ext_i;
1167 HYPRE_BigInt *big_A_ext_j;
1168 int *A_ext_j;
1169
1170 int *fine_to_coarse = NULL;
1171 HYPRE_BigInt *fine_to_coarse_offd = NULL;
1172 HYPRE_BigInt *found;
1173
1174 int num_cols_P_offd;
1175 int newoff, loc_col;
1176 int A_ext_rows, full_off_procNodes;
1177
1178 hypre_BigCSRMatrix *Sop;
1179 int *Sop_i;
1180 HYPRE_BigInt *big_Sop_j;
1181 int *Sop_j;
1182
1183 int Soprows, sgn;
1184
1185 /* Variables to keep count of interpolatory points */
1186 int *jj_count, *jj_count_offd;
1187 int jj_counter, jj_counter_offd;
1188 int jj_begin_row, jj_end_row;
1189 int jj_begin_row_offd = 0;
1190 int jj_end_row_offd = 0;
1191 int *coarse_counter, coarse_counter_offd;
1192
1193 /* Interpolation weight variables */
1194 double sum, diagonal, distribute;
1195 int strong_f_marker = -2;
1196
1197 /* Loop variables */
1198 int index;
1199 int start_indexing = 0;
1200 int i, i1, i2, j, jj, kk, k1, jj1;
1201
1202 /* Definitions */
1203 double zero = 0.0;
1204 double one = 1.0;
1205 double wall_time;
1206 int size, rest, ns, ne;
1207 int num_threads;
1208 int coarse_shift;
1209
1210
1211 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
1212
1213 if (debug_flag==4) wall_time = time_getWallclockSeconds();
1214
1215 /* BEGIN */
1216 MPI_Comm_size(comm, &num_procs);
1217 MPI_Comm_rank(comm,&my_id);
1218 num_threads = hypre_NumThreads();
1219
1220
1221#ifdef HYPRE_NO_GLOBAL_PARTITION
1222 my_first_cpt = num_cpts_global[0];
1223 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
1224 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
1225#else
1226 my_first_cpt = num_cpts_global[my_id];
1227 total_global_cpts = num_cpts_global[num_procs];
1228#endif
1229
1230 if (!comm_pkg)
1231 {
1232#ifdef HYPRE_NO_GLOBAL_PARTITION
1233 hypre_NewCommPkgCreate(A);
1234#else
1235 hypre_MatvecCommPkgCreate(A);
1236#endif
1237 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1238 }
1239
1240 /* Set up off processor information (specifically for neighbors of
1241 * neighbors */
1242 newoff = 0;
1243 full_off_procNodes = 0;
1244 if (num_procs > 1)
1245 {
1246 /*----------------------------------------------------------------------
1247 * Get the off processors rows for A and S, associated with columns in
1248 * A_offd and S_offd.
1249 *---------------------------------------------------------------------*/
1250 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
1251 A_ext_i = hypre_BigCSRMatrixI(A_ext);
1252 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
1253 A_ext_data = hypre_BigCSRMatrixData(A_ext);
1254 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
1255
1256 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
1257 Sop_i = hypre_BigCSRMatrixI(Sop);
1258 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
1259 Soprows = hypre_BigCSRMatrixNumRows(Sop);
1260
1261 /* Find nodes that are neighbors of neighbors, not found in offd */
1262 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
1263 &A_ext_j,
1264 Soprows, col_map_offd, col_1, col_n,
1265 Sop_i, big_Sop_j, &Sop_j, CF_marker, comm_pkg);
1266
1267 hypre_BigCSRMatrixJ(A_ext) = NULL;
1268 hypre_BigCSRMatrixJ(Sop) = NULL;
1269
1270 if(newoff >= 0)
1271 full_off_procNodes = newoff + num_cols_A_offd;
1272 else
1273 return(1);
1274
1275 /* Possibly add new points and new processors to the comm_pkg, all
1276 * processors need new_comm_pkg */
1277
1278 /* AHB - create a new comm package just for extended info -
1279 this will work better with the assumed partition*/
1280 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
1281 &extend_comm_pkg);
1282
1283 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
1284
1285 if (num_functions > 1 && full_off_procNodes > 0)
1286 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
1287
1288 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
1289 full_off_procNodes, CF_marker_offd);
1290
1291 if(num_functions > 1)
1292 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
1293 full_off_procNodes, dof_func_offd);
1294 }
1295
1296
1297 /*-----------------------------------------------------------------------
1298 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
1299 *-----------------------------------------------------------------------*/
1300
1301 /*-----------------------------------------------------------------------
1302 * Intialize counters and allocate mapping vector.
1303 *-----------------------------------------------------------------------*/
1304 P_diag_i = hypre_CTAlloc(int, n_fine+1);
1305 P_offd_i = hypre_CTAlloc(int, n_fine+1);
1306
1307 if (n_fine)
1308 {
1309 fine_to_coarse = hypre_CTAlloc(int, n_fine);
1310 }
1311
1312 if (full_off_procNodes)
1313 {
1314 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
1315 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
1316 }
1317
1318 for (i=0; i < full_off_procNodes; i++)
1319 {
1320 tmp_CF_marker_offd[i] = -1;
1321 fine_to_coarse_offd[i] = -1;
1322 }
1323
1324 coarse_counter = hypre_CTAlloc(int, num_threads);
1325 jj_count = hypre_CTAlloc(int, num_threads);
1326 jj_count_offd = hypre_CTAlloc(int, num_threads);
1327
1328 jj_counter = start_indexing;
1329 jj_counter_offd = start_indexing;
1330 coarse_counter_offd = 0;
1331
1332 /*-----------------------------------------------------------------------
1333 * Loop over fine grid.
1334 *-----------------------------------------------------------------------*/
1335#define HYPRE_SMP_PRIVATE i,j,i1,jj,kk,k1,ns,ne,size,rest,P_marker,P_marker_offd, loc_col
1336#include "../utilities/hypre_smp_forloop.h"
1337 for (j = 0; j < num_threads; j++)
1338 {
1339 size = n_fine/num_threads;
1340 rest = n_fine - size*num_threads;
1341 jj_count[j] = 0;
1342 jj_count_offd[j] = 0;
1343 coarse_counter[j] = 0;
1344 if (n_fine)
1345 {
1346 P_marker = hypre_CTAlloc(int, n_fine);
1347 for (i=0; i < n_fine; i++)
1348 P_marker[i] = -1;
1349 }
1350
1351 if (full_off_procNodes)
1352 {
1353 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
1354 for (i=0; i < full_off_procNodes; i++)
1355 P_marker_offd[i] = -1;
1356 }
1357
1358 if (j < rest)
1359 {
1360 ns = j*size+j;
1361 ne = (j+1)*size+j+1;
1362 }
1363 else
1364 {
1365 ns = j*size+rest;
1366 ne = (j+1)*size+rest;
1367 }
1368 for (i = ns; i < ne; i++)
1369 {
1370 P_diag_i[i] = jj_count[j];
1371 if (num_procs > 1)
1372 P_offd_i[i] = jj_count_offd[j];
1373
1374 fine_to_coarse[i] = -1;
1375 if (CF_marker[i] >= 0)
1376 {
1377 jj_count[j]++;
1378 fine_to_coarse[i] = coarse_counter[j];
1379 coarse_counter[j]++;
1380 }
1381
1382 /*--------------------------------------------------------------------
1383 * If i is an F-point, interpolation is from the C-points that
1384 * strongly influence i, or C-points that stronly influence F-points
1385 * that strongly influence i.
1386 *--------------------------------------------------------------------*/
1387 else if (CF_marker[i] != -3)
1388 {
1389 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
1390 {
1391 i1 = S_diag_j[jj];
1392 if (CF_marker[i1] > 0)
1393 { /* i1 is a C point */
1394 if (P_marker[i1] < P_diag_i[i])
1395 {
1396 P_marker[i1] = jj_count[j];
1397 jj_count[j]++;
1398 }
1399 }
1400 else if (CF_marker[i1] != -3)
1401 { /* i1 is a F point, loop through it's strong neighbors */
1402 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
1403 {
1404 k1 = S_diag_j[kk];
1405 if (CF_marker[k1] > 0)
1406 {
1407 if(P_marker[k1] < P_diag_i[i])
1408 {
1409 P_marker[k1] = jj_count[j];
1410 jj_count[j]++;
1411 }
1412 }
1413 }
1414 if(num_procs > 1)
1415 {
1416 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
1417 {
1418 if(col_offd_S_to_A)
1419 k1 = col_offd_S_to_A[S_offd_j[kk]];
1420 else
1421 k1 = S_offd_j[kk];
1422 if (CF_marker_offd[k1] > 0)
1423 {
1424 if(P_marker_offd[k1] < P_offd_i[i])
1425 {
1426 tmp_CF_marker_offd[k1] = 1;
1427 P_marker_offd[k1] = jj_count_offd[j];
1428 jj_count_offd[j]++;
1429 }
1430 }
1431 }
1432 }
1433 }
1434 }
1435 /* Look at off diag strong connections of i */
1436 if (num_procs > 1)
1437 {
1438 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
1439 {
1440 i1 = S_offd_j[jj];
1441 if(col_offd_S_to_A)
1442 i1 = col_offd_S_to_A[i1];
1443 if (CF_marker_offd[i1] > 0)
1444 {
1445 if(P_marker_offd[i1] < P_offd_i[i])
1446 {
1447 tmp_CF_marker_offd[i1] = 1;
1448 P_marker_offd[i1] = jj_count_offd[j];
1449 jj_count_offd[j]++;
1450 }
1451 }
1452 else if (CF_marker_offd[i1] != -3)
1453 { /* F point; look at neighbors of i1. Sop contains global col
1454 * numbers and entries that could be in S_diag or S_offd or
1455 * neither. */
1456 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
1457 {
1458 loc_col = Sop_j[kk];
1459 if(loc_col > -1)
1460 { /* In S_diag */
1461 if(CF_marker[loc_col] > 0)
1462 {
1463 if(P_marker[loc_col] < P_diag_i[i])
1464 {
1465 P_marker[loc_col] = jj_count[j];
1466 jj_count[j]++;
1467 }
1468 }
1469 }
1470 else
1471 {
1472 loc_col = -loc_col - 1;
1473 if(CF_marker_offd[loc_col] > 0)
1474 {
1475 if(P_marker_offd[loc_col] < P_offd_i[i])
1476 {
1477 P_marker_offd[loc_col] = jj_count_offd[j];
1478 tmp_CF_marker_offd[loc_col] = 1;
1479 jj_count_offd[j]++;
1480 }
1481 }
1482 }
1483 }
1484 }
1485 }
1486 }
1487 }
1488 }
1489 if (n_fine) hypre_TFree(P_marker);
1490 if (full_off_procNodes) hypre_TFree(P_marker_offd);
1491 }
1492
1493 for (i=0; i < num_threads-1; i++)
1494 {
1495 coarse_counter[i+1] += coarse_counter[i];
1496 jj_count[i+1] += jj_count[i];
1497 jj_count_offd[i+1] += jj_count_offd[i];
1498 }
1499 i = num_threads-1;
1500 jj_counter = jj_count[i];
1501 jj_counter_offd = jj_count_offd[i];
1502
1503 if (debug_flag==4)
1504 {
1505 wall_time = time_getWallclockSeconds() - wall_time;
1506 printf("Proc = %d determine structure %f\n",
1507 my_id, wall_time);
1508 fflush(NULL);
1509 }
1510 /*-----------------------------------------------------------------------
1511 * Allocate arrays.
1512 *-----------------------------------------------------------------------*/
1513
1514 if (debug_flag== 4) wall_time = time_getWallclockSeconds();
1515
1516 P_diag_size = jj_counter;
1517 P_offd_size = jj_counter_offd;
1518 if (P_diag_size)
1519 {
1520 P_diag_j = hypre_CTAlloc(int, P_diag_size);
1521 P_diag_data = hypre_CTAlloc(double, P_diag_size);
1522 }
1523
1524 if (P_offd_size)
1525 {
1526 P_offd_j = hypre_CTAlloc(int, P_offd_size);
1527 P_offd_data = hypre_CTAlloc(double, P_offd_size);
1528 }
1529
1530 P_diag_i[n_fine] = jj_counter;
1531 P_offd_i[n_fine] = jj_counter_offd;
1532
1533 jj_counter = start_indexing;
1534 jj_counter_offd = start_indexing;
1535
1536#define HYPRE_SMP_PRIVATE i,j,ns,ne,size,rest,coarse_shift
1537#include "../utilities/hypre_smp_forloop.h"
1538 for (j = 0; j < num_threads; j++)
1539 {
1540 coarse_shift = 0;
1541 if (j > 0) coarse_shift = coarse_counter[j-1];
1542 size = n_fine/num_threads;
1543 rest = n_fine - size*num_threads;
1544 if (j < rest)
1545 {
1546 ns = j*size+j;
1547 ne = (j+1)*size+j+1;
1548 }
1549 else
1550 {
1551 ns = j*size+rest;
1552 ne = (j+1)*size+rest;
1553 }
1554 for (i = ns; i < ne; i++)
1555 fine_to_coarse[i] += coarse_shift;
1556 }
1557
1558 /* Fine to coarse mapping */
1559 if(num_procs > 1)
1560 {
1561 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
1562 full_off_procNodes, my_first_cpt,
1563 fine_to_coarse_offd);
1564 }
1565
1566 /*-----------------------------------------------------------------------
1567 * Loop over fine grid points.
1568 *-----------------------------------------------------------------------*/
1569#define HYPRE_SMP_PRIVATE i,j,i1,i2,jj,jj1,kk,k1,ns,ne,size,rest,P_marker,P_marker_offd,loc_col,jj_counter,jj_counter_offd,strong_f_marker,sum,sgn,distribute,diagonal,jj_begin_row,jj_end_row,jj_begin_row_offd,jj_end_row_offd
1570#include "../utilities/hypre_smp_forloop.h"
1571 for (j = 0; j < num_threads; j++)
1572 {
1573 size = n_fine/num_threads;
1574 rest = n_fine - size*num_threads;
1575 strong_f_marker = -2;
1576 if (n_fine)
1577 {
1578 P_marker = hypre_CTAlloc(int, n_fine);
1579 for (i=0; i < n_fine; i++)
1580 P_marker[i] = -1;
1581 }
1582
1583 jj_counter = 0;
1584 if (j > 0) jj_counter = jj_count[j-1];
1585 jj_counter_offd = 0;
1586 if (j > 0) jj_counter_offd = jj_count_offd[j-1];
1587
1588 if (full_off_procNodes)
1589 {
1590 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
1591 for (i=0; i < full_off_procNodes; i++)
1592 P_marker_offd[i] = -1;
1593 }
1594
1595 if (j < rest)
1596 {
1597 ns = j*size+j;
1598 ne = (j+1)*size+j+1;
1599 }
1600 else
1601 {
1602 ns = j*size+rest;
1603 ne = (j+1)*size+rest;
1604 }
1605 for (i = ns; i < ne; i++)
1606 {
1607 /*for (i = 0; i < n_fine; i++)
1608 {*/
1609 jj_begin_row = jj_counter;
1610 jj_begin_row_offd = jj_counter_offd;
1611 P_diag_i[i] = jj_counter;
1612 if (num_procs > 1)
1613 P_offd_i[i] = jj_counter_offd;
1614
1615 /*--------------------------------------------------------------------
1616 * If i is a c-point, interpolation is the identity.
1617 *--------------------------------------------------------------------*/
1618
1619 if (CF_marker[i] >= 0)
1620 {
1621 P_diag_j[jj_counter] = fine_to_coarse[i];
1622 P_diag_data[jj_counter] = one;
1623 jj_counter++;
1624 }
1625
1626 /*--------------------------------------------------------------------
1627 * If i is an F-point, build interpolation.
1628 *--------------------------------------------------------------------*/
1629
1630 else if (CF_marker[i] != -3)
1631 {
1632 strong_f_marker--;
1633 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
1634 {
1635 i1 = S_diag_j[jj];
1636
1637 /*--------------------------------------------------------------
1638 * If neighbor i1 is a C-point, set column number in P_diag_j
1639 * and initialize interpolation weight to zero.
1640 *--------------------------------------------------------------*/
1641
1642 if (CF_marker[i1] >= 0)
1643 {
1644 if (P_marker[i1] < jj_begin_row)
1645 {
1646 P_marker[i1] = jj_counter;
1647 P_diag_j[jj_counter] = fine_to_coarse[i1];
1648 P_diag_data[jj_counter] = zero;
1649 jj_counter++;
1650 }
1651 }
1652 else if (CF_marker[i1] != -3)
1653 {
1654 P_marker[i1] = strong_f_marker;
1655 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
1656 {
1657 k1 = S_diag_j[kk];
1658 if (CF_marker[k1] > 0)
1659 {
1660 if(P_marker[k1] < jj_begin_row)
1661 {
1662 P_marker[k1] = jj_counter;
1663 P_diag_j[jj_counter] = fine_to_coarse[k1];
1664 P_diag_data[jj_counter] = zero;
1665 jj_counter++;
1666 }
1667 }
1668 }
1669 if(num_procs > 1)
1670 {
1671 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
1672 {
1673 if(col_offd_S_to_A)
1674 k1 = col_offd_S_to_A[S_offd_j[kk]];
1675 else
1676 k1 = S_offd_j[kk];
1677 if(CF_marker_offd[k1] > 0)
1678 {
1679 if(P_marker_offd[k1] < jj_begin_row_offd)
1680 {
1681 P_marker_offd[k1] = jj_counter_offd;
1682 P_offd_j[jj_counter_offd] = k1;
1683 P_offd_data[jj_counter_offd] = zero;
1684 jj_counter_offd++;
1685 }
1686 }
1687 }
1688 }
1689 }
1690 }
1691
1692 if ( num_procs > 1)
1693 {
1694 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
1695 {
1696 i1 = S_offd_j[jj];
1697 if(col_offd_S_to_A)
1698 i1 = col_offd_S_to_A[i1];
1699 if ( CF_marker_offd[i1] > 0)
1700 {
1701 if(P_marker_offd[i1] < jj_begin_row_offd)
1702 {
1703 P_marker_offd[i1] = jj_counter_offd;
1704 P_offd_j[jj_counter_offd] = i1;
1705 P_offd_data[jj_counter_offd] = zero;
1706 jj_counter_offd++;
1707 }
1708 }
1709 else if (CF_marker_offd[i1] != -3)
1710 {
1711 P_marker_offd[i1] = strong_f_marker;
1712 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
1713 {
1714 loc_col = Sop_j[kk];
1715 if(loc_col > -1)
1716 {
1717 if(CF_marker[loc_col] > 0)
1718 {
1719 if(P_marker[loc_col] < jj_begin_row)
1720 {
1721 P_marker[loc_col] = jj_counter;
1722 P_diag_j[jj_counter] = fine_to_coarse[loc_col];
1723 P_diag_data[jj_counter] = zero;
1724 jj_counter++;
1725 }
1726 }
1727 }
1728 else
1729 {
1730 loc_col = - loc_col - 1;
1731 if(CF_marker_offd[loc_col] > 0)
1732 {
1733 if(P_marker_offd[loc_col] < jj_begin_row_offd)
1734 {
1735 P_marker_offd[loc_col] = jj_counter_offd;
1736 P_offd_j[jj_counter_offd]=loc_col;
1737 P_offd_data[jj_counter_offd] = zero;
1738 jj_counter_offd++;
1739 }
1740 }
1741 }
1742 }
1743 }
1744 }
1745 }
1746
1747 jj_end_row = jj_counter;
1748 jj_end_row_offd = jj_counter_offd;
1749
1750 diagonal = A_diag_data[A_diag_i[i]];
1751
1752 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
1753 { /* i1 is a c-point and strongly influences i, accumulate
1754 * a_(i,i1) into interpolation weight */
1755 i1 = A_diag_j[jj];
1756 if (P_marker[i1] >= jj_begin_row)
1757 {
1758 P_diag_data[P_marker[i1]] += A_diag_data[jj];
1759 }
1760 else if(P_marker[i1] == strong_f_marker)
1761 {
1762 sum = zero;
1763 sgn = 1;
1764 if(A_diag_data[A_diag_i[i1]] < 0) sgn = -1;
1765 /* Loop over row of A for point i1 and calculate the sum
1766 * of the connections to c-points that strongly incluence i. */
1767 for(jj1 = A_diag_i[i1]+1; jj1 < A_diag_i[i1+1]; jj1++)
1768 {
1769 i2 = A_diag_j[jj1];
1770 if((P_marker[i2] >= jj_begin_row || i2 == i) && (sgn*A_diag_data[jj1]) < 0)
1771 sum += A_diag_data[jj1];
1772 }
1773 if(num_procs > 1)
1774 {
1775 for(jj1 = A_offd_i[i1]; jj1< A_offd_i[i1+1]; jj1++)
1776 {
1777 i2 = A_offd_j[jj1];
1778 if(P_marker_offd[i2] >= jj_begin_row_offd &&
1779 (sgn*A_offd_data[jj1]) < 0)
1780 sum += A_offd_data[jj1];
1781 }
1782 }
1783 if(sum != 0)
1784 {
1785 distribute = A_diag_data[jj]/sum;
1786 /* Loop over row of A for point i1 and do the distribution */
1787 for(jj1 = A_diag_i[i1]+1; jj1 < A_diag_i[i1+1]; jj1++)
1788 {
1789 i2 = A_diag_j[jj1];
1790 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
1791 P_diag_data[P_marker[i2]] +=
1792 distribute*A_diag_data[jj1];
1793 if(i2 == i && (sgn*A_diag_data[jj1]) < 0)
1794 diagonal += distribute*A_diag_data[jj1];
1795 }
1796 if(num_procs > 1)
1797 {
1798 for(jj1 = A_offd_i[i1]; jj1 < A_offd_i[i1+1]; jj1++)
1799 {
1800 i2 = A_offd_j[jj1];
1801 if(P_marker_offd[i2] >= jj_begin_row_offd &&
1802 (sgn*A_offd_data[jj1]) < 0)
1803 P_offd_data[P_marker_offd[i2]] +=
1804 distribute*A_offd_data[jj1];
1805 }
1806 }
1807 }
1808 else
1809 {
1810 diagonal += A_diag_data[jj];
1811 }
1812 }
1813 /* neighbor i1 weakly influences i, accumulate a_(i,i1) into
1814 * diagonal */
1815 else if (CF_marker[i1] != -3)
1816 {
1817 if(num_functions == 1 || dof_func[i] == dof_func[i1])
1818 diagonal += A_diag_data[jj];
1819 }
1820 }
1821 if(num_procs > 1)
1822 {
1823 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
1824 {
1825 i1 = A_offd_j[jj];
1826 if(P_marker_offd[i1] >= jj_begin_row_offd)
1827 P_offd_data[P_marker_offd[i1]] += A_offd_data[jj];
1828 else if(P_marker_offd[i1] == strong_f_marker)
1829 {
1830 sum = zero;
1831 sgn = 1;
1832 if(A_ext_data[A_ext_i[i1]] < 0) sgn = -1;
1833 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
1834 {
1835 loc_col = A_ext_j[jj1];
1836 if(loc_col > -1)
1837 { /* diag */
1838 if((P_marker[loc_col] >= jj_begin_row || loc_col == i)
1839 && (sgn*A_ext_data[jj1]) < 0)
1840 sum += A_ext_data[jj1];
1841 }
1842 else
1843 {
1844 loc_col = - loc_col - 1;
1845 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
1846 (sgn*A_ext_data[jj1]) < 0)
1847 sum += A_ext_data[jj1];
1848 }
1849 }
1850 if(sum != 0)
1851 {
1852 distribute = A_offd_data[jj] / sum;
1853 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
1854 {
1855 loc_col = A_ext_j[jj1];
1856 if(loc_col > -1)
1857 { /* diag */
1858 if(P_marker[loc_col] >= jj_begin_row &&
1859 (sgn*A_ext_data[jj1]) < 0)
1860 P_diag_data[P_marker[loc_col]] += distribute*
1861 A_ext_data[jj1];
1862 if(loc_col == i && (sgn*A_ext_data[jj1]) < 0)
1863 diagonal += distribute*A_ext_data[jj1];
1864 }
1865 else
1866 {
1867 loc_col = - loc_col - 1;
1868 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
1869 (sgn*A_ext_data[jj1]) < 0)
1870 P_offd_data[P_marker_offd[loc_col]] += distribute*
1871 A_ext_data[jj1];
1872 }
1873 }
1874 }
1875 else
1876 {
1877 diagonal += A_offd_data[jj];
1878 }
1879 }
1880 else if (CF_marker_offd[i1] != -3)
1881 {
1882 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
1883 diagonal += A_offd_data[jj];
1884 }
1885 }
1886 }
1887 for(jj = jj_begin_row; jj < jj_end_row; jj++)
1888 P_diag_data[jj] /= -diagonal;
1889 for(jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
1890 P_offd_data[jj] /= -diagonal;
1891 }
1892 strong_f_marker--;
1893 }
1894 if (n_fine) hypre_TFree(P_marker);
1895 if (full_off_procNodes) hypre_TFree(P_marker_offd);
1896 }
1897
1898 if (debug_flag==4)
1899 {
1900 wall_time = time_getWallclockSeconds() - wall_time;
1901 printf("Proc = %d fill structure %f\n",
1902 my_id, wall_time);
1903 fflush(NULL);
1904 }
1905 /*-----------------------------------------------------------------------
1906 * Allocate arrays.
1907 *-----------------------------------------------------------------------*/
1908
1909 P = hypre_ParCSRMatrixCreate(comm,
1910 hypre_ParCSRMatrixGlobalNumRows(A),
1911 total_global_cpts,
1912 hypre_ParCSRMatrixColStarts(A),
1913 num_cpts_global,
1914 0,
1915 P_diag_i[n_fine],
1916 P_offd_i[n_fine]);
1917
1918 P_diag = hypre_ParCSRMatrixDiag(P);
1919 hypre_CSRMatrixData(P_diag) = P_diag_data;
1920 hypre_CSRMatrixI(P_diag) = P_diag_i;
1921 hypre_CSRMatrixJ(P_diag) = P_diag_j;
1922 P_offd = hypre_ParCSRMatrixOffd(P);
1923 hypre_CSRMatrixData(P_offd) = P_offd_data;
1924 hypre_CSRMatrixI(P_offd) = P_offd_i;
1925 hypre_CSRMatrixJ(P_offd) = P_offd_j;
1926 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
1927
1928 /* Compress P, removing coefficients smaller than trunc_factor * Max */
1929 if (trunc_factor != 0.0 || max_elmts > 0)
1930 {
1931 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
1932 P_diag_data = hypre_CSRMatrixData(P_diag);
1933 P_diag_i = hypre_CSRMatrixI(P_diag);
1934 P_diag_j = hypre_CSRMatrixJ(P_diag);
1935 P_offd_data = hypre_CSRMatrixData(P_offd);
1936 P_offd_i = hypre_CSRMatrixI(P_offd);
1937 P_offd_j = hypre_CSRMatrixJ(P_offd);
1938 P_diag_size = P_diag_i[n_fine];
1939 P_offd_size = P_offd_i[n_fine];
1940 }
1941
1942 /* This builds col_map, col_map should be monotone increasing and contain
1943 * global numbers. */
1944 num_cols_P_offd = 0;
1945 if(P_offd_size)
1946 {
1947 num_cols_P_offd = 0;
1948 for (i=0; i < P_offd_size; i++)
1949 {
1950 index = P_offd_j[i];
1951 if (tmp_CF_marker_offd[index] == 1)
1952 {
1953 num_cols_P_offd++;
1954 tmp_CF_marker_offd[index] = 2;
1955 }
1956 }
1957
1958 if (num_cols_P_offd)
1959 {
1960 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
1961 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
1962 }
1963
1964 index = 0;
1965 for(i = 0; i < num_cols_P_offd; i++)
1966 {
1967 while( tmp_CF_marker_offd[index] != 2) index++;
1968 /*printf("index %d tmp_CF %d\n", index, tmp_CF_marker_offd[index]);*/
1969 col_map_offd_P[i] = fine_to_coarse_offd[index];
1970 tmp_map_offd[i] = index++;
1971 }
1972
1973 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
1974
1975 for (i = 0; i < num_cols_P_offd; i++)
1976 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
1977
1978 for(i = 0; i < P_offd_size; i++)
1979 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
1980 }
1981
1982 if (num_cols_P_offd)
1983 {
1984 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
1985 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
1986 hypre_TFree(tmp_map_offd);
1987 }
1988
1989#ifdef HYPRE_NO_GLOBAL_PARTITION
1990 hypre_NewCommPkgCreate(P);
1991#else
1992 hypre_MatvecCommPkgCreate(P);
1993#endif
1994
1995 for (i=0; i < n_fine; i++)
1996 if (CF_marker[i] == -3) CF_marker[i] = -1;
1997
1998 *P_ptr = P;
1999
2000 /* Deallocate memory */
2001 hypre_TFree(fine_to_coarse);
2002 hypre_TFree(jj_count);
2003 hypre_TFree(jj_count_offd);
2004 hypre_TFree(coarse_counter);
2005
2006 if (num_procs > 1)
2007 {
2008 hypre_BigCSRMatrixDestroy(Sop);
2009 hypre_BigCSRMatrixDestroy(A_ext);
2010 hypre_TFree(fine_to_coarse_offd);
2011 hypre_TFree(CF_marker_offd);
2012 hypre_TFree(tmp_CF_marker_offd);
2013 if(num_functions > 1)
2014 hypre_TFree(dof_func_offd);
2015 hypre_TFree(found);
2016
2017
2018 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
2019
2020
2021 }
2022
2023 return hypre_error_flag;
2024}
2025
2026/*---------------------------------------------------------------------------
2027 * hypre_BoomerAMGBuildExtPICCInterp
2028 * Comment: Only use FF when there is no common c point.
2029 *--------------------------------------------------------------------------*/
2030int
2031hypre_BoomerAMGBuildExtPICCInterp(hypre_ParCSRMatrix *A, int *CF_marker,
2032 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
2033 int num_functions, int *dof_func, int debug_flag,
2034 double trunc_factor, int max_elmts,
2035 int *col_offd_S_to_A,
2036 hypre_ParCSRMatrix **P_ptr)
2037{
2038 /* Communication Variables */
2039 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
2040 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
2041
2042
2043 int my_id, num_procs;
2044
2045 /* Variables to store input variables */
2046 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
2047 double *A_diag_data = hypre_CSRMatrixData(A_diag);
2048 int *A_diag_i = hypre_CSRMatrixI(A_diag);
2049 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
2050
2051 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
2052 double *A_offd_data = hypre_CSRMatrixData(A_offd);
2053 int *A_offd_i = hypre_CSRMatrixI(A_offd);
2054 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
2055
2056 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
2057 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
2058 int n_fine = hypre_CSRMatrixNumRows(A_diag);
2059 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
2060 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
2061 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt) local_numrows;
2062 HYPRE_BigInt total_global_cpts, my_first_cpt;
2063
2064 /* Variables to store strong connection matrix info */
2065 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
2066 int *S_diag_i = hypre_CSRMatrixI(S_diag);
2067 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
2068
2069 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
2070 int *S_offd_i = hypre_CSRMatrixI(S_offd);
2071 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
2072
2073 /* Interpolation matrix P */
2074 hypre_ParCSRMatrix *P;
2075 hypre_CSRMatrix *P_diag;
2076 hypre_CSRMatrix *P_offd;
2077
2078 double *P_diag_data = NULL;
2079 int *P_diag_i, *P_diag_j = NULL;
2080 double *P_offd_data = NULL;
2081 int *P_offd_i, *P_offd_j = NULL;
2082
2083 HYPRE_BigInt *col_map_offd_P;
2084 int *tmp_map_offd;
2085 int P_diag_size;
2086 int P_offd_size;
2087 int *P_marker = NULL;
2088 int *P_marker_offd = NULL;
2089 int *CF_marker_offd = NULL;
2090 int *tmp_CF_marker_offd = NULL;
2091 int *dof_func_offd = NULL;
2092 /*int **ext_p, **ext_p_offd;*/
2093 int ccounter_offd;
2094 /*int *clist_offd;*/
2095 int common_c;
2096
2097 /* Full row information for columns of A that are off diag*/
2098 hypre_BigCSRMatrix *A_ext;
2099 double *A_ext_data;
2100 int *A_ext_i;
2101 HYPRE_BigInt *big_A_ext_j;
2102 int *A_ext_j;
2103
2104 int *fine_to_coarse = NULL;
2105 HYPRE_BigInt *fine_to_coarse_offd = NULL;
2106 HYPRE_BigInt *found;
2107
2108 int num_cols_P_offd;
2109 int newoff, loc_col;
2110 int A_ext_rows, full_off_procNodes;
2111
2112 hypre_BigCSRMatrix *Sop;
2113 int *Sop_i;
2114 HYPRE_BigInt *big_Sop_j;
2115 int *Sop_j;
2116
2117 int Soprows, sgn;
2118
2119 /* Variables to keep count of interpolatory points */
2120 int jj_counter, jj_counter_offd;
2121 int jj_begin_row, jj_end_row;
2122 int jj_begin_row_offd = 0;
2123 int jj_end_row_offd = 0;
2124 int coarse_counter, coarse_counter_offd;
2125
2126 /* Interpolation weight variables */
2127 double sum, diagonal, distribute;
2128 int strong_f_marker = -2;
2129
2130 /* Loop variables */
2131 int index;
2132 int start_indexing = 0;
2133 int i, i1, i2, jj, kk, k1, jj1;
2134 int ccounter;
2135 /*int *clist, ccounter;*/
2136
2137 /* Definitions */
2138 double zero = 0.0;
2139 double one = 1.0;
2140
2141 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
2142
2143 /* BEGIN */
2144 MPI_Comm_size(comm, &num_procs);
2145 MPI_Comm_rank(comm,&my_id);
2146
2147#ifdef HYPRE_NO_GLOBAL_PARTITION
2148 my_first_cpt = num_cpts_global[0];
2149 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
2150 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
2151#else
2152 my_first_cpt = num_cpts_global[my_id];
2153 total_global_cpts = num_cpts_global[num_procs];
2154#endif
2155
2156 if (!comm_pkg)
2157 {
2158#ifdef HYPRE_NO_GLOBAL_PARTITION
2159 hypre_NewCommPkgCreate(A);
2160#else
2161 hypre_MatvecCommPkgCreate(A);
2162#endif
2163 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
2164 }
2165
2166 /* Set up off processor information (specifically for neighbors of
2167 * neighbors */
2168 newoff = 0;
2169 full_off_procNodes = 0;
2170 if (num_procs > 1)
2171 {
2172 /*----------------------------------------------------------------------
2173 * Get the off processors rows for A and S, associated with columns in
2174 * A_offd and S_offd.
2175 *---------------------------------------------------------------------*/
2176 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
2177 A_ext_i = hypre_BigCSRMatrixI(A_ext);
2178 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
2179 A_ext_data = hypre_BigCSRMatrixData(A_ext);
2180 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
2181
2182 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
2183 Sop_i = hypre_BigCSRMatrixI(Sop);
2184 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
2185 Soprows = hypre_BigCSRMatrixNumRows(Sop);
2186
2187 /* Find nodes that are neighbors of neighbors, not found in offd */
2188 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
2189 &A_ext_j,
2190 Soprows, col_map_offd, col_1, col_n,
2191 Sop_i, big_Sop_j, &Sop_j, CF_marker, comm_pkg);
2192
2193 hypre_BigCSRMatrixJ(A_ext) = NULL;
2194 hypre_BigCSRMatrixJ(Sop) = NULL;
2195
2196 if(newoff >= 0)
2197 full_off_procNodes = newoff + num_cols_A_offd;
2198 else
2199 return(1);
2200
2201 /* Possibly add new points and new processors to the comm_pkg, all
2202 * processors need new_comm_pkg */
2203
2204 /* AHB - create a new comm package just for extended info -
2205 this will work better with the assumed partition*/
2206 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
2207 &extend_comm_pkg);
2208
2209 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
2210
2211 if (num_functions > 1 && full_off_procNodes > 0)
2212 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
2213
2214 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
2215 full_off_procNodes, CF_marker_offd);
2216
2217 if(num_functions > 1)
2218 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
2219 full_off_procNodes, dof_func_offd);
2220
2221
2222 }
2223
2224 /*-----------------------------------------------------------------------
2225 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
2226 *-----------------------------------------------------------------------*/
2227
2228 /*-----------------------------------------------------------------------
2229 * Intialize counters and allocate mapping vector.
2230 *-----------------------------------------------------------------------*/
2231 P_diag_i = hypre_CTAlloc(int, n_fine+1);
2232 P_offd_i = hypre_CTAlloc(int, n_fine+1);
2233
2234 if (n_fine)
2235 {
2236 fine_to_coarse = hypre_CTAlloc(int, n_fine);
2237 P_marker = hypre_CTAlloc(int, n_fine);
2238 }
2239
2240 if (full_off_procNodes)
2241 {
2242 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
2243 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
2244 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
2245 }
2246
2247 for (i=0; i < full_off_procNodes; i++)
2248 {
2249 P_marker_offd[i] = -1;
2250 tmp_CF_marker_offd[i] = -1;
2251 fine_to_coarse_offd[i] = -1;
2252 }
2253
2254 for (i=0; i < n_fine; i++)
2255 {
2256 P_marker[i] = -1;
2257 fine_to_coarse[i] = -1;
2258 }
2259
2260 jj_counter = start_indexing;
2261 jj_counter_offd = start_indexing;
2262 coarse_counter = 0;
2263 coarse_counter_offd = 0;
2264
2265 /*-----------------------------------------------------------------------
2266 * Loop over fine grid.
2267 *-----------------------------------------------------------------------*/
2268 for (i = 0; i < n_fine; i++)
2269 {
2270 P_diag_i[i] = jj_counter;
2271 if (num_procs > 1)
2272 P_offd_i[i] = jj_counter_offd;
2273
2274 if (CF_marker[i] > 0)
2275 {
2276 jj_counter++;
2277 fine_to_coarse[i] = coarse_counter;
2278 coarse_counter++;
2279 }
2280
2281 /*--------------------------------------------------------------------
2282 * If i is an F-point, interpolation is from the C-points that
2283 * strongly influence i, or C-points that stronly influence F-points
2284 * that strongly influence i.
2285 *--------------------------------------------------------------------*/
2286 else if (CF_marker[i] != -3)
2287 {
2288 /* Initialize ccounter for each f point */
2289 ccounter = 0;
2290 ccounter_offd = 0;
2291 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2292 { /* search through diag to find all c neighbors */
2293 i1 = S_diag_j[jj];
2294 if (CF_marker[i1] > 0)
2295 { /* i1 is a C point */
2296 CF_marker[i1] = 2;
2297 if (P_marker[i1] < P_diag_i[i])
2298 {
2299 /*clist[ccounter++] = i1;*/
2300 P_marker[i1] = jj_counter;
2301 jj_counter++;
2302 }
2303 }
2304 }
2305 /*qsort0(clist,0,ccounter-1);*/
2306 if(num_procs > 1)
2307 {
2308 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2309 { /* search through offd to find all c neighbors */
2310 if(col_offd_S_to_A)
2311 i1 = col_offd_S_to_A[S_offd_j[jj]];
2312 else
2313 i1 = S_offd_j[jj];
2314 if(CF_marker_offd[i1] > 0)
2315 { /* i1 is a C point direct neighbor */
2316 CF_marker_offd[i1] = 2;
2317 if(P_marker_offd[i1] < P_offd_i[i])
2318 {
2319 /*clist_offd[ccounter_offd++] = i1;*/
2320 tmp_CF_marker_offd[i1] = 1;
2321 P_marker_offd[i1] = jj_counter_offd;
2322 jj_counter_offd++;
2323 }
2324 }
2325 }
2326 /*qsort0(clist_offd,0,ccounter_offd-1);*/
2327 }
2328 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2329 { /* Search diag to find f neighbors and determine if common c point */
2330 i1 = S_diag_j[jj];
2331 if (CF_marker[i1] == -1)
2332 { /* i1 is a F point, loop through it's strong neighbors */
2333 common_c = 0;
2334 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
2335 {
2336 k1 = S_diag_j[kk];
2337 if (CF_marker[k1] == 2)
2338 {
2339 /*if(hypre_BinarySearch(clist,k1,ccounter) >= 0)
2340 {*/
2341 common_c = 1;
2342 break;
2343 /*kk = S_diag_i[i1+1];
2344 }*/
2345 }
2346 }
2347 if(num_procs > 1 && common_c == 0)
2348 { /* no common c point yet, check offd */
2349 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
2350 {
2351 if(col_offd_S_to_A)
2352 k1 = col_offd_S_to_A[S_offd_j[kk]];
2353 else
2354 k1 = S_offd_j[kk];
2355
2356 if (CF_marker_offd[k1] == 2)
2357 { /* k1 is a c point check if it is common */
2358 /*if(hypre_BinarySearch(clist_offd,k1,ccounter_offd) >= 0)
2359 {*/
2360 common_c = 1;
2361 break;
2362 /*kk = S_offd_i[i1+1];
2363 }*/
2364 }
2365 }
2366 }
2367 if(!common_c)
2368 { /* No common c point, extend the interp set */
2369 for(kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
2370 {
2371 k1 = S_diag_j[kk];
2372 if(CF_marker[k1] > 0)
2373 {
2374 if(P_marker[k1] < P_diag_i[i])
2375 {
2376 P_marker[k1] = jj_counter;
2377 jj_counter++;
2378 /*break;*/
2379 }
2380 }
2381 }
2382 if(num_procs > 1)
2383 {
2384 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
2385 {
2386 if(col_offd_S_to_A)
2387 k1 = col_offd_S_to_A[S_offd_j[kk]];
2388 else
2389 k1 = S_offd_j[kk];
2390 if (CF_marker_offd[k1] > 0)
2391 {
2392 if(P_marker_offd[k1] < P_offd_i[i])
2393 {
2394 tmp_CF_marker_offd[k1] = 1;
2395 P_marker_offd[k1] = jj_counter_offd;
2396 jj_counter_offd++;
2397 /*break;*/
2398 }
2399 }
2400 }
2401 }
2402 }
2403 }
2404 }
2405 /* Look at off diag strong connections of i */
2406 if (num_procs > 1)
2407 {
2408 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2409 {
2410 i1 = S_offd_j[jj];
2411 if(col_offd_S_to_A)
2412 i1 = col_offd_S_to_A[i1];
2413 if (CF_marker_offd[i1] == -1)
2414 { /* F point; look at neighbors of i1. Sop contains global col
2415 * numbers and entries that could be in S_diag or S_offd or
2416 * neither. */
2417 common_c = 0;
2418 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
2419 { /* Check if common c */
2420 loc_col = Sop_j[kk];
2421 if(loc_col > -1)
2422 { /* In S_diag */
2423 if(CF_marker[loc_col] == 2)
2424 {
2425 common_c = 1;
2426 break;
2427 }
2428 }
2429 else
2430 {
2431 loc_col = -loc_col - 1;
2432 if(CF_marker_offd[loc_col] == 2)
2433 {
2434 common_c = 1;
2435 break;
2436 }
2437 }
2438 }
2439 if(!common_c)
2440 {
2441 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
2442 { /* Check if common c */
2443 loc_col = Sop_j[kk];
2444 if(loc_col > -1)
2445 { /* In S_diag */
2446 if(CF_marker[loc_col] > 0)
2447 {
2448 if(P_marker[loc_col] < P_diag_i[i])
2449 {
2450 P_marker[loc_col] = jj_counter;
2451 jj_counter++;
2452 }
2453 }
2454 }
2455 else
2456 {
2457 loc_col = - loc_col - 1;
2458 if(CF_marker_offd[loc_col] > 0)
2459 {
2460 if(P_marker_offd[loc_col] < P_offd_i[i])
2461 {
2462 P_marker_offd[loc_col] = jj_counter_offd;
2463 tmp_CF_marker_offd[loc_col] = 1;
2464 jj_counter_offd++;
2465 }
2466 }
2467 }
2468 }
2469 }
2470 }
2471 }
2472 }
2473 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2474 { /* search through diag to find all c neighbors */
2475 i1 = S_diag_j[jj];
2476 if (CF_marker[i1] == 2)
2477 CF_marker[i1] = 1;
2478 }
2479 if(num_procs > 1)
2480 {
2481 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2482 { /* search through offd to find all c neighbors */
2483 if(col_offd_S_to_A)
2484 i1 = col_offd_S_to_A[S_offd_j[jj]];
2485 else
2486 i1 = S_offd_j[jj];
2487 if(CF_marker_offd[i1] == 2)
2488 { /* i1 is a C point direct neighbor */
2489 CF_marker_offd[i1] = 1;
2490 }
2491 }
2492 }
2493 }
2494 }
2495
2496 /*-----------------------------------------------------------------------
2497 * Allocate arrays.
2498 *-----------------------------------------------------------------------*/
2499
2500 P_diag_size = jj_counter;
2501 P_offd_size = jj_counter_offd;
2502 if (P_diag_size)
2503 {
2504 P_diag_j = hypre_CTAlloc(int, P_diag_size);
2505 P_diag_data = hypre_CTAlloc(double, P_diag_size);
2506 }
2507
2508 if (P_offd_size)
2509 {
2510 P_offd_j = hypre_CTAlloc(int, P_offd_size);
2511 P_offd_data = hypre_CTAlloc(double, P_offd_size);
2512 }
2513 P_diag_i[n_fine] = jj_counter;
2514 P_offd_i[n_fine] = jj_counter_offd;
2515
2516 jj_counter = start_indexing;
2517 jj_counter_offd = start_indexing;
2518 ccounter = start_indexing;
2519 ccounter_offd = start_indexing;
2520
2521 /* Fine to coarse mapping */
2522 if(num_procs > 1)
2523 {
2524 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
2525 full_off_procNodes, my_first_cpt,
2526 fine_to_coarse_offd);
2527 }
2528
2529 for (i = 0; i < n_fine; i++)
2530 P_marker[i] = -1;
2531
2532 for (i = 0; i < full_off_procNodes; i++)
2533 P_marker_offd[i] = -1;
2534
2535 /*-----------------------------------------------------------------------
2536 * Loop over fine grid points.
2537 *-----------------------------------------------------------------------*/
2538 for (i = 0; i < n_fine; i++)
2539 {
2540 jj_begin_row = jj_counter;
2541 if(num_procs > 1)
2542 jj_begin_row_offd = jj_counter_offd;
2543
2544 /*--------------------------------------------------------------------
2545 * If i is a c-point, interpolation is the identity.
2546 *--------------------------------------------------------------------*/
2547
2548 if (CF_marker[i] >= 0)
2549 {
2550 P_diag_j[jj_counter] = fine_to_coarse[i];
2551 P_diag_data[jj_counter] = one;
2552 jj_counter++;
2553 }
2554
2555 /*--------------------------------------------------------------------
2556 * If i is an F-point, build interpolation.
2557 *--------------------------------------------------------------------*/
2558
2559 else if (CF_marker[i] != -3)
2560 {
2561 ccounter = 0;
2562 ccounter_offd = 0;
2563 strong_f_marker--;
2564
2565 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2566 { /* Search C points only */
2567 i1 = S_diag_j[jj];
2568
2569 /*--------------------------------------------------------------
2570 * If neighbor i1 is a C-point, set column number in P_diag_j
2571 * and initialize interpolation weight to zero.
2572 *--------------------------------------------------------------*/
2573
2574 if (CF_marker[i1] > 0)
2575 {
2576 CF_marker[i1] = 2;
2577 if (P_marker[i1] < jj_begin_row)
2578 {
2579 P_marker[i1] = jj_counter;
2580 P_diag_j[jj_counter] = fine_to_coarse[i1];
2581 P_diag_data[jj_counter] = zero;
2582 jj_counter++;
2583 }
2584 }
2585 }
2586 if ( num_procs > 1)
2587 {
2588 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2589 {
2590 if(col_offd_S_to_A)
2591 i1 = col_offd_S_to_A[S_offd_j[jj]];
2592 else
2593 i1 = S_offd_j[jj];
2594 if ( CF_marker_offd[i1] > 0)
2595 {
2596 CF_marker_offd[i1] = 2;
2597 if(P_marker_offd[i1] < jj_begin_row_offd)
2598 {
2599 P_marker_offd[i1] = jj_counter_offd;
2600 P_offd_j[jj_counter_offd] = i1;
2601 P_offd_data[jj_counter_offd] = zero;
2602 jj_counter_offd++;
2603 }
2604 }
2605 }
2606 }
2607
2608 for(jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2609 { /* Search through F points */
2610 i1 = S_diag_j[jj];
2611 if(CF_marker[i1] == -1)
2612 {
2613 P_marker[i1] = strong_f_marker;
2614 common_c = 0;
2615 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
2616 {
2617 k1 = S_diag_j[kk];
2618 if (CF_marker[k1] == 2)
2619 {
2620 common_c = 1;
2621 break;
2622 }
2623 }
2624 if(num_procs > 1 && common_c == 0)
2625 { /* no common c point yet, check offd */
2626 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
2627 {
2628 if(col_offd_S_to_A)
2629 k1 = col_offd_S_to_A[S_offd_j[kk]];
2630 else
2631 k1 = S_offd_j[kk];
2632
2633 if (CF_marker_offd[k1] == 2)
2634 { /* k1 is a c point check if it is common */
2635 common_c = 1;
2636 break;
2637 }
2638 }
2639 }
2640 if(!common_c)
2641 { /* No common c point, extend the interp set */
2642 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
2643 {
2644 k1 = S_diag_j[kk];
2645 if (CF_marker[k1] >= 0)
2646 {
2647 if(P_marker[k1] < jj_begin_row)
2648 {
2649 P_marker[k1] = jj_counter;
2650 P_diag_j[jj_counter] = fine_to_coarse[k1];
2651 P_diag_data[jj_counter] = zero;
2652 jj_counter++;
2653 }
2654 }
2655 }
2656 if(num_procs > 1)
2657 {
2658 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
2659 {
2660 if(col_offd_S_to_A)
2661 k1 = col_offd_S_to_A[S_offd_j[kk]];
2662 else
2663 k1 = S_offd_j[kk];
2664 if(CF_marker_offd[k1] >= 0)
2665 {
2666 if(P_marker_offd[k1] < jj_begin_row_offd)
2667 {
2668 P_marker_offd[k1] = jj_counter_offd;
2669 P_offd_j[jj_counter_offd] = k1;
2670 P_offd_data[jj_counter_offd] = zero;
2671 jj_counter_offd++;
2672 }
2673 }
2674 }
2675 }
2676 }
2677 }
2678 }
2679 if ( num_procs > 1)
2680 {
2681 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2682 {
2683 i1 = S_offd_j[jj];
2684 if(col_offd_S_to_A)
2685 i1 = col_offd_S_to_A[i1];
2686 if(CF_marker_offd[i1] == -1)
2687 { /* F points that are off proc */
2688 P_marker_offd[i1] = strong_f_marker;
2689 common_c = 0;
2690 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
2691 { /* Check if common c */
2692 loc_col = Sop_j[kk];
2693 if(loc_col > -1)
2694 { /* In S_diag */
2695 if(CF_marker[loc_col] == 2)
2696 {
2697 common_c = 1;
2698 break;
2699 }
2700 }
2701 else
2702 {
2703 loc_col = -loc_col - 1;
2704 if(CF_marker_offd[loc_col] == 2)
2705 {
2706 common_c = 1;
2707 break;
2708 }
2709 }
2710 }
2711 if(!common_c)
2712 {
2713 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
2714 {
2715 loc_col = Sop_j[kk];
2716 /* Find local col number */
2717 if(loc_col > -1)
2718 {
2719 if(CF_marker[loc_col] > 0)
2720 {
2721 if(P_marker[loc_col] < jj_begin_row)
2722 {
2723 P_marker[loc_col] = jj_counter;
2724 P_diag_j[jj_counter] = fine_to_coarse[loc_col];
2725 P_diag_data[jj_counter] = zero;
2726 jj_counter++;
2727 }
2728 }
2729 }
2730 else
2731 {
2732 loc_col = -loc_col - 1;
2733 if(CF_marker_offd[loc_col] > 0)
2734 {
2735 if(P_marker_offd[loc_col] < jj_begin_row_offd)
2736 {
2737 P_marker_offd[loc_col] = jj_counter_offd;
2738 P_offd_j[jj_counter_offd]=loc_col;
2739 P_offd_data[jj_counter_offd] = zero;
2740 jj_counter_offd++;
2741 }
2742 }
2743 }
2744 }
2745 }
2746 }
2747 }
2748 }
2749 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
2750 { /* Search C points only */
2751 i1 = S_diag_j[jj];
2752
2753 /*--------------------------------------------------------------
2754 * If neighbor i1 is a C-point, set column number in P_diag_j
2755 * and initialize interpolation weight to zero.
2756 *--------------------------------------------------------------*/
2757
2758 if (CF_marker[i1] == 2)
2759 {
2760 CF_marker[i1] = 1;
2761 }
2762 }
2763 if ( num_procs > 1)
2764 {
2765 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
2766 {
2767 if(col_offd_S_to_A)
2768 i1 = col_offd_S_to_A[S_offd_j[jj]];
2769 else
2770 i1 = S_offd_j[jj];
2771 if ( CF_marker_offd[i1] == 2)
2772 {
2773 CF_marker_offd[i1] = 1;
2774 }
2775 }
2776 }
2777
2778
2779 jj_end_row = jj_counter;
2780 jj_end_row_offd = jj_counter_offd;
2781
2782 diagonal = A_diag_data[A_diag_i[i]];
2783 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
2784 { /* i1 is a c-point and strongly influences i, accumulate
2785 * a_(i,i1) into interpolation weight */
2786 i1 = A_diag_j[jj];
2787 if (P_marker[i1] >= jj_begin_row)
2788 {
2789 P_diag_data[P_marker[i1]] += A_diag_data[jj];
2790 }
2791 else if(P_marker[i1] == strong_f_marker)
2792 {
2793 sum = zero;
2794 sgn = 1;
2795 if(A_diag_data[A_diag_i[i1]] < 0) sgn = -1;
2796 for(jj1 = A_diag_i[i1]+1; jj1 < A_diag_i[i1+1]; jj1++)
2797 {
2798 i2 = A_diag_j[jj1];
2799 if((P_marker[i2] >= jj_begin_row || i2 == i) && (sgn*A_diag_data[jj1]) < 0)
2800 sum += A_diag_data[jj1];
2801 }
2802 if(num_procs > 1)
2803 {
2804 for(jj1 = A_offd_i[i1]; jj1< A_offd_i[i1+1]; jj1++)
2805 {
2806 i2 = A_offd_j[jj1];
2807 if(P_marker_offd[i2] >= jj_begin_row_offd &&
2808 (sgn*A_offd_data[jj1]) < 0)
2809 sum += A_offd_data[jj1];
2810 }
2811 }
2812 if(sum != 0)
2813 {
2814 distribute = A_diag_data[jj]/sum;
2815 /* Loop over row of A for point i1 and do the distribution */
2816 for(jj1 = A_diag_i[i1]; jj1 < A_diag_i[i1+1]; jj1++)
2817 {
2818 i2 = A_diag_j[jj1];
2819 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
2820 P_diag_data[P_marker[i2]] +=
2821 distribute*A_diag_data[jj1];
2822 if(i2 == i && (sgn*A_diag_data[jj1]) < 0)
2823 diagonal += distribute*A_diag_data[jj1];
2824 }
2825 if(num_procs > 1)
2826 {
2827 for(jj1 = A_offd_i[i1]; jj1 < A_offd_i[i1+1]; jj1++)
2828 {
2829 i2 = A_offd_j[jj1];
2830 if(P_marker_offd[i2] >= jj_begin_row_offd &&
2831 (sgn*A_offd_data[jj1]) < 0)
2832 P_offd_data[P_marker_offd[i2]] +=
2833 distribute*A_offd_data[jj1];
2834 }
2835 }
2836 }
2837 else
2838 diagonal += A_diag_data[jj];
2839 }
2840 /* neighbor i1 weakly influences i, accumulate a_(i,i1) into
2841 * diagonal */
2842 else if (CF_marker[i1] != -3)
2843 {
2844 if(num_functions == 1 || dof_func[i] == dof_func[i1])
2845 diagonal += A_diag_data[jj];
2846 }
2847 }
2848 if(num_procs > 1)
2849 {
2850 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
2851 {
2852 i1 = A_offd_j[jj];
2853 if(P_marker_offd[i1] >= jj_begin_row_offd)
2854 P_offd_data[P_marker_offd[i1]] += A_offd_data[jj];
2855 else if(P_marker_offd[i1] == strong_f_marker)
2856 {
2857 sum = zero;
2858 sgn = 1;
2859 if(A_ext_data[A_ext_i[i1]] < 0) sgn = -1;
2860 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
2861 {
2862 loc_col = A_ext_j[jj1];
2863 if(loc_col > -1)
2864 { /* diag */
2865 if((P_marker[loc_col] >= jj_begin_row || loc_col == i)
2866 && (sgn*A_ext_data[jj1]) < 0)
2867 sum += A_ext_data[jj1];
2868 }
2869 else
2870 {
2871 loc_col = -loc_col - 1;
2872 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
2873 (sgn*A_ext_data[jj1]) < 0)
2874 sum += A_ext_data[jj1];
2875 }
2876 }
2877 if(sum != 0)
2878 {
2879 distribute = A_offd_data[jj] / sum;
2880 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
2881 {
2882 loc_col = A_ext_j[jj1];
2883 if(loc_col > -1)
2884 { /* diag */
2885 if(P_marker[loc_col] >= jj_begin_row &&
2886 (sgn*A_ext_data[jj1]) < 0)
2887 P_diag_data[P_marker[loc_col]] += distribute*
2888 A_ext_data[jj1];
2889 if(loc_col == i && (sgn*A_ext_data[jj1]) < 0)
2890 diagonal += distribute*A_ext_data[jj1];
2891 }
2892 else
2893 {
2894 loc_col = -loc_col - 1;
2895 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
2896 (sgn*A_ext_data[jj1]) < 0)
2897 P_offd_data[P_marker_offd[loc_col]] += distribute*
2898 A_ext_data[jj1];
2899 }
2900 }
2901 }
2902 else
2903 diagonal += A_offd_data[jj];
2904 }
2905 else if (CF_marker_offd[i1] != -3)
2906 {
2907 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
2908 diagonal += A_offd_data[jj];
2909 }
2910 }
2911 }
2912 for(jj = jj_begin_row; jj < jj_end_row; jj++)
2913 P_diag_data[jj] /= -diagonal;
2914 for(jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
2915 P_offd_data[jj] /= -diagonal;
2916 }
2917 strong_f_marker--;
2918 }
2919
2920 P = hypre_ParCSRMatrixCreate(comm,
2921 hypre_ParCSRMatrixGlobalNumRows(A),
2922 total_global_cpts,
2923 hypre_ParCSRMatrixColStarts(A),
2924 num_cpts_global,
2925 0,
2926 P_diag_i[n_fine],
2927 P_offd_i[n_fine]);
2928
2929 P_diag = hypre_ParCSRMatrixDiag(P);
2930 hypre_CSRMatrixData(P_diag) = P_diag_data;
2931 hypre_CSRMatrixI(P_diag) = P_diag_i;
2932 hypre_CSRMatrixJ(P_diag) = P_diag_j;
2933 P_offd = hypre_ParCSRMatrixOffd(P);
2934 hypre_CSRMatrixData(P_offd) = P_offd_data;
2935 hypre_CSRMatrixI(P_offd) = P_offd_i;
2936 hypre_CSRMatrixJ(P_offd) = P_offd_j;
2937 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
2938
2939 /* Compress P, removing coefficients smaller than trunc_factor * Max */
2940 if (trunc_factor != 0.0 || max_elmts > 0)
2941 {
2942 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
2943 P_diag_data = hypre_CSRMatrixData(P_diag);
2944 P_diag_i = hypre_CSRMatrixI(P_diag);
2945 P_diag_j = hypre_CSRMatrixJ(P_diag);
2946 P_offd_data = hypre_CSRMatrixData(P_offd);
2947 P_offd_i = hypre_CSRMatrixI(P_offd);
2948 P_offd_j = hypre_CSRMatrixJ(P_offd);
2949 P_diag_size = P_diag_i[n_fine];
2950 P_offd_size = P_offd_i[n_fine];
2951 }
2952
2953 /* This builds col_map, col_map should be monotone increasing and contain
2954 * global numbers. */
2955 num_cols_P_offd = 0;
2956 if(P_offd_size)
2957 {
2958 num_cols_P_offd = 0;
2959 for (i=0; i < P_offd_size; i++)
2960 {
2961 index = P_offd_j[i];
2962 if(tmp_CF_marker_offd[index] == 1)
2963 {
2964 num_cols_P_offd++;
2965 tmp_CF_marker_offd[index] = 2;
2966 }
2967 }
2968
2969 if (num_cols_P_offd)
2970 {
2971 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
2972 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
2973 }
2974
2975 index = 0;
2976 for(i = 0; i < num_cols_P_offd; i++)
2977 {
2978 while( tmp_CF_marker_offd[index] != 2) index++;
2979 col_map_offd_P[i] = fine_to_coarse_offd[index];
2980 tmp_map_offd[i] = index++;
2981 }
2982
2983 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
2984
2985 for (i = 0; i < num_cols_P_offd; i++)
2986 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
2987
2988 for(i = 0; i < P_offd_size; i++)
2989 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
2990 }
2991
2992 if (num_cols_P_offd)
2993 {
2994 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
2995 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
2996 hypre_TFree(tmp_map_offd);
2997 }
2998
2999#ifdef HYPRE_NO_GLOBAL_PARTITION
3000 hypre_NewCommPkgCreate(P);
3001#else
3002 hypre_MatvecCommPkgCreate(P);
3003#endif
3004
3005 for (i=0; i < n_fine; i++)
3006 if (CF_marker[i] == -3) CF_marker[i] = -1;
3007
3008 *P_ptr = P;
3009
3010 /* Deallocate memory */
3011 hypre_TFree(fine_to_coarse);
3012 hypre_TFree(P_marker);
3013
3014 if (num_procs > 1)
3015 {
3016 /*hypre_TFree(clist_offd);*/
3017 hypre_BigCSRMatrixDestroy(Sop);
3018 hypre_BigCSRMatrixDestroy(A_ext);
3019 hypre_TFree(fine_to_coarse_offd);
3020 hypre_TFree(P_marker_offd);
3021 hypre_TFree(CF_marker_offd);
3022 hypre_TFree(tmp_CF_marker_offd);
3023 if(num_functions > 1)
3024 hypre_TFree(dof_func_offd);
3025 hypre_TFree(found);
3026
3027 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
3028 }
3029
3030 return hypre_error_flag;
3031}
3032
3033/*---------------------------------------------------------------------------
3034 * hypre_BoomerAMGBuildFFInterp
3035 * Comment: Only use FF when there is no common c point.
3036 *--------------------------------------------------------------------------*/
3037int
3038hypre_BoomerAMGBuildFFInterp(hypre_ParCSRMatrix *A, int *CF_marker,
3039 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
3040 int num_functions, int *dof_func, int debug_flag,
3041 double trunc_factor, int max_elmts,
3042 int *col_offd_S_to_A,
3043 hypre_ParCSRMatrix **P_ptr)
3044{
3045 /* Communication Variables */
3046 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
3047 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
3048
3049 int my_id, num_procs;
3050
3051 /* Variables to store input variables */
3052 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
3053 double *A_diag_data = hypre_CSRMatrixData(A_diag);
3054 int *A_diag_i = hypre_CSRMatrixI(A_diag);
3055 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
3056
3057 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
3058 double *A_offd_data = hypre_CSRMatrixData(A_offd);
3059 int *A_offd_i = hypre_CSRMatrixI(A_offd);
3060 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
3061
3062 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
3063 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
3064 int n_fine = hypre_CSRMatrixNumRows(A_diag);
3065 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
3066 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
3067 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt) local_numrows;
3068 HYPRE_BigInt total_global_cpts, my_first_cpt;
3069
3070 /* Variables to store strong connection matrix info */
3071 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
3072 int *S_diag_i = hypre_CSRMatrixI(S_diag);
3073 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
3074
3075 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
3076 int *S_offd_i = hypre_CSRMatrixI(S_offd);
3077 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
3078
3079 /* Interpolation matrix P */
3080 hypre_ParCSRMatrix *P;
3081 hypre_CSRMatrix *P_diag;
3082 hypre_CSRMatrix *P_offd;
3083
3084 double *P_diag_data = NULL;
3085 int *P_diag_i, *P_diag_j = NULL;
3086 double *P_offd_data = NULL;
3087 int *P_offd_i, *P_offd_j = NULL;
3088
3089 HYPRE_BigInt *col_map_offd_P;
3090 int *tmp_map_offd;
3091 int P_diag_size;
3092 int P_offd_size;
3093 int *P_marker = NULL;
3094 int *P_marker_offd = NULL;
3095 int *CF_marker_offd = NULL;
3096 int *tmp_CF_marker_offd = NULL;
3097 int *dof_func_offd = NULL;
3098 int ccounter_offd;
3099 int common_c;
3100
3101 /* Full row information for columns of A that are off diag*/
3102 hypre_BigCSRMatrix *A_ext;
3103 double *A_ext_data;
3104 int *A_ext_i;
3105 HYPRE_BigInt *big_A_ext_j;
3106 int *A_ext_j;
3107
3108 int *fine_to_coarse = NULL;
3109 HYPRE_BigInt *fine_to_coarse_offd = NULL;
3110 HYPRE_BigInt *found;
3111
3112 int num_cols_P_offd;
3113 int newoff, loc_col;
3114 int A_ext_rows, full_off_procNodes;
3115
3116 hypre_BigCSRMatrix *Sop;
3117 int *Sop_i;
3118 HYPRE_BigInt *big_Sop_j;
3119 int *Sop_j;
3120
3121 int Soprows;
3122
3123 /* Variables to keep count of interpolatory points */
3124 int jj_counter, jj_counter_offd;
3125 int jj_begin_row, jj_end_row;
3126 int jj_begin_row_offd = 0;
3127 int jj_end_row_offd = 0;
3128 int coarse_counter, coarse_counter_offd;
3129
3130 /* Interpolation weight variables */
3131 double sum, diagonal, distribute;
3132 int strong_f_marker = -2;
3133 int sgn;
3134
3135 /* Loop variables */
3136 int index;
3137 int start_indexing = 0;
3138 int i, i1, i2, jj, kk, k1, jj1;
3139 int ccounter;
3140
3141 /* Definitions */
3142 double zero = 0.0;
3143 double one = 1.0;
3144
3145 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
3146
3147 /* BEGIN */
3148 MPI_Comm_size(comm, &num_procs);
3149 MPI_Comm_rank(comm,&my_id);
3150
3151#ifdef HYPRE_NO_GLOBAL_PARTITION
3152 my_first_cpt = num_cpts_global[0];
3153 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
3154 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
3155#else
3156 my_first_cpt = num_cpts_global[my_id];
3157 total_global_cpts = num_cpts_global[num_procs];
3158#endif
3159
3160 if (!comm_pkg)
3161 {
3162#ifdef HYPRE_NO_GLOBAL_PARTITION
3163 hypre_NewCommPkgCreate(A);
3164#else
3165 hypre_MatvecCommPkgCreate(A);
3166#endif
3167 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
3168 }
3169
3170 /* Set up off processor information (specifically for neighbors of
3171 * neighbors */
3172 newoff = 0;
3173 full_off_procNodes = 0;
3174 if (num_procs > 1)
3175 {
3176 /*----------------------------------------------------------------------
3177 * Get the off processors rows for A and S, associated with columns in
3178 * A_offd and S_offd.
3179 *---------------------------------------------------------------------*/
3180 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
3181 A_ext_i = hypre_BigCSRMatrixI(A_ext);
3182 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
3183 A_ext_data = hypre_BigCSRMatrixData(A_ext);
3184 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
3185
3186 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
3187 Sop_i = hypre_BigCSRMatrixI(Sop);
3188 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
3189 Soprows = hypre_BigCSRMatrixNumRows(Sop);
3190
3191 /* Find nodes that are neighbors of neighbors, not found in offd */
3192 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
3193 &A_ext_j,
3194 Soprows, col_map_offd, col_1, col_n,
3195 Sop_i, big_Sop_j, &Sop_j, CF_marker, comm_pkg);
3196
3197 hypre_BigCSRMatrixJ(A_ext) = NULL;
3198 hypre_BigCSRMatrixJ(Sop) = NULL;
3199
3200 if(newoff >= 0)
3201 full_off_procNodes = newoff + num_cols_A_offd;
3202 else
3203 return(1);
3204
3205 /* Possibly add new points and new processors to the comm_pkg, all
3206 * processors need new_comm_pkg */
3207
3208 /* AHB - create a new comm package just for extended info -
3209 this will work better with the assumed partition*/
3210 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
3211 &extend_comm_pkg);
3212
3213 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
3214
3215 if (num_functions > 1 && full_off_procNodes > 0)
3216 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
3217
3218 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
3219 full_off_procNodes, CF_marker_offd);
3220
3221 if(num_functions > 1)
3222 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
3223 full_off_procNodes, dof_func_offd);
3224
3225
3226 }
3227
3228 /*-----------------------------------------------------------------------
3229 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
3230 *-----------------------------------------------------------------------*/
3231
3232 /*-----------------------------------------------------------------------
3233 * Intialize counters and allocate mapping vector.
3234 *-----------------------------------------------------------------------*/
3235 P_diag_i = hypre_CTAlloc(int, n_fine+1);
3236 P_offd_i = hypre_CTAlloc(int, n_fine+1);
3237
3238 if (n_fine)
3239 {
3240 fine_to_coarse = hypre_CTAlloc(int, n_fine);
3241 P_marker = hypre_CTAlloc(int, n_fine);
3242 }
3243
3244 if (full_off_procNodes)
3245 {
3246 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
3247 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
3248 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
3249 }
3250
3251 for (i=0; i < full_off_procNodes; i++)
3252 {
3253 P_marker_offd[i] = -1;
3254 tmp_CF_marker_offd[i] = -1;
3255 fine_to_coarse_offd[i] = -1;
3256 }
3257
3258 for (i=0; i < n_fine; i++)
3259 {
3260 P_marker[i] = -1;
3261 fine_to_coarse[i] = -1;
3262 }
3263
3264 jj_counter = start_indexing;
3265 jj_counter_offd = start_indexing;
3266 coarse_counter = 0;
3267 coarse_counter_offd = 0;
3268
3269 /*-----------------------------------------------------------------------
3270 * Loop over fine grid.
3271 *-----------------------------------------------------------------------*/
3272 for (i = 0; i < n_fine; i++)
3273 {
3274 P_diag_i[i] = jj_counter;
3275 if (num_procs > 1)
3276 P_offd_i[i] = jj_counter_offd;
3277
3278 if (CF_marker[i] > 0)
3279 {
3280 jj_counter++;
3281 fine_to_coarse[i] = coarse_counter;
3282 coarse_counter++;
3283 }
3284
3285 /*--------------------------------------------------------------------
3286 * If i is an F-point, interpolation is from the C-points that
3287 * strongly influence i, or C-points that stronly influence F-points
3288 * that strongly influence i.
3289 *--------------------------------------------------------------------*/
3290 else
3291 {
3292 /* Initialize ccounter for each f point */
3293 ccounter = 0;
3294 ccounter_offd = 0;
3295 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3296 { /* search through diag to find all c neighbors */
3297 i1 = S_diag_j[jj];
3298 if (CF_marker[i1] > 0)
3299 { /* i1 is a C point */
3300 CF_marker[i1] = 2;
3301 if (P_marker[i1] < P_diag_i[i])
3302 {
3303 P_marker[i1] = jj_counter;
3304 jj_counter++;
3305 }
3306 }
3307 }
3308 if(num_procs > 1)
3309 {
3310 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3311 { /* search through offd to find all c neighbors */
3312 if(col_offd_S_to_A)
3313 i1 = col_offd_S_to_A[S_offd_j[jj]];
3314 else
3315 i1 = S_offd_j[jj];
3316 if(CF_marker_offd[i1] > 0)
3317 { /* i1 is a C point direct neighbor */
3318 CF_marker_offd[i1] = 2;
3319 if(P_marker_offd[i1] < P_offd_i[i])
3320 {
3321 tmp_CF_marker_offd[i1] = 1;
3322 P_marker_offd[i1] = jj_counter_offd;
3323 jj_counter_offd++;
3324 }
3325 }
3326 }
3327 }
3328 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3329 { /* Search diag to find f neighbors and determine if common c point */
3330 i1 = S_diag_j[jj];
3331 if (CF_marker[i1] < 0)
3332 { /* i1 is a F point, loop through it's strong neighbors */
3333 common_c = 0;
3334 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
3335 {
3336 k1 = S_diag_j[kk];
3337 if (CF_marker[k1] == 2)
3338 {
3339 common_c = 1;
3340 break;
3341 }
3342 }
3343 if(num_procs > 1 && common_c == 0)
3344 { /* no common c point yet, check offd */
3345 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
3346 {
3347 if(col_offd_S_to_A)
3348 k1 = col_offd_S_to_A[S_offd_j[kk]];
3349 else
3350 k1 = S_offd_j[kk];
3351
3352 if (CF_marker_offd[k1] == 2)
3353 {
3354 common_c = 1;
3355 break;
3356 }
3357 }
3358 }
3359 if(!common_c)
3360 { /* No common c point, extend the interp set */
3361 for(kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
3362 {
3363 k1 = S_diag_j[kk];
3364 if(CF_marker[k1] > 0)
3365 {
3366 if(P_marker[k1] < P_diag_i[i])
3367 {
3368 P_marker[k1] = jj_counter;
3369 jj_counter++;
3370 }
3371 }
3372 }
3373 if(num_procs > 1)
3374 {
3375 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
3376 {
3377 if(col_offd_S_to_A)
3378 k1 = col_offd_S_to_A[S_offd_j[kk]];
3379 else
3380 k1 = S_offd_j[kk];
3381 if (CF_marker_offd[k1] > 0)
3382 {
3383 if(P_marker_offd[k1] < P_offd_i[i])
3384 {
3385 tmp_CF_marker_offd[k1] = 1;
3386 P_marker_offd[k1] = jj_counter_offd;
3387 jj_counter_offd++;
3388 }
3389 }
3390 }
3391 }
3392 }
3393 }
3394 }
3395 /* Look at off diag strong connections of i */
3396 if (num_procs > 1)
3397 {
3398 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3399 {
3400 i1 = S_offd_j[jj];
3401 if(col_offd_S_to_A)
3402 i1 = col_offd_S_to_A[i1];
3403 if (CF_marker_offd[i1] < 0)
3404 { /* F point; look at neighbors of i1. Sop contains global col
3405 * numbers and entries that could be in S_diag or S_offd or
3406 * neither. */
3407 common_c = 0;
3408 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
3409 { /* Check if common c */
3410 loc_col = Sop_j[kk];
3411 if(loc_col > -1)
3412 { /* In S_diag */
3413 if(CF_marker[loc_col] == 2)
3414 {
3415 common_c = 1;
3416 break;
3417 }
3418 }
3419 else
3420 {
3421 loc_col = -loc_col - 1;
3422 if(CF_marker_offd[loc_col] == 2)
3423 {
3424 common_c = 1;
3425 break;
3426 }
3427 }
3428 }
3429 if(!common_c)
3430 {
3431 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
3432 { /* Check if common c */
3433 loc_col = Sop_j[kk];
3434 if(loc_col > -1)
3435 { /* In S_diag */
3436 if(CF_marker[loc_col] > 0)
3437 {
3438 if(P_marker[loc_col] < P_diag_i[i])
3439 {
3440 P_marker[loc_col] = jj_counter;
3441 jj_counter++;
3442 }
3443 }
3444 }
3445 else
3446 {
3447 loc_col = -loc_col - 1;
3448 if(CF_marker_offd[loc_col] > 0)
3449 {
3450 if(P_marker_offd[loc_col] < P_offd_i[i])
3451 {
3452 P_marker_offd[loc_col] = jj_counter_offd;
3453 tmp_CF_marker_offd[loc_col] = 1;
3454 jj_counter_offd++;
3455 }
3456 }
3457 }
3458 }
3459 }
3460 }
3461 }
3462 }
3463 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3464 { /* search through diag to find all c neighbors */
3465 i1 = S_diag_j[jj];
3466 if (CF_marker[i1] == 2)
3467 CF_marker[i1] = 1;
3468 }
3469 if(num_procs > 1)
3470 {
3471 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3472 { /* search through offd to find all c neighbors */
3473 if(col_offd_S_to_A)
3474 i1 = col_offd_S_to_A[S_offd_j[jj]];
3475 else
3476 i1 = S_offd_j[jj];
3477 if(CF_marker_offd[i1] == 2)
3478 { /* i1 is a C point direct neighbor */
3479 CF_marker_offd[i1] = 1;
3480 }
3481 }
3482 }
3483 }
3484 }
3485
3486 /*-----------------------------------------------------------------------
3487 * Allocate arrays.
3488 *-----------------------------------------------------------------------*/
3489
3490 P_diag_size = jj_counter;
3491 P_offd_size = jj_counter_offd;
3492
3493 if (P_diag_size)
3494 {
3495 P_diag_j = hypre_CTAlloc(int, P_diag_size);
3496 P_diag_data = hypre_CTAlloc(double, P_diag_size);
3497 }
3498
3499 if (P_offd_size)
3500 {
3501 P_offd_j = hypre_CTAlloc(int, P_offd_size);
3502 P_offd_data = hypre_CTAlloc(double, P_offd_size);
3503 }
3504
3505 P_diag_i[n_fine] = jj_counter;
3506 P_offd_i[n_fine] = jj_counter_offd;
3507
3508 jj_counter = start_indexing;
3509 jj_counter_offd = start_indexing;
3510 ccounter = start_indexing;
3511 ccounter_offd = start_indexing;
3512
3513 /* Fine to coarse mapping */
3514 if(num_procs > 1)
3515 {
3516 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
3517 full_off_procNodes, my_first_cpt,
3518 fine_to_coarse_offd);
3519 }
3520
3521 for (i = 0; i < n_fine; i++)
3522 P_marker[i] = -1;
3523
3524 for (i = 0; i < full_off_procNodes; i++)
3525 P_marker_offd[i] = -1;
3526
3527 /*-----------------------------------------------------------------------
3528 * Loop over fine grid points.
3529 *-----------------------------------------------------------------------*/
3530 jj_begin_row_offd = 0;
3531 for (i = 0; i < n_fine; i++)
3532 {
3533 jj_begin_row = jj_counter;
3534 if(num_procs > 1)
3535 jj_begin_row_offd = jj_counter_offd;
3536
3537 /*--------------------------------------------------------------------
3538 * If i is a c-point, interpolation is the identity.
3539 *--------------------------------------------------------------------*/
3540
3541 if (CF_marker[i] > 0)
3542 {
3543 P_diag_j[jj_counter] = fine_to_coarse[i];
3544 P_diag_data[jj_counter] = one;
3545 jj_counter++;
3546 }
3547
3548 /*--------------------------------------------------------------------
3549 * If i is an F-point, build interpolation.
3550 *--------------------------------------------------------------------*/
3551
3552 else if (CF_marker[i] != -3)
3553 {
3554 ccounter = 0;
3555 ccounter_offd = 0;
3556 strong_f_marker--;
3557
3558 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3559 { /* Search C points only */
3560 i1 = S_diag_j[jj];
3561
3562 /*--------------------------------------------------------------
3563 * If neighbor i1 is a C-point, set column number in P_diag_j
3564 * and initialize interpolation weight to zero.
3565 *--------------------------------------------------------------*/
3566
3567 if (CF_marker[i1] > 0)
3568 {
3569 CF_marker[i1] = 2;
3570 if (P_marker[i1] < jj_begin_row)
3571 {
3572 P_marker[i1] = jj_counter;
3573 P_diag_j[jj_counter] = fine_to_coarse[i1];
3574 P_diag_data[jj_counter] = zero;
3575 jj_counter++;
3576 }
3577 }
3578 }
3579 if ( num_procs > 1)
3580 {
3581 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3582 {
3583 if(col_offd_S_to_A)
3584 i1 = col_offd_S_to_A[S_offd_j[jj]];
3585 else
3586 i1 = S_offd_j[jj];
3587 if ( CF_marker_offd[i1] > 0)
3588 {
3589 CF_marker_offd[i1] = 2;
3590 if(P_marker_offd[i1] < jj_begin_row_offd)
3591 {
3592 P_marker_offd[i1] = jj_counter_offd;
3593 P_offd_j[jj_counter_offd] = i1;
3594 P_offd_data[jj_counter_offd] = zero;
3595 jj_counter_offd++;
3596 }
3597 }
3598 }
3599 }
3600
3601 for(jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3602 { /* Search through F points */
3603 i1 = S_diag_j[jj];
3604 if(CF_marker[i1] == -1)
3605 {
3606 P_marker[i1] = strong_f_marker;
3607 common_c = 0;
3608 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
3609 {
3610 k1 = S_diag_j[kk];
3611 if (CF_marker[k1] == 2)
3612 {
3613 common_c = 1;
3614 break;
3615 }
3616 }
3617 if(num_procs > 1 && common_c == 0)
3618 { /* no common c point yet, check offd */
3619 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
3620 {
3621 if(col_offd_S_to_A)
3622 k1 = col_offd_S_to_A[S_offd_j[kk]];
3623 else
3624 k1 = S_offd_j[kk];
3625
3626 if (CF_marker_offd[k1] == 2)
3627 {
3628 common_c = 1;
3629 break;
3630 }
3631 }
3632 }
3633 if(!common_c)
3634 { /* No common c point, extend the interp set */
3635 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
3636 {
3637 k1 = S_diag_j[kk];
3638 if (CF_marker[k1] >= 0)
3639 {
3640 if(P_marker[k1] < jj_begin_row)
3641 {
3642 P_marker[k1] = jj_counter;
3643 P_diag_j[jj_counter] = fine_to_coarse[k1];
3644 P_diag_data[jj_counter] = zero;
3645 jj_counter++;
3646 }
3647 }
3648 }
3649 if(num_procs > 1)
3650 {
3651 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
3652 {
3653 if(col_offd_S_to_A)
3654 k1 = col_offd_S_to_A[S_offd_j[kk]];
3655 else
3656 k1 = S_offd_j[kk];
3657 if(CF_marker_offd[k1] >= 0)
3658 {
3659 if(P_marker_offd[k1] < jj_begin_row_offd)
3660 {
3661 P_marker_offd[k1] = jj_counter_offd;
3662 P_offd_j[jj_counter_offd] = k1;
3663 P_offd_data[jj_counter_offd] = zero;
3664 jj_counter_offd++;
3665 }
3666 }
3667 }
3668 }
3669 }
3670 }
3671 }
3672 if ( num_procs > 1)
3673 {
3674 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3675 {
3676 i1 = S_offd_j[jj];
3677 if(col_offd_S_to_A)
3678 i1 = col_offd_S_to_A[i1];
3679 if(CF_marker_offd[i1] == -1)
3680 { /* F points that are off proc */
3681 P_marker_offd[i1] = strong_f_marker;
3682 common_c = 0;
3683 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
3684 { /* Check if common c */
3685 loc_col = Sop_j[kk];
3686 if(loc_col > -1)
3687 { /* In S_diag */
3688 if(CF_marker[loc_col] == 2)
3689 {
3690 common_c = 1;
3691 break;
3692 }
3693 }
3694 else
3695 {
3696 loc_col = -loc_col - 1;
3697 if(CF_marker_offd[loc_col] == 2)
3698 {
3699 common_c = 1;
3700 break;
3701 }
3702 }
3703 }
3704 if(!common_c)
3705 {
3706 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
3707 {
3708 loc_col = Sop_j[kk];
3709 /* Find local col number */
3710 if(loc_col > -1)
3711 {
3712 if(CF_marker[loc_col] > 0)
3713 {
3714 if(P_marker[loc_col] < jj_begin_row)
3715 {
3716 P_marker[loc_col] = jj_counter;
3717 P_diag_j[jj_counter] = fine_to_coarse[loc_col];
3718 P_diag_data[jj_counter] = zero;
3719 jj_counter++;
3720 }
3721 }
3722 }
3723 else
3724 {
3725 loc_col = -loc_col - 1;
3726 if(CF_marker_offd[loc_col] > 0)
3727 {
3728 if(P_marker_offd[loc_col] < jj_begin_row_offd)
3729 {
3730 P_marker_offd[loc_col] = jj_counter_offd;
3731 P_offd_j[jj_counter_offd]=loc_col;
3732 P_offd_data[jj_counter_offd] = zero;
3733 jj_counter_offd++;
3734 }
3735 }
3736 }
3737 }
3738 }
3739 }
3740 }
3741 }
3742 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
3743 { /* Search C points only */
3744 i1 = S_diag_j[jj];
3745
3746 /*--------------------------------------------------------------
3747 * If neighbor i1 is a C-point, set column number in P_diag_j
3748 * and initialize interpolation weight to zero.
3749 *--------------------------------------------------------------*/
3750
3751 if (CF_marker[i1] == 2)
3752 {
3753 CF_marker[i1] = 1;
3754 }
3755 }
3756 if ( num_procs > 1)
3757 {
3758 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
3759 {
3760 if(col_offd_S_to_A)
3761 i1 = col_offd_S_to_A[S_offd_j[jj]];
3762 else
3763 i1 = S_offd_j[jj];
3764 if ( CF_marker_offd[i1] == 2)
3765 {
3766 CF_marker_offd[i1] = 1;
3767 }
3768 }
3769 }
3770
3771
3772 jj_end_row = jj_counter;
3773 jj_end_row_offd = jj_counter_offd;
3774
3775 diagonal = A_diag_data[A_diag_i[i]];
3776 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
3777 { /* i1 is a c-point and strongly influences i, accumulate
3778 * a_(i,i1) into interpolation weight */
3779 i1 = A_diag_j[jj];
3780 if (P_marker[i1] >= jj_begin_row)
3781 {
3782 P_diag_data[P_marker[i1]] += A_diag_data[jj];
3783 }
3784 else if(P_marker[i1] == strong_f_marker)
3785 {
3786 sum = zero;
3787 if(A_diag_data[A_diag_i[i1]] < 0) sgn = -1;
3788 /* Loop over row of A for point i1 and calculate the sum
3789 * of the connections to c-points that strongly incluence i. */
3790 for(jj1 = A_diag_i[i1]; jj1 < A_diag_i[i1+1]; jj1++)
3791 {
3792 i2 = A_diag_j[jj1];
3793 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
3794 sum += A_diag_data[jj1];
3795 }
3796 if(num_procs > 1)
3797 {
3798 for(jj1 = A_offd_i[i1]; jj1< A_offd_i[i1+1]; jj1++)
3799 {
3800 i2 = A_offd_j[jj1];
3801 if(P_marker_offd[i2] >= jj_begin_row_offd &&
3802 (sgn*A_offd_data[jj1]) < 0)
3803 sum += A_offd_data[jj1];
3804 }
3805 }
3806 if(sum != 0)
3807 {
3808 distribute = A_diag_data[jj]/sum;
3809 /* Loop over row of A for point i1 and do the distribution */
3810 for(jj1 = A_diag_i[i1]; jj1 < A_diag_i[i1+1]; jj1++)
3811 {
3812 i2 = A_diag_j[jj1];
3813 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
3814 P_diag_data[P_marker[i2]] +=
3815 distribute*A_diag_data[jj1];
3816 }
3817 if(num_procs > 1)
3818 {
3819 for(jj1 = A_offd_i[i1]; jj1 < A_offd_i[i1+1]; jj1++)
3820 {
3821 i2 = A_offd_j[jj1];
3822 if(P_marker_offd[i2] >= jj_begin_row_offd &&
3823 (sgn*A_offd_data[jj1]) < 0)
3824 P_offd_data[P_marker_offd[i2]] +=
3825 distribute*A_offd_data[jj1];
3826 }
3827 }
3828 }
3829 else
3830 diagonal += A_diag_data[jj];
3831 }
3832 /* neighbor i1 weakly influences i, accumulate a_(i,i1) into
3833 * diagonal */
3834 else if (CF_marker[i1] != -3)
3835 {
3836 if(num_functions == 1 || dof_func[i] == dof_func[i1])
3837 diagonal += A_diag_data[jj];
3838 }
3839 }
3840 if(num_procs > 1)
3841 {
3842 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
3843 {
3844 i1 = A_offd_j[jj];
3845 if(P_marker_offd[i1] >= jj_begin_row_offd)
3846 P_offd_data[P_marker_offd[i1]] += A_offd_data[jj];
3847 else if(P_marker_offd[i1] == strong_f_marker)
3848 {
3849 sum = zero;
3850 sgn = 1;
3851 if(A_ext_data[A_ext_i[i1]] < 0) sgn = -1;
3852 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
3853 {
3854 loc_col = A_ext_j[jj1];
3855 if(loc_col > -1)
3856 { /* diag */
3857 if(P_marker[loc_col] >= jj_begin_row && (sgn*A_ext_data[jj1])<0)
3858 sum += A_ext_data[jj1];
3859 }
3860 else
3861 {
3862 loc_col = -loc_col - 1;
3863 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
3864 (sgn*A_ext_data[jj1]) < 0)
3865 sum += A_ext_data[jj1];
3866 }
3867 }
3868 if(sum != 0)
3869 {
3870 distribute = A_offd_data[jj] / sum;
3871 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
3872 {
3873 loc_col = A_ext_j[jj1];
3874 if(loc_col > -1)
3875 { /* diag */
3876 if(P_marker[loc_col] >= jj_begin_row &&
3877 (sgn*A_ext_data[jj1]) < 0)
3878 P_diag_data[P_marker[loc_col]] += distribute*
3879 A_ext_data[jj1];
3880 }
3881 else
3882 {
3883 loc_col = -loc_col - 1;
3884 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
3885 (sgn*A_ext_data[jj1]) < 0)
3886 P_offd_data[P_marker_offd[loc_col]] += distribute*
3887 A_ext_data[jj1];
3888 }
3889 }
3890 }
3891 else
3892 diagonal += A_offd_data[jj];
3893 }
3894 else if (CF_marker_offd[i1] != -3)
3895 {
3896 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
3897 diagonal += A_offd_data[jj];
3898 }
3899 }
3900 }
3901 for(jj = jj_begin_row; jj < jj_end_row; jj++)
3902 P_diag_data[jj] /= -diagonal;
3903 for(jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
3904 P_offd_data[jj] /= -diagonal;
3905 }
3906 strong_f_marker--;
3907 }
3908
3909 P = hypre_ParCSRMatrixCreate(comm,
3910 hypre_ParCSRMatrixGlobalNumRows(A),
3911 total_global_cpts,
3912 hypre_ParCSRMatrixColStarts(A),
3913 num_cpts_global,
3914 0,
3915 P_diag_i[n_fine],
3916 P_offd_i[n_fine]);
3917
3918 P_diag = hypre_ParCSRMatrixDiag(P);
3919 hypre_CSRMatrixData(P_diag) = P_diag_data;
3920 hypre_CSRMatrixI(P_diag) = P_diag_i;
3921 hypre_CSRMatrixJ(P_diag) = P_diag_j;
3922 P_offd = hypre_ParCSRMatrixOffd(P);
3923 hypre_CSRMatrixData(P_offd) = P_offd_data;
3924 hypre_CSRMatrixI(P_offd) = P_offd_i;
3925 hypre_CSRMatrixJ(P_offd) = P_offd_j;
3926 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
3927
3928 /* Compress P, removing coefficients smaller than trunc_factor * Max */
3929 if (trunc_factor != 0.0 || max_elmts > 0)
3930 {
3931 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
3932 P_diag_data = hypre_CSRMatrixData(P_diag);
3933 P_diag_i = hypre_CSRMatrixI(P_diag);
3934 P_diag_j = hypre_CSRMatrixJ(P_diag);
3935 P_offd_data = hypre_CSRMatrixData(P_offd);
3936 P_offd_i = hypre_CSRMatrixI(P_offd);
3937 P_offd_j = hypre_CSRMatrixJ(P_offd);
3938 P_diag_size = P_diag_i[n_fine];
3939 P_offd_size = P_offd_i[n_fine];
3940 }
3941
3942 /* This builds col_map, col_map should be monotone increasing and contain
3943 * global numbers. */
3944 num_cols_P_offd = 0;
3945 if(P_offd_size)
3946 {
3947 num_cols_P_offd = 0;
3948 for (i=0; i < P_offd_size; i++)
3949 {
3950 index = P_offd_j[i];
3951 if(tmp_CF_marker_offd[index] == 1)
3952 {
3953 num_cols_P_offd++;
3954 tmp_CF_marker_offd[index] = 2;
3955 }
3956 }
3957
3958 if (num_cols_P_offd)
3959 {
3960 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
3961 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
3962 }
3963
3964 index = 0;
3965 for(i = 0; i < num_cols_P_offd; i++)
3966 {
3967 while( tmp_CF_marker_offd[index] == -1) index++;
3968 col_map_offd_P[i] = fine_to_coarse_offd[index];
3969 tmp_map_offd[i] = index++;
3970 }
3971
3972 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
3973
3974 for (i = 0; i < num_cols_P_offd; i++)
3975 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
3976
3977 for(i = 0; i < P_offd_size; i++)
3978 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
3979 }
3980
3981 if (num_cols_P_offd)
3982 {
3983 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
3984 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
3985 hypre_TFree(tmp_map_offd);
3986 }
3987
3988#ifdef HYPRE_NO_GLOBAL_PARTITION
3989 hypre_NewCommPkgCreate(P);
3990#else
3991 hypre_MatvecCommPkgCreate(P);
3992#endif
3993
3994 for (i=0; i < n_fine; i++)
3995 if (CF_marker[i] == -3) CF_marker[i] = -1;
3996
3997 *P_ptr = P;
3998
3999 /* Deallocate memory */
4000 hypre_TFree(fine_to_coarse);
4001 hypre_TFree(P_marker);
4002
4003 if (num_procs > 1)
4004 {
4005 hypre_BigCSRMatrixDestroy(Sop);
4006 hypre_BigCSRMatrixDestroy(A_ext);
4007 hypre_TFree(fine_to_coarse_offd);
4008 hypre_TFree(P_marker_offd);
4009 hypre_TFree(CF_marker_offd);
4010 hypre_TFree(tmp_CF_marker_offd);
4011 if(num_functions > 1)
4012 hypre_TFree(dof_func_offd);
4013 hypre_TFree(found);
4014
4015 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
4016
4017 }
4018
4019 return hypre_error_flag;
4020}
4021/*---------------------------------------------------------------------------
4022 * hypre_BoomerAMGBuildFF1Interp
4023 * Comment: Only use FF when there is no common c point.
4024 *--------------------------------------------------------------------------*/
4025int
4026hypre_BoomerAMGBuildFF1Interp(hypre_ParCSRMatrix *A, int *CF_marker,
4027 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
4028 int num_functions, int *dof_func, int debug_flag,
4029 double trunc_factor, int max_elmts,
4030 int *col_offd_S_to_A,
4031 hypre_ParCSRMatrix **P_ptr)
4032{
4033 /* Communication Variables */
4034 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
4035 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
4036
4037 int my_id, num_procs;
4038
4039 /* Variables to store input variables */
4040 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
4041 double *A_diag_data = hypre_CSRMatrixData(A_diag);
4042 int *A_diag_i = hypre_CSRMatrixI(A_diag);
4043 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
4044
4045 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
4046 double *A_offd_data = hypre_CSRMatrixData(A_offd);
4047 int *A_offd_i = hypre_CSRMatrixI(A_offd);
4048 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
4049
4050 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
4051 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
4052 int n_fine = hypre_CSRMatrixNumRows(A_diag);
4053 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
4054 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
4055 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt) local_numrows;
4056 HYPRE_BigInt total_global_cpts, my_first_cpt;
4057
4058 /* Variables to store strong connection matrix info */
4059 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
4060 int *S_diag_i = hypre_CSRMatrixI(S_diag);
4061 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
4062
4063 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
4064 int *S_offd_i = hypre_CSRMatrixI(S_offd);
4065 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
4066
4067 /* Interpolation matrix P */
4068 hypre_ParCSRMatrix *P;
4069 hypre_CSRMatrix *P_diag;
4070 hypre_CSRMatrix *P_offd;
4071
4072 double *P_diag_data = NULL;
4073 int *P_diag_i, *P_diag_j = NULL;
4074 double *P_offd_data = NULL;
4075 int *P_offd_i, *P_offd_j = NULL;
4076
4077 HYPRE_BigInt *col_map_offd_P;
4078 int *tmp_map_offd;
4079 int P_diag_size;
4080 int P_offd_size;
4081 int *P_marker = NULL;
4082 int *P_marker_offd = NULL;
4083 int *CF_marker_offd = NULL;
4084 int *tmp_CF_marker_offd = NULL;
4085 int *dof_func_offd = NULL;
4086 int ccounter_offd;
4087 int common_c;
4088
4089 /* Full row information for columns of A that are off diag*/
4090 hypre_BigCSRMatrix *A_ext;
4091 double *A_ext_data;
4092 int *A_ext_i;
4093 HYPRE_BigInt *big_A_ext_j;
4094 int *A_ext_j;
4095
4096 int *fine_to_coarse = NULL;
4097 HYPRE_BigInt *fine_to_coarse_offd = NULL;
4098 HYPRE_BigInt *found;
4099
4100 int num_cols_P_offd;
4101 int newoff, loc_col;
4102 int A_ext_rows, full_off_procNodes;
4103
4104 hypre_BigCSRMatrix *Sop;
4105 int *Sop_i;
4106 HYPRE_BigInt *big_Sop_j;
4107 int *Sop_j;
4108
4109 int Soprows;
4110
4111 /* Variables to keep count of interpolatory points */
4112 int jj_counter, jj_counter_offd;
4113 int jj_begin_row, jj_end_row;
4114 int jj_begin_row_offd = 0;
4115 int jj_end_row_offd = 0;
4116 int coarse_counter, coarse_counter_offd;
4117
4118 /* Interpolation weight variables */
4119 double sum, diagonal, distribute;
4120 int strong_f_marker = -2;
4121 int sgn;
4122
4123 /* Loop variables */
4124 int index;
4125 int start_indexing = 0;
4126 int i, i1, i2, jj, kk, k1, jj1;
4127 int ccounter;
4128 int found_c = 0;
4129
4130 /* Definitions */
4131 double zero = 0.0;
4132 double one = 1.0;
4133
4134 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
4135 /* BEGIN */
4136 MPI_Comm_size(comm, &num_procs);
4137 MPI_Comm_rank(comm,&my_id);
4138
4139#ifdef HYPRE_NO_GLOBAL_PARTITION
4140 my_first_cpt = num_cpts_global[0];
4141 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
4142 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
4143#else
4144 my_first_cpt = num_cpts_global[my_id];
4145 total_global_cpts = num_cpts_global[num_procs];
4146#endif
4147
4148 if (!comm_pkg)
4149 {
4150#ifdef HYPRE_NO_GLOBAL_PARTITION
4151 hypre_NewCommPkgCreate(A);
4152#else
4153 hypre_MatvecCommPkgCreate(A);
4154#endif
4155 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
4156 }
4157
4158 /* Set up off processor information (specifically for neighbors of
4159 * neighbors */
4160 newoff = 0;
4161 full_off_procNodes = 0;
4162 if (num_procs > 1)
4163 {
4164 /*----------------------------------------------------------------------
4165 * Get the off processors rows for A and S, associated with columns in
4166 * A_offd and S_offd.
4167 *---------------------------------------------------------------------*/
4168 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
4169 A_ext_i = hypre_BigCSRMatrixI(A_ext);
4170 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
4171 A_ext_data = hypre_BigCSRMatrixData(A_ext);
4172 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
4173
4174 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
4175 Sop_i = hypre_BigCSRMatrixI(Sop);
4176 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
4177 Soprows = hypre_BigCSRMatrixNumRows(Sop);
4178
4179 /* Find nodes that are neighbors of neighbors, not found in offd */
4180 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
4181 &A_ext_j,
4182 Soprows, col_map_offd, col_1, col_n,
4183 Sop_i, big_Sop_j, &Sop_j, CF_marker, comm_pkg);
4184
4185 hypre_BigCSRMatrixJ(A_ext) = NULL;
4186 hypre_BigCSRMatrixJ(Sop) = NULL;
4187
4188 if(newoff >= 0)
4189 full_off_procNodes = newoff + num_cols_A_offd;
4190 else
4191 return(1);
4192
4193 /* Possibly add new points and new processors to the comm_pkg, all
4194 * processors need new_comm_pkg */
4195
4196 /* AHB - create a new comm package just for extended info -
4197 this will work better with the assumed partition*/
4198 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
4199 &extend_comm_pkg);
4200
4201 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
4202
4203 if (num_functions > 1 && full_off_procNodes > 0)
4204 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
4205
4206 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
4207 full_off_procNodes, CF_marker_offd);
4208
4209 if(num_functions > 1)
4210 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
4211 full_off_procNodes, dof_func_offd);
4212 }
4213
4214 /*-----------------------------------------------------------------------
4215 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
4216 *-----------------------------------------------------------------------*/
4217
4218 /*-----------------------------------------------------------------------
4219 * Intialize counters and allocate mapping vector.
4220 *-----------------------------------------------------------------------*/
4221 P_diag_i = hypre_CTAlloc(int, n_fine+1);
4222 P_offd_i = hypre_CTAlloc(int, n_fine+1);
4223
4224 if (n_fine)
4225 {
4226 fine_to_coarse = hypre_CTAlloc(int, n_fine);
4227 P_marker = hypre_CTAlloc(int, n_fine);
4228 }
4229
4230 if (full_off_procNodes)
4231 {
4232 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
4233 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
4234 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
4235 }
4236
4237 for (i=0; i < full_off_procNodes; i++)
4238 {
4239 P_marker_offd[i] = -1;
4240 tmp_CF_marker_offd[i] = -1;
4241 fine_to_coarse_offd[i] = -1;
4242 }
4243
4244 for (i=0; i < n_fine; i++)
4245 {
4246 P_marker[i] = -1;
4247 fine_to_coarse[i] = -1;
4248 }
4249
4250 jj_counter = start_indexing;
4251 jj_counter_offd = start_indexing;
4252 coarse_counter = 0;
4253 coarse_counter_offd = 0;
4254
4255 /*-----------------------------------------------------------------------
4256 * Loop over fine grid.
4257 *-----------------------------------------------------------------------*/
4258 for (i = 0; i < n_fine; i++)
4259 {
4260 P_diag_i[i] = jj_counter;
4261 if (num_procs > 1)
4262 P_offd_i[i] = jj_counter_offd;
4263
4264 if (CF_marker[i] >= 0)
4265 {
4266 jj_counter++;
4267 fine_to_coarse[i] = coarse_counter;
4268 coarse_counter++;
4269 }
4270
4271 /*--------------------------------------------------------------------
4272 * If i is an F-point, interpolation is from the C-points that
4273 * strongly influence i, or C-points that stronly influence F-points
4274 * that strongly influence i.
4275 *--------------------------------------------------------------------*/
4276 else
4277 {
4278 /* Initialize ccounter for each f point */
4279 ccounter = 0;
4280 ccounter_offd = 0;
4281 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4282 { /* search through diag to find all c neighbors */
4283 i1 = S_diag_j[jj];
4284 if (CF_marker[i1] > 0)
4285 { /* i1 is a C point */
4286 CF_marker[i1] = 2;
4287 if (P_marker[i1] < P_diag_i[i])
4288 {
4289 P_marker[i1] = jj_counter;
4290 jj_counter++;
4291 }
4292 }
4293 }
4294 if(num_procs > 1)
4295 {
4296 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4297 { /* search through offd to find all c neighbors */
4298 if(col_offd_S_to_A)
4299 i1 = col_offd_S_to_A[S_offd_j[jj]];
4300 else
4301 i1 = S_offd_j[jj];
4302 if(CF_marker_offd[i1] > 0)
4303 { /* i1 is a C point direct neighbor */
4304 CF_marker_offd[i1] = 2;
4305 if(P_marker_offd[i1] < P_offd_i[i])
4306 {
4307 tmp_CF_marker_offd[i1] = 1;
4308 P_marker_offd[i1] = jj_counter_offd;
4309 jj_counter_offd++;
4310 }
4311 }
4312 }
4313 }
4314 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4315 { /* Search diag to find f neighbors and determine if common c point */
4316 i1 = S_diag_j[jj];
4317 if (CF_marker[i1] < 0)
4318 { /* i1 is a F point, loop through it's strong neighbors */
4319 common_c = 0;
4320 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
4321 {
4322 k1 = S_diag_j[kk];
4323 if (CF_marker[k1] == 2)
4324 {
4325 common_c = 1;
4326 break;
4327 }
4328 }
4329 if(num_procs > 1 && common_c == 0)
4330 { /* no common c point yet, check offd */
4331 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
4332 {
4333 if(col_offd_S_to_A)
4334 k1 = col_offd_S_to_A[S_offd_j[kk]];
4335 else
4336 k1 = S_offd_j[kk];
4337
4338 if (CF_marker_offd[k1] == 2)
4339 { /* k1 is a c point check if it is common */
4340 common_c = 1;
4341 break;
4342 }
4343 }
4344 }
4345 if(!common_c)
4346 { /* No common c point, extend the interp set */
4347 found_c = 0;
4348 for(kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
4349 {
4350 k1 = S_diag_j[kk];
4351 if(CF_marker[k1] > 0)
4352 {
4353 if(P_marker[k1] < P_diag_i[i])
4354 {
4355 P_marker[k1] = jj_counter;
4356 jj_counter++;
4357 found_c = 1;
4358 break;
4359 }
4360 }
4361 }
4362 if(num_procs > 1 && !found_c)
4363 {
4364 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
4365 {
4366 if(col_offd_S_to_A)
4367 k1 = col_offd_S_to_A[S_offd_j[kk]];
4368 else
4369 k1 = S_offd_j[kk];
4370 if (CF_marker_offd[k1] > 0)
4371 {
4372 if(P_marker_offd[k1] < P_offd_i[i])
4373 {
4374 tmp_CF_marker_offd[k1] = 1;
4375 P_marker_offd[k1] = jj_counter_offd;
4376 jj_counter_offd++;
4377 break;
4378 }
4379 }
4380 }
4381 }
4382 }
4383 }
4384 }
4385 /* Look at off diag strong connections of i */
4386 if (num_procs > 1)
4387 {
4388 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4389 {
4390 i1 = S_offd_j[jj];
4391 if(col_offd_S_to_A)
4392 i1 = col_offd_S_to_A[i1];
4393 if (CF_marker_offd[i1] < 0)
4394 { /* F point; look at neighbors of i1. Sop contains global col
4395 * numbers and entries that could be in S_diag or S_offd or
4396 * neither. */
4397 common_c = 0;
4398 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
4399 { /* Check if common c */
4400 loc_col = Sop_j[kk];
4401 if(loc_col > -1)
4402 { /* In S_diag */
4403 if(CF_marker[loc_col] == 2)
4404 {
4405 common_c = 1;
4406 break;
4407 }
4408 }
4409 else
4410 {
4411 loc_col = -loc_col - 1;
4412 if(CF_marker_offd[loc_col] == 2)
4413 {
4414 common_c = 1;
4415 break;
4416 }
4417 }
4418 }
4419 if(!common_c)
4420 {
4421 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
4422 { /* Check if common c */
4423 loc_col = Sop_j[kk];
4424 if(loc_col > -1)
4425 { /* In S_diag */
4426 if(CF_marker[loc_col] > 0)
4427 {
4428 if(P_marker[loc_col] < P_diag_i[i])
4429 {
4430 P_marker[loc_col] = jj_counter;
4431 jj_counter++;
4432 break;
4433 }
4434 }
4435 }
4436 else
4437 {
4438 loc_col = -loc_col - 1;
4439 if(CF_marker_offd[loc_col] > 0)
4440 {
4441 if(P_marker_offd[loc_col] < P_offd_i[i])
4442 {
4443 P_marker_offd[loc_col] = jj_counter_offd;
4444 tmp_CF_marker_offd[loc_col] = 1;
4445 jj_counter_offd++;
4446 break;
4447 }
4448 }
4449 }
4450 }
4451 }
4452 }
4453 }
4454 }
4455 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4456 { /* search through diag to find all c neighbors */
4457 i1 = S_diag_j[jj];
4458 if (CF_marker[i1] == 2)
4459 CF_marker[i1] = 1;
4460 }
4461 if(num_procs > 1)
4462 {
4463 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4464 { /* search through offd to find all c neighbors */
4465 if(col_offd_S_to_A)
4466 i1 = col_offd_S_to_A[S_offd_j[jj]];
4467 else
4468 i1 = S_offd_j[jj];
4469 if(CF_marker_offd[i1] == 2)
4470 { /* i1 is a C point direct neighbor */
4471 CF_marker_offd[i1] = 1;
4472 }
4473 }
4474 }
4475 }
4476 }
4477
4478 /*-----------------------------------------------------------------------
4479 * Allocate arrays.
4480 *-----------------------------------------------------------------------*/
4481
4482 P_diag_size = jj_counter;
4483 P_offd_size = jj_counter_offd;
4484 if (P_diag_size)
4485 {
4486 P_diag_j = hypre_CTAlloc(int, P_diag_size);
4487 P_diag_data = hypre_CTAlloc(double, P_diag_size);
4488 }
4489
4490 if (P_offd_size)
4491 {
4492 P_offd_j = hypre_CTAlloc(int, P_offd_size);
4493 P_offd_data = hypre_CTAlloc(double, P_offd_size);
4494 }
4495
4496 P_diag_i[n_fine] = jj_counter;
4497 P_offd_i[n_fine] = jj_counter_offd;
4498
4499 jj_counter = start_indexing;
4500 jj_counter_offd = start_indexing;
4501 ccounter = start_indexing;
4502 ccounter_offd = start_indexing;
4503
4504 /* Fine to coarse mapping */
4505 if(num_procs > 1)
4506 {
4507 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
4508 full_off_procNodes, my_first_cpt,
4509 fine_to_coarse_offd);
4510 }
4511
4512 for (i = 0; i < n_fine; i++)
4513 P_marker[i] = -1;
4514
4515 for (i = 0; i < full_off_procNodes; i++)
4516 P_marker_offd[i] = -1;
4517
4518 /*-----------------------------------------------------------------------
4519 * Loop over fine grid points.
4520 *-----------------------------------------------------------------------*/
4521 jj_begin_row_offd = 0;
4522 for (i = 0; i < n_fine; i++)
4523 {
4524 jj_begin_row = jj_counter;
4525 if(num_procs > 1)
4526 jj_begin_row_offd = jj_counter_offd;
4527
4528 /*--------------------------------------------------------------------
4529 * If i is a c-point, interpolation is the identity.
4530 *--------------------------------------------------------------------*/
4531
4532 if (CF_marker[i] > 0)
4533 {
4534 P_diag_j[jj_counter] = fine_to_coarse[i];
4535 P_diag_data[jj_counter] = one;
4536 jj_counter++;
4537 }
4538
4539 /*--------------------------------------------------------------------
4540 * If i is an F-point, build interpolation.
4541 *--------------------------------------------------------------------*/
4542
4543 else if (CF_marker[i] != -3)
4544 {
4545 ccounter = 0;
4546 ccounter_offd = 0;
4547 strong_f_marker--;
4548
4549 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4550 { /* Search C points only */
4551 i1 = S_diag_j[jj];
4552
4553 /*--------------------------------------------------------------
4554 * If neighbor i1 is a C-point, set column number in P_diag_j
4555 * and initialize interpolation weight to zero.
4556 *--------------------------------------------------------------*/
4557
4558 if (CF_marker[i1] > 0)
4559 {
4560 CF_marker[i1] = 2;
4561 if (P_marker[i1] < jj_begin_row)
4562 {
4563 P_marker[i1] = jj_counter;
4564 P_diag_j[jj_counter] = fine_to_coarse[i1];
4565 P_diag_data[jj_counter] = zero;
4566 jj_counter++;
4567 }
4568 }
4569 }
4570 if ( num_procs > 1)
4571 {
4572 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4573 {
4574 if(col_offd_S_to_A)
4575 i1 = col_offd_S_to_A[S_offd_j[jj]];
4576 else
4577 i1 = S_offd_j[jj];
4578 if ( CF_marker_offd[i1] > 0)
4579 {
4580 CF_marker_offd[i1] = 2;
4581 if(P_marker_offd[i1] < jj_begin_row_offd)
4582 {
4583 P_marker_offd[i1] = jj_counter_offd;
4584 P_offd_j[jj_counter_offd] = i1;
4585 P_offd_data[jj_counter_offd] = zero;
4586 jj_counter_offd++;
4587 }
4588 }
4589 }
4590 }
4591
4592 for(jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4593 { /* Search through F points */
4594 i1 = S_diag_j[jj];
4595 if(CF_marker[i1] == -1)
4596 {
4597 P_marker[i1] = strong_f_marker;
4598 common_c = 0;
4599 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
4600 {
4601 k1 = S_diag_j[kk];
4602 if (CF_marker[k1] == 2)
4603 {
4604 common_c = 1;
4605 break;
4606 }
4607 }
4608 if(num_procs > 1 && common_c == 0)
4609 { /* no common c point yet, check offd */
4610 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
4611 {
4612 if(col_offd_S_to_A)
4613 k1 = col_offd_S_to_A[S_offd_j[kk]];
4614 else
4615 k1 = S_offd_j[kk];
4616
4617 if (CF_marker_offd[k1] == 2)
4618 { /* k1 is a c point check if it is common */
4619 common_c = 1;
4620 break;
4621 }
4622 }
4623 }
4624 if(!common_c)
4625 { /* No common c point, extend the interp set */
4626 found_c = 0;
4627 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
4628 {
4629 k1 = S_diag_j[kk];
4630 if (CF_marker[k1] >= 0)
4631 {
4632 if(P_marker[k1] < jj_begin_row)
4633 {
4634 P_marker[k1] = jj_counter;
4635 P_diag_j[jj_counter] = fine_to_coarse[k1];
4636 P_diag_data[jj_counter] = zero;
4637 jj_counter++;
4638 found_c = 1;
4639 break;
4640 }
4641 }
4642 }
4643 if(num_procs > 1 && !found_c)
4644 {
4645 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
4646 {
4647 if(col_offd_S_to_A)
4648 k1 = col_offd_S_to_A[S_offd_j[kk]];
4649 else
4650 k1 = S_offd_j[kk];
4651 if(CF_marker_offd[k1] >= 0)
4652 {
4653 if(P_marker_offd[k1] < jj_begin_row_offd)
4654 {
4655 P_marker_offd[k1] = jj_counter_offd;
4656 P_offd_j[jj_counter_offd] = k1;
4657 P_offd_data[jj_counter_offd] = zero;
4658 jj_counter_offd++;
4659 break;
4660 }
4661 }
4662 }
4663 }
4664 }
4665 }
4666 }
4667 if ( num_procs > 1)
4668 {
4669 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4670 {
4671 i1 = S_offd_j[jj];
4672 if(col_offd_S_to_A)
4673 i1 = col_offd_S_to_A[i1];
4674 if(CF_marker_offd[i1] == -1)
4675 { /* F points that are off proc */
4676 P_marker_offd[i1] = strong_f_marker;
4677 common_c = 0;
4678 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
4679 { /* Check if common c */
4680 loc_col = Sop_j[kk];
4681 if(loc_col > -1)
4682 { /* In S_diag */
4683 if(CF_marker[loc_col] == 2)
4684 {
4685 common_c = 1;
4686 break;
4687 }
4688 }
4689 else
4690 {
4691 loc_col = -loc_col - 1;
4692 if(CF_marker_offd[loc_col] == 2)
4693 {
4694 common_c = 1;
4695 break;
4696 }
4697 }
4698 }
4699 if(!common_c)
4700 {
4701 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
4702 {
4703 loc_col = Sop_j[kk];
4704 /* Find local col number */
4705 if(loc_col > -1)
4706 {
4707 if(CF_marker[loc_col] > 0)
4708 {
4709 if(P_marker[loc_col] < jj_begin_row)
4710 {
4711 P_marker[loc_col] = jj_counter;
4712 P_diag_j[jj_counter] = fine_to_coarse[loc_col];
4713 P_diag_data[jj_counter] = zero;
4714 jj_counter++;
4715 break;
4716 }
4717 }
4718 }
4719 else
4720 {
4721 loc_col = -loc_col - 1;
4722 if(CF_marker_offd[loc_col] > 0)
4723 {
4724 if(P_marker_offd[loc_col] < jj_begin_row_offd)
4725 {
4726 P_marker_offd[loc_col] = jj_counter_offd;
4727 P_offd_j[jj_counter_offd]=loc_col;
4728 P_offd_data[jj_counter_offd] = zero;
4729 jj_counter_offd++;
4730 break;
4731 }
4732 }
4733 }
4734 }
4735 }
4736 }
4737 }
4738 }
4739 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
4740 { /* Search C points only */
4741 i1 = S_diag_j[jj];
4742
4743 /*--------------------------------------------------------------
4744 * If neighbor i1 is a C-point, set column number in P_diag_j
4745 * and initialize interpolation weight to zero.
4746 *--------------------------------------------------------------*/
4747
4748 if (CF_marker[i1] == 2)
4749 {
4750 CF_marker[i1] = 1;
4751 }
4752 }
4753 if ( num_procs > 1)
4754 {
4755 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
4756 {
4757 if(col_offd_S_to_A)
4758 i1 = col_offd_S_to_A[S_offd_j[jj]];
4759 else
4760 i1 = S_offd_j[jj];
4761 if ( CF_marker_offd[i1] == 2)
4762 {
4763 CF_marker_offd[i1] = 1;
4764 }
4765 }
4766 }
4767
4768
4769 jj_end_row = jj_counter;
4770 jj_end_row_offd = jj_counter_offd;
4771
4772 diagonal = A_diag_data[A_diag_i[i]];
4773 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
4774 { /* i1 is a c-point and strongly influences i, accumulate
4775 * a_(i,i1) into interpolation weight */
4776 i1 = A_diag_j[jj];
4777 if (P_marker[i1] >= jj_begin_row)
4778 {
4779 P_diag_data[P_marker[i1]] += A_diag_data[jj];
4780 }
4781 else if(P_marker[i1] == strong_f_marker)
4782 {
4783 sum = zero;
4784 if(A_diag_data[A_diag_i[i1]] < 0) sgn = -1;
4785 /* Loop over row of A for point i1 and calculate the sum
4786 * of the connections to c-points that strongly incluence i. */
4787 for(jj1 = A_diag_i[i1]; jj1 < A_diag_i[i1+1]; jj1++)
4788 {
4789 i2 = A_diag_j[jj1];
4790 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
4791 sum += A_diag_data[jj1];
4792 }
4793 if(num_procs > 1)
4794 {
4795 for(jj1 = A_offd_i[i1]; jj1< A_offd_i[i1+1]; jj1++)
4796 {
4797 i2 = A_offd_j[jj1];
4798 if(P_marker_offd[i2] >= jj_begin_row_offd &&
4799 (sgn*A_offd_data[jj1]) < 0)
4800 sum += A_offd_data[jj1];
4801 }
4802 }
4803 if(sum != 0)
4804 {
4805 distribute = A_diag_data[jj]/sum;
4806 /* Loop over row of A for point i1 and do the distribution */
4807 for(jj1 = A_diag_i[i1]; jj1 < A_diag_i[i1+1]; jj1++)
4808 {
4809 i2 = A_diag_j[jj1];
4810 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
4811 P_diag_data[P_marker[i2]] +=
4812 distribute*A_diag_data[jj1];
4813 }
4814 if(num_procs > 1)
4815 {
4816 for(jj1 = A_offd_i[i1]; jj1 < A_offd_i[i1+1]; jj1++)
4817 {
4818 i2 = A_offd_j[jj1];
4819 if(P_marker_offd[i2] >= jj_begin_row_offd &&
4820 (sgn*A_offd_data[jj1]) < 0)
4821 P_offd_data[P_marker_offd[i2]] +=
4822 distribute*A_offd_data[jj1];
4823 }
4824 }
4825 }
4826 else
4827 diagonal += A_diag_data[jj];
4828 }
4829 /* neighbor i1 weakly influences i, accumulate a_(i,i1) into
4830 * diagonal */
4831 else if (CF_marker[i1] != -3)
4832 {
4833 if(num_functions == 1 || dof_func[i] == dof_func[i1])
4834 diagonal += A_diag_data[jj];
4835 }
4836 }
4837 if(num_procs > 1)
4838 {
4839 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
4840 {
4841 i1 = A_offd_j[jj];
4842 if(P_marker_offd[i1] >= jj_begin_row_offd)
4843 P_offd_data[P_marker_offd[i1]] += A_offd_data[jj];
4844 else if(P_marker_offd[i1] == strong_f_marker)
4845 {
4846 sum = zero;
4847 sgn = 1;
4848 if(A_ext_data[A_ext_i[i1]] < 0) sgn = -1;
4849 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
4850 {
4851 loc_col = A_ext_j[jj1];
4852 if(loc_col > -1)
4853 { /* diag */
4854 if(P_marker[loc_col] >= jj_begin_row && (sgn*A_ext_data[jj1])<0)
4855 sum += A_ext_data[jj1];
4856 }
4857 else
4858 {
4859 loc_col = -loc_col - 1;
4860 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
4861 (sgn*A_ext_data[jj1]) < 0)
4862 sum += A_ext_data[jj1];
4863 }
4864 }
4865 if(sum != 0)
4866 {
4867 distribute = A_offd_data[jj] / sum;
4868 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
4869 {
4870 loc_col = A_ext_j[jj1];
4871 if(loc_col > -1)
4872 { /* diag */
4873 if(P_marker[loc_col] >= jj_begin_row &&
4874 (sgn*A_ext_data[jj1]) < 0)
4875 P_diag_data[P_marker[loc_col]] += distribute*
4876 A_ext_data[jj1];
4877 }
4878 else
4879 {
4880 loc_col = -loc_col - 1;
4881 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
4882 (sgn*A_ext_data[jj1]) < 0)
4883 P_offd_data[P_marker_offd[loc_col]] += distribute*
4884 A_ext_data[jj1];
4885 }
4886 }
4887 }
4888 else
4889 diagonal += A_offd_data[jj];
4890 }
4891 else if (CF_marker_offd[i1] != -3)
4892 {
4893 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
4894 diagonal += A_offd_data[jj];
4895 }
4896 }
4897 }
4898 for(jj = jj_begin_row; jj < jj_end_row; jj++)
4899 P_diag_data[jj] /= -diagonal;
4900 for(jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
4901 P_offd_data[jj] /= -diagonal;
4902 }
4903 strong_f_marker--;
4904 }
4905
4906 P = hypre_ParCSRMatrixCreate(comm,
4907 hypre_ParCSRMatrixGlobalNumRows(A),
4908 total_global_cpts,
4909 hypre_ParCSRMatrixColStarts(A),
4910 num_cpts_global,
4911 0,
4912 P_diag_i[n_fine],
4913 P_offd_i[n_fine]);
4914
4915 P_diag = hypre_ParCSRMatrixDiag(P);
4916 hypre_CSRMatrixData(P_diag) = P_diag_data;
4917 hypre_CSRMatrixI(P_diag) = P_diag_i;
4918 hypre_CSRMatrixJ(P_diag) = P_diag_j;
4919 P_offd = hypre_ParCSRMatrixOffd(P);
4920 hypre_CSRMatrixData(P_offd) = P_offd_data;
4921 hypre_CSRMatrixI(P_offd) = P_offd_i;
4922 hypre_CSRMatrixJ(P_offd) = P_offd_j;
4923 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
4924
4925 /* Compress P, removing coefficients smaller than trunc_factor * Max */
4926 if (trunc_factor != 0.0 || max_elmts > 0)
4927 {
4928 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
4929 P_diag_data = hypre_CSRMatrixData(P_diag);
4930 P_diag_i = hypre_CSRMatrixI(P_diag);
4931 P_diag_j = hypre_CSRMatrixJ(P_diag);
4932 P_offd_data = hypre_CSRMatrixData(P_offd);
4933 P_offd_i = hypre_CSRMatrixI(P_offd);
4934 P_offd_j = hypre_CSRMatrixJ(P_offd);
4935 P_diag_size = P_diag_i[n_fine];
4936 P_offd_size = P_offd_i[n_fine];
4937 }
4938
4939 /* This builds col_map, col_map should be monotone increasing and contain
4940 * global numbers. */
4941 num_cols_P_offd = 0;
4942 if(P_offd_size)
4943 {
4944 num_cols_P_offd = 0;
4945 for (i=0; i < P_offd_size; i++)
4946 {
4947 index = P_offd_j[i];
4948 if(tmp_CF_marker_offd[index] == 1)
4949 {
4950 num_cols_P_offd++;
4951 tmp_CF_marker_offd[index] = 2;
4952 }
4953 }
4954 if (num_cols_P_offd)
4955 {
4956 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
4957 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
4958 }
4959
4960 index = 0;
4961 for(i = 0; i < num_cols_P_offd; i++)
4962 {
4963 while( tmp_CF_marker_offd[index] == -1) index++;
4964 col_map_offd_P[i] = fine_to_coarse_offd[index];
4965 tmp_map_offd[i] = index++;
4966 }
4967
4968 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
4969
4970 for (i = 0; i < num_cols_P_offd; i++)
4971 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
4972
4973 for(i = 0; i < P_offd_size; i++)
4974 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
4975 }
4976
4977 if (num_cols_P_offd)
4978 {
4979 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
4980 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
4981 hypre_TFree(tmp_map_offd);
4982 }
4983
4984#ifdef HYPRE_NO_GLOBAL_PARTITION
4985 hypre_NewCommPkgCreate(P);
4986#else
4987 hypre_MatvecCommPkgCreate(P);
4988#endif
4989
4990 for (i=0; i < n_fine; i++)
4991 if (CF_marker[i] == -3) CF_marker[i] = -1;
4992
4993 *P_ptr = P;
4994
4995 /* Deallocate memory */
4996 hypre_TFree(fine_to_coarse);
4997 hypre_TFree(P_marker);
4998
4999 if (num_procs > 1)
5000 {
5001
5002 hypre_BigCSRMatrixDestroy(Sop);
5003 hypre_BigCSRMatrixDestroy(A_ext);
5004 hypre_TFree(fine_to_coarse_offd);
5005 hypre_TFree(P_marker_offd);
5006 hypre_TFree(CF_marker_offd);
5007 hypre_TFree(tmp_CF_marker_offd);
5008 if(num_functions > 1)
5009 hypre_TFree(dof_func_offd);
5010 hypre_TFree(found);
5011
5012 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
5013
5014 }
5015
5016 return hypre_error_flag;
5017}
5018
5019/*---------------------------------------------------------------------------
5020 * hypre_BoomerAMGBuildExtInterp
5021 * Comment:
5022 *--------------------------------------------------------------------------*/
5023int
5024hypre_BoomerAMGBuildExtInterp(hypre_ParCSRMatrix *A, int *CF_marker,
5025 hypre_ParCSRMatrix *S, HYPRE_BigInt *num_cpts_global,
5026 int num_functions, int *dof_func, int debug_flag,
5027 double trunc_factor, int max_elmts,
5028 int *col_offd_S_to_A,
5029 hypre_ParCSRMatrix **P_ptr)
5030{
5031 /* Communication Variables */
5032 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
5033 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
5034
5035
5036 int my_id, num_procs;
5037
5038 /* Variables to store input variables */
5039 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
5040 double *A_diag_data = hypre_CSRMatrixData(A_diag);
5041 int *A_diag_i = hypre_CSRMatrixI(A_diag);
5042 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
5043
5044 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
5045 double *A_offd_data = hypre_CSRMatrixData(A_offd);
5046 int *A_offd_i = hypre_CSRMatrixI(A_offd);
5047 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
5048
5049 int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd);
5050 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
5051 int n_fine = hypre_CSRMatrixNumRows(A_diag);
5052 HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstRowIndex(A);
5053 int local_numrows = hypre_CSRMatrixNumRows(A_diag);
5054 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt) local_numrows;
5055 HYPRE_BigInt total_global_cpts, my_first_cpt;
5056
5057 /* Variables to store strong connection matrix info */
5058 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
5059 int *S_diag_i = hypre_CSRMatrixI(S_diag);
5060 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
5061
5062 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
5063 int *S_offd_i = hypre_CSRMatrixI(S_offd);
5064 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
5065
5066 /* Interpolation matrix P */
5067 hypre_ParCSRMatrix *P;
5068 hypre_CSRMatrix *P_diag;
5069 hypre_CSRMatrix *P_offd;
5070
5071 double *P_diag_data = NULL;
5072 int *P_diag_i, *P_diag_j = NULL;
5073 double *P_offd_data = NULL;
5074 int *P_offd_i, *P_offd_j = NULL;
5075
5076 HYPRE_BigInt *col_map_offd_P;
5077 int *tmp_map_offd;
5078 int P_diag_size;
5079 int P_offd_size;
5080 int *P_marker = NULL;
5081 int *P_marker_offd = NULL;
5082 int *CF_marker_offd = NULL;
5083 int *tmp_CF_marker_offd = NULL;
5084 int *dof_func_offd = NULL;
5085
5086 /* Full row information for columns of A that are off diag*/
5087 hypre_BigCSRMatrix *A_ext;
5088 double *A_ext_data;
5089 int *A_ext_i;
5090 HYPRE_BigInt *big_A_ext_j;
5091 int *A_ext_j;
5092
5093 int *fine_to_coarse = NULL;
5094 HYPRE_BigInt *fine_to_coarse_offd = NULL;
5095 HYPRE_BigInt *found;
5096
5097 int num_cols_P_offd;
5098 int newoff, loc_col;
5099 int A_ext_rows, full_off_procNodes;
5100
5101 hypre_BigCSRMatrix *Sop;
5102 int *Sop_i;
5103 HYPRE_BigInt *big_Sop_j;
5104 int *Sop_j;
5105
5106 int Soprows, sgn;
5107
5108 /* Variables to keep count of interpolatory points */
5109 int jj_counter, jj_counter_offd;
5110 int jj_begin_row, jj_end_row;
5111 int jj_begin_row_offd = 0;
5112 int jj_end_row_offd = 0;
5113 int coarse_counter, coarse_counter_offd;
5114
5115 /* Interpolation weight variables */
5116 double sum, diagonal, distribute;
5117 int strong_f_marker = -2;
5118
5119 /* Loop variables */
5120 int index;
5121 int start_indexing = 0;
5122 int i, i1, i2, jj, kk, k1, jj1;
5123
5124 /* Definitions */
5125 double zero = 0.0;
5126 double one = 1.0;
5127 double wall_time;
5128
5129
5130 hypre_ParCSRCommPkg *extend_comm_pkg = NULL;
5131
5132 if (debug_flag==4) wall_time = time_getWallclockSeconds();
5133
5134 /* BEGIN */
5135 MPI_Comm_size(comm, &num_procs);
5136 MPI_Comm_rank(comm,&my_id);
5137
5138#ifdef HYPRE_NO_GLOBAL_PARTITION
5139 my_first_cpt = num_cpts_global[0];
5140 if (my_id == (num_procs -1)) total_global_cpts = num_cpts_global[1];
5141 MPI_Bcast(&total_global_cpts, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
5142#else
5143 my_first_cpt = num_cpts_global[my_id];
5144 total_global_cpts = num_cpts_global[num_procs];
5145#endif
5146
5147 if (!comm_pkg)
5148 {
5149#ifdef HYPRE_NO_GLOBAL_PARTITION
5150 hypre_NewCommPkgCreate(A);
5151#else
5152 hypre_MatvecCommPkgCreate(A);
5153#endif
5154 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
5155 }
5156
5157 /* Set up off processor information (specifically for neighbors of
5158 * neighbors */
5159 newoff = 0;
5160 full_off_procNodes = 0;
5161 if (num_procs > 1)
5162 {
5163 /*----------------------------------------------------------------------
5164 * Get the off processors rows for A and S, associated with columns in
5165 * A_offd and S_offd.
5166 *---------------------------------------------------------------------*/
5167 A_ext = hypre_ParCSRMatrixExtractBigExt(A,A,1);
5168 A_ext_i = hypre_BigCSRMatrixI(A_ext);
5169 big_A_ext_j = hypre_BigCSRMatrixJ(A_ext);
5170 A_ext_data = hypre_BigCSRMatrixData(A_ext);
5171 A_ext_rows = hypre_BigCSRMatrixNumRows(A_ext);
5172
5173 Sop = hypre_ParCSRMatrixExtractBigExt(S,A,0);
5174 Sop_i = hypre_BigCSRMatrixI(Sop);
5175 big_Sop_j = hypre_BigCSRMatrixJ(Sop);
5176 Soprows = hypre_BigCSRMatrixNumRows(Sop);
5177
5178 /* Find nodes that are neighbors of neighbors, not found in offd */
5179 newoff = new_offd_nodes(&found, A_ext_rows, A_ext_i, big_A_ext_j,
5180 &A_ext_j,
5181 Soprows, col_map_offd, col_1, col_n,
5182 Sop_i, big_Sop_j, &Sop_j, CF_marker, comm_pkg);
5183
5184 hypre_BigCSRMatrixJ(A_ext) = NULL;
5185 hypre_BigCSRMatrixJ(Sop) = NULL;
5186
5187 if(newoff >= 0)
5188 full_off_procNodes = newoff + num_cols_A_offd;
5189 else
5190 return(1);
5191
5192 /* Possibly add new points and new processors to the comm_pkg, all
5193 * processors need new_comm_pkg */
5194
5195 /* AHB - create a new comm package just for extended info -
5196 this will work better with the assumed partition*/
5197 hypre_ParCSRFindExtendCommPkg(A, newoff, found,
5198 &extend_comm_pkg);
5199
5200 CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
5201
5202 if (num_functions > 1 && full_off_procNodes > 0)
5203 dof_func_offd = hypre_CTAlloc(int, full_off_procNodes);
5204
5205 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, CF_marker,
5206 full_off_procNodes, CF_marker_offd);
5207
5208 if(num_functions > 1)
5209 alt_insert_new_nodes(comm_pkg, extend_comm_pkg, dof_func,
5210 full_off_procNodes, dof_func_offd);
5211 }
5212
5213
5214 /*-----------------------------------------------------------------------
5215 * First Pass: Determine size of P and fill in fine_to_coarse mapping.
5216 *-----------------------------------------------------------------------*/
5217
5218 /*-----------------------------------------------------------------------
5219 * Intialize counters and allocate mapping vector.
5220 *-----------------------------------------------------------------------*/
5221
5222 P_diag_i = hypre_CTAlloc(int, n_fine+1);
5223 P_offd_i = hypre_CTAlloc(int, n_fine+1);
5224
5225 if (n_fine)
5226 {
5227 fine_to_coarse = hypre_CTAlloc(int, n_fine);
5228 P_marker = hypre_CTAlloc(int, n_fine);
5229 }
5230
5231 if (full_off_procNodes)
5232 {
5233 P_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
5234 fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, full_off_procNodes);
5235 tmp_CF_marker_offd = hypre_CTAlloc(int, full_off_procNodes);
5236 }
5237
5238 for (i=0; i < full_off_procNodes; i++)
5239 {
5240 P_marker_offd[i] = -1;
5241 tmp_CF_marker_offd[i] = -1;
5242 fine_to_coarse_offd[i] = -1;
5243 }
5244
5245 for (i=0; i < n_fine; i++)
5246 {
5247 P_marker[i] = -1;
5248 fine_to_coarse[i] = -1;
5249 }
5250
5251 jj_counter = start_indexing;
5252 jj_counter_offd = start_indexing;
5253 coarse_counter = 0;
5254 coarse_counter_offd = 0;
5255
5256 /*-----------------------------------------------------------------------
5257 * Loop over fine grid.
5258 *-----------------------------------------------------------------------*/
5259 for (i = 0; i < n_fine; i++)
5260 {
5261 P_diag_i[i] = jj_counter;
5262 if (num_procs > 1)
5263 P_offd_i[i] = jj_counter_offd;
5264
5265 if (CF_marker[i] >= 0)
5266 {
5267 jj_counter++;
5268 fine_to_coarse[i] = coarse_counter;
5269 coarse_counter++;
5270 }
5271
5272 /*--------------------------------------------------------------------
5273 * If i is an F-point, interpolation is from the C-points that
5274 * strongly influence i, or C-points that stronly influence F-points
5275 * that strongly influence i.
5276 *--------------------------------------------------------------------*/
5277 else if (CF_marker[i] != -3)
5278 {
5279 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
5280 {
5281 i1 = S_diag_j[jj];
5282 if (CF_marker[i1] > 0)
5283 { /* i1 is a C point */
5284 if (P_marker[i1] < P_diag_i[i])
5285 {
5286 P_marker[i1] = jj_counter;
5287 jj_counter++;
5288 }
5289 }
5290 else if (CF_marker[i1] != -3)
5291 { /* i1 is a F point, loop through it's strong neighbors */
5292 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
5293 {
5294 k1 = S_diag_j[kk];
5295 if (CF_marker[k1] > 0)
5296 {
5297 if(P_marker[k1] < P_diag_i[i])
5298 {
5299 P_marker[k1] = jj_counter;
5300 jj_counter++;
5301 }
5302 }
5303 }
5304 if(num_procs > 1)
5305 {
5306 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
5307 {
5308 if(col_offd_S_to_A)
5309 k1 = col_offd_S_to_A[S_offd_j[kk]];
5310 else
5311 k1 = S_offd_j[kk];
5312 if (CF_marker_offd[k1] > 0)
5313 {
5314 if(P_marker_offd[k1] < P_offd_i[i])
5315 {
5316 tmp_CF_marker_offd[k1] = 1;
5317 P_marker_offd[k1] = jj_counter_offd;
5318 jj_counter_offd++;
5319 }
5320 }
5321 }
5322 }
5323 }
5324 }
5325 /* Look at off diag strong connections of i */
5326 if (num_procs > 1)
5327 {
5328 for (jj = S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
5329 {
5330 i1 = S_offd_j[jj];
5331 if(col_offd_S_to_A)
5332 i1 = col_offd_S_to_A[i1];
5333 if (CF_marker_offd[i1] > 0)
5334 {
5335 if(P_marker_offd[i1] < P_offd_i[i])
5336 {
5337 tmp_CF_marker_offd[i1] = 1;
5338 P_marker_offd[i1] = jj_counter_offd;
5339 jj_counter_offd++;
5340 }
5341 }
5342 else if (CF_marker_offd[i1] != -3)
5343 { /* F point; look at neighbors of i1. Sop contains global col
5344 * numbers and entries that could be in S_diag or S_offd or
5345 * neither. */
5346 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
5347 {
5348 loc_col = Sop_j[kk];
5349 if(loc_col > -1)
5350 { /* In S_diag */
5351 if(CF_marker[loc_col] > 0)
5352 {
5353 if(P_marker[loc_col] < P_diag_i[i])
5354 {
5355 P_marker[loc_col] = jj_counter;
5356 jj_counter++;
5357 }
5358 }
5359 }
5360 else
5361 {
5362 loc_col = -loc_col - 1;
5363 if(CF_marker_offd[loc_col] > 0)
5364 {
5365 if(P_marker_offd[loc_col] < P_offd_i[i])
5366 {
5367 P_marker_offd[loc_col] = jj_counter_offd;
5368 tmp_CF_marker_offd[loc_col] = 1;
5369 jj_counter_offd++;
5370 }
5371 }
5372 }
5373 }
5374 }
5375 }
5376 }
5377 }
5378 }
5379
5380 if (debug_flag==4)
5381 {
5382 wall_time = time_getWallclockSeconds() - wall_time;
5383 printf("Proc = %d determine structure %f\n",
5384 my_id, wall_time);
5385 fflush(NULL);
5386 }
5387 /*-----------------------------------------------------------------------
5388 * Allocate arrays.
5389 *-----------------------------------------------------------------------*/
5390
5391 if (debug_flag== 4) wall_time = time_getWallclockSeconds();
5392
5393 P_diag_size = jj_counter;
5394 P_offd_size = jj_counter_offd;
5395 if (P_diag_size)
5396 {
5397 P_diag_j = hypre_CTAlloc(int, P_diag_size);
5398 P_diag_data = hypre_CTAlloc(double, P_diag_size);
5399 }
5400
5401 if (P_offd_size)
5402 {
5403 P_offd_j = hypre_CTAlloc(int, P_offd_size);
5404 P_offd_data = hypre_CTAlloc(double, P_offd_size);
5405 }
5406 P_diag_i[n_fine] = jj_counter;
5407 P_offd_i[n_fine] = jj_counter_offd;
5408
5409 jj_counter = start_indexing;
5410 jj_counter_offd = start_indexing;
5411
5412 /* Fine to coarse mapping */
5413 if(num_procs > 1)
5414 {
5415 big_insert_new_nodes(comm_pkg, extend_comm_pkg, fine_to_coarse,
5416 full_off_procNodes, my_first_cpt,
5417 fine_to_coarse_offd);
5418 }
5419
5420 for (i = 0; i < n_fine; i++)
5421 P_marker[i] = -1;
5422
5423 for (i = 0; i < full_off_procNodes; i++)
5424 P_marker_offd[i] = -1;
5425
5426 /*-----------------------------------------------------------------------
5427 * Loop over fine grid points.
5428 *-----------------------------------------------------------------------*/
5429 for (i = 0; i < n_fine; i++)
5430 {
5431 jj_begin_row = jj_counter;
5432 jj_begin_row_offd = jj_counter_offd;
5433
5434 /*--------------------------------------------------------------------
5435 * If i is a c-point, interpolation is the identity.
5436 *--------------------------------------------------------------------*/
5437
5438 if (CF_marker[i] >= 0)
5439 {
5440 P_diag_j[jj_counter] = fine_to_coarse[i];
5441 P_diag_data[jj_counter] = one;
5442 jj_counter++;
5443 }
5444
5445 /*--------------------------------------------------------------------
5446 * If i is an F-point, build interpolation.
5447 *--------------------------------------------------------------------*/
5448
5449 else if (CF_marker[i] != -3)
5450 {
5451 strong_f_marker--;
5452 for (jj = S_diag_i[i]; jj < S_diag_i[i+1]; jj++)
5453 {
5454 i1 = S_diag_j[jj];
5455
5456 /*--------------------------------------------------------------
5457 * If neighbor i1 is a C-point, set column number in P_diag_j
5458 * and initialize interpolation weight to zero.
5459 *--------------------------------------------------------------*/
5460
5461 if (CF_marker[i1] >= 0)
5462 {
5463 if (P_marker[i1] < jj_begin_row)
5464 {
5465 P_marker[i1] = jj_counter;
5466 P_diag_j[jj_counter] = fine_to_coarse[i1];
5467 P_diag_data[jj_counter] = zero;
5468 jj_counter++;
5469 }
5470 }
5471 else if (CF_marker[i1] != -3)
5472 {
5473 P_marker[i1] = strong_f_marker;
5474 for (kk = S_diag_i[i1]; kk < S_diag_i[i1+1]; kk++)
5475 {
5476 k1 = S_diag_j[kk];
5477 if (CF_marker[k1] > 0)
5478 {
5479 if(P_marker[k1] < jj_begin_row)
5480 {
5481 P_marker[k1] = jj_counter;
5482 P_diag_j[jj_counter] = fine_to_coarse[k1];
5483 P_diag_data[jj_counter] = zero;
5484 jj_counter++;
5485 }
5486 }
5487 }
5488 if(num_procs > 1)
5489 {
5490 for (kk = S_offd_i[i1]; kk < S_offd_i[i1+1]; kk++)
5491 {
5492 if(col_offd_S_to_A)
5493 k1 = col_offd_S_to_A[S_offd_j[kk]];
5494 else
5495 k1 = S_offd_j[kk];
5496 if(CF_marker_offd[k1] > 0)
5497 {
5498 if(P_marker_offd[k1] < jj_begin_row_offd)
5499 {
5500 P_marker_offd[k1] = jj_counter_offd;
5501 P_offd_j[jj_counter_offd] = k1;
5502 P_offd_data[jj_counter_offd] = zero;
5503 jj_counter_offd++;
5504 }
5505 }
5506 }
5507 }
5508 }
5509 }
5510
5511 if ( num_procs > 1)
5512 {
5513 for (jj=S_offd_i[i]; jj < S_offd_i[i+1]; jj++)
5514 {
5515 i1 = S_offd_j[jj];
5516 if(col_offd_S_to_A)
5517 i1 = col_offd_S_to_A[i1];
5518 if ( CF_marker_offd[i1] > 0)
5519 {
5520 if(P_marker_offd[i1] < jj_begin_row_offd)
5521 {
5522 P_marker_offd[i1] = jj_counter_offd;
5523 P_offd_j[jj_counter_offd] = i1;
5524 P_offd_data[jj_counter_offd] = zero;
5525 jj_counter_offd++;
5526 }
5527 }
5528 else if (CF_marker_offd[i1] != -3)
5529 {
5530 P_marker_offd[i1] = strong_f_marker;
5531 for(kk = Sop_i[i1]; kk < Sop_i[i1+1]; kk++)
5532 {
5533 loc_col = Sop_j[kk];
5534 if(loc_col > -1)
5535 {
5536 if(CF_marker[loc_col] > 0)
5537 {
5538 if(P_marker[loc_col] < jj_begin_row)
5539 {
5540 P_marker[loc_col] = jj_counter;
5541 P_diag_j[jj_counter] = fine_to_coarse[loc_col];
5542 P_diag_data[jj_counter] = zero;
5543 jj_counter++;
5544 }
5545 }
5546 }
5547 else
5548 {
5549 loc_col = - loc_col - 1;
5550 if(CF_marker_offd[loc_col] > 0)
5551 {
5552 if(P_marker_offd[loc_col] < jj_begin_row_offd)
5553 {
5554 P_marker_offd[loc_col] = jj_counter_offd;
5555 P_offd_j[jj_counter_offd]=loc_col;
5556 P_offd_data[jj_counter_offd] = zero;
5557 jj_counter_offd++;
5558 }
5559 }
5560 }
5561 }
5562 }
5563 }
5564 }
5565
5566 jj_end_row = jj_counter;
5567 jj_end_row_offd = jj_counter_offd;
5568
5569 diagonal = A_diag_data[A_diag_i[i]];
5570
5571 for (jj = A_diag_i[i]+1; jj < A_diag_i[i+1]; jj++)
5572 { /* i1 is a c-point and strongly influences i, accumulate
5573 * a_(i,i1) into interpolation weight */
5574 i1 = A_diag_j[jj];
5575 if (P_marker[i1] >= jj_begin_row)
5576 {
5577 P_diag_data[P_marker[i1]] += A_diag_data[jj];
5578 }
5579 else if(P_marker[i1] == strong_f_marker)
5580 {
5581 sum = zero;
5582 sgn = 1;
5583 if(A_diag_data[A_diag_i[i1]] < 0) sgn = -1;
5584 /* Loop over row of A for point i1 and calculate the sum
5585 * of the connections to c-points that strongly incluence i. */
5586 for(jj1 = A_diag_i[i1]+1; jj1 < A_diag_i[i1+1]; jj1++)
5587 {
5588 i2 = A_diag_j[jj1];
5589 if((P_marker[i2] >= jj_begin_row) && (sgn*A_diag_data[jj1]) < 0)
5590 sum += A_diag_data[jj1];
5591 }
5592 if(num_procs > 1)
5593 {
5594 for(jj1 = A_offd_i[i1]; jj1< A_offd_i[i1+1]; jj1++)
5595 {
5596 i2 = A_offd_j[jj1];
5597 if(P_marker_offd[i2] >= jj_begin_row_offd &&
5598 (sgn*A_offd_data[jj1]) < 0)
5599 sum += A_offd_data[jj1];
5600 }
5601 }
5602 if(sum != 0)
5603 {
5604 distribute = A_diag_data[jj]/sum;
5605 /* Loop over row of A for point i1 and do the distribution */
5606 for(jj1 = A_diag_i[i1]+1; jj1 < A_diag_i[i1+1]; jj1++)
5607 {
5608 i2 = A_diag_j[jj1];
5609 if(P_marker[i2] >= jj_begin_row && (sgn*A_diag_data[jj1]) < 0)
5610 P_diag_data[P_marker[i2]] +=
5611 distribute*A_diag_data[jj1];
5612 }
5613 if(num_procs > 1)
5614 {
5615 for(jj1 = A_offd_i[i1]; jj1 < A_offd_i[i1+1]; jj1++)
5616 {
5617 i2 = A_offd_j[jj1];
5618 if(P_marker_offd[i2] >= jj_begin_row_offd &&
5619 (sgn*A_offd_data[jj1]) < 0)
5620 P_offd_data[P_marker_offd[i2]] +=
5621 distribute*A_offd_data[jj1];
5622 }
5623 }
5624 }
5625 else
5626 {
5627 diagonal += A_diag_data[jj];
5628 }
5629 }
5630 /* neighbor i1 weakly influences i, accumulate a_(i,i1) into
5631 * diagonal */
5632 else if (CF_marker[i1] != -3)
5633 {
5634 if(num_functions == 1 || dof_func[i] == dof_func[i1])
5635 diagonal += A_diag_data[jj];
5636 }
5637 }
5638 if(num_procs > 1)
5639 {
5640 for(jj = A_offd_i[i]; jj < A_offd_i[i+1]; jj++)
5641 {
5642 i1 = A_offd_j[jj];
5643 if(P_marker_offd[i1] >= jj_begin_row_offd)
5644 P_offd_data[P_marker_offd[i1]] += A_offd_data[jj];
5645 else if(P_marker_offd[i1] == strong_f_marker)
5646 {
5647 sum = zero;
5648 sgn = 1;
5649 if(A_ext_data[A_ext_i[i1]] < 0) sgn = -1;
5650 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
5651 {
5652 loc_col = A_ext_j[jj1];
5653 if(loc_col > -1)
5654 { /* diag */
5655 if((P_marker[loc_col] >= jj_begin_row)
5656 && (sgn*A_ext_data[jj1]) < 0)
5657 sum += A_ext_data[jj1];
5658 }
5659 else
5660 {
5661 loc_col = - loc_col - 1;
5662 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
5663 (sgn*A_ext_data[jj1]) < 0)
5664 sum += A_ext_data[jj1];
5665 }
5666 }
5667 if(sum != 0)
5668 {
5669 distribute = A_offd_data[jj] / sum;
5670 for(jj1 = A_ext_i[i1]+1; jj1 < A_ext_i[i1+1]; jj1++)
5671 {
5672 loc_col = A_ext_j[jj1];
5673 if(loc_col > -1)
5674 { /* diag */
5675 if(P_marker[loc_col] >= jj_begin_row &&
5676 (sgn*A_ext_data[jj1]) < 0)
5677 P_diag_data[P_marker[loc_col]] += distribute*
5678 A_ext_data[jj1];
5679 }
5680 else
5681 {
5682 loc_col = - loc_col - 1;
5683 if(P_marker_offd[loc_col] >= jj_begin_row_offd &&
5684 (sgn*A_ext_data[jj1]) < 0)
5685 P_offd_data[P_marker_offd[loc_col]] += distribute*
5686 A_ext_data[jj1];
5687 }
5688 }
5689 }
5690 else
5691 {
5692 diagonal += A_offd_data[jj];
5693 }
5694 }
5695 else if (CF_marker_offd[i1] != -3)
5696 {
5697 if(num_functions == 1 || dof_func[i] == dof_func_offd[i1])
5698 diagonal += A_offd_data[jj];
5699 }
5700 }
5701 }
5702 for(jj = jj_begin_row; jj < jj_end_row; jj++)
5703 P_diag_data[jj] /= -diagonal;
5704 for(jj = jj_begin_row_offd; jj < jj_end_row_offd; jj++)
5705 P_offd_data[jj] /= -diagonal;
5706 }
5707 strong_f_marker--;
5708 }
5709
5710 if (debug_flag==4)
5711 {
5712 wall_time = time_getWallclockSeconds() - wall_time;
5713 printf("Proc = %d fill structure %f\n",
5714 my_id, wall_time);
5715 fflush(NULL);
5716 }
5717 /*-----------------------------------------------------------------------
5718 * Allocate arrays.
5719 *-----------------------------------------------------------------------*/
5720
5721 P = hypre_ParCSRMatrixCreate(comm,
5722 hypre_ParCSRMatrixGlobalNumRows(A),
5723 total_global_cpts,
5724 hypre_ParCSRMatrixColStarts(A),
5725 num_cpts_global,
5726 0,
5727 P_diag_i[n_fine],
5728 P_offd_i[n_fine]);
5729
5730 P_diag = hypre_ParCSRMatrixDiag(P);
5731 hypre_CSRMatrixData(P_diag) = P_diag_data;
5732 hypre_CSRMatrixI(P_diag) = P_diag_i;
5733 hypre_CSRMatrixJ(P_diag) = P_diag_j;
5734 P_offd = hypre_ParCSRMatrixOffd(P);
5735 hypre_CSRMatrixData(P_offd) = P_offd_data;
5736 hypre_CSRMatrixI(P_offd) = P_offd_i;
5737 hypre_CSRMatrixJ(P_offd) = P_offd_j;
5738 hypre_ParCSRMatrixOwnsRowStarts(P) = 0;
5739
5740 /* Compress P, removing coefficients smaller than trunc_factor * Max */
5741 if (trunc_factor != 0.0 || max_elmts > 0)
5742 {
5743 hypre_BoomerAMGInterpTruncation(P, trunc_factor, max_elmts);
5744 P_diag_data = hypre_CSRMatrixData(P_diag);
5745 P_diag_i = hypre_CSRMatrixI(P_diag);
5746 P_diag_j = hypre_CSRMatrixJ(P_diag);
5747 P_offd_data = hypre_CSRMatrixData(P_offd);
5748 P_offd_i = hypre_CSRMatrixI(P_offd);
5749 P_offd_j = hypre_CSRMatrixJ(P_offd);
5750 P_diag_size = P_diag_i[n_fine];
5751 P_offd_size = P_offd_i[n_fine];
5752 }
5753
5754 /* This builds col_map, col_map should be monotone increasing and contain
5755 * global numbers. */
5756 num_cols_P_offd = 0;
5757 if(P_offd_size)
5758 {
5759 num_cols_P_offd = 0;
5760 for (i=0; i < P_offd_size; i++)
5761 {
5762 index = P_offd_j[i];
5763 if (tmp_CF_marker_offd[index] == 1)
5764 {
5765 num_cols_P_offd++;
5766 tmp_CF_marker_offd[index] = 2;
5767 }
5768 }
5769
5770 if (num_cols_P_offd)
5771 {
5772 col_map_offd_P = hypre_CTAlloc(HYPRE_BigInt, num_cols_P_offd);
5773 tmp_map_offd = hypre_CTAlloc(int, num_cols_P_offd);
5774 }
5775
5776 index = 0;
5777 for(i = 0; i < num_cols_P_offd; i++)
5778 {
5779 while( tmp_CF_marker_offd[index] != 2) index++;
5780 /*printf("index %d tmp_CF %d\n", index, tmp_CF_marker_offd[index]);*/
5781 col_map_offd_P[i] = fine_to_coarse_offd[index];
5782 tmp_map_offd[i] = index++;
5783 }
5784
5785 hypre_BigQsortbi(col_map_offd_P, tmp_map_offd, 0, num_cols_P_offd-1);
5786
5787 for (i = 0; i < num_cols_P_offd; i++)
5788 tmp_CF_marker_offd[tmp_map_offd[i]] = i;
5789
5790 for(i = 0; i < P_offd_size; i++)
5791 P_offd_j[i] = tmp_CF_marker_offd[P_offd_j[i]];
5792 }
5793
5794 if (num_cols_P_offd)
5795 {
5796 hypre_ParCSRMatrixColMapOffd(P) = col_map_offd_P;
5797 hypre_CSRMatrixNumCols(P_offd) = num_cols_P_offd;
5798 hypre_TFree(tmp_map_offd);
5799 }
5800
5801#ifdef HYPRE_NO_GLOBAL_PARTITION
5802 hypre_NewCommPkgCreate(P);
5803#else
5804 hypre_MatvecCommPkgCreate(P);
5805#endif
5806
5807 for (i=0; i < n_fine; i++)
5808 if (CF_marker[i] == -3) CF_marker[i] = -1;
5809
5810 *P_ptr = P;
5811
5812 /* Deallocate memory */
5813 hypre_TFree(fine_to_coarse);
5814 hypre_TFree(P_marker);
5815
5816 if (num_procs > 1)
5817 {
5818 hypre_BigCSRMatrixDestroy(Sop);
5819 hypre_BigCSRMatrixDestroy(A_ext);
5820 hypre_TFree(fine_to_coarse_offd);
5821 hypre_TFree(P_marker_offd);
5822 hypre_TFree(CF_marker_offd);
5823 hypre_TFree(tmp_CF_marker_offd);
5824 if(num_functions > 1)
5825 hypre_TFree(dof_func_offd);
5826 hypre_TFree(found);
5827
5828
5829 hypre_MatvecCommPkgDestroy(extend_comm_pkg);
5830
5831
5832 }
5833
5834 return hypre_error_flag;
5835}
5836
Note: See TracBrowser for help on using the repository browser.