source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_coarsen.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: 78.9 KB
Line 
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
15/******************************************************************************
16 *
17 *****************************************************************************/
18
19/* following should be in a header file */
20
21
22#include "headers.h"
23
24
25
26/*==========================================================================*/
27/*==========================================================================*/
28/**
29 Selects a coarse "grid" based on the graph of a matrix.
30
31 Notes:
32 \begin{itemize}
33 \item The underlying matrix storage scheme is a hypre_ParCSR matrix.
34 \item The routine returns the following:
35 \begin{itemize}
36 \item S - a ParCSR matrix representing the "strength matrix". This is
37 used in the "build interpolation" routine.
38 \item CF\_marker - an array indicating both C-pts (value = 1) and
39 F-pts (value = -1)
40 \end{itemize}
41 \item We define the following temporary storage:
42 \begin{itemize}
43 \item measure\_array - an array containing the "measures" for each
44 of the fine-grid points
45 \item graph\_array - an array containing the list of points in the
46 "current subgraph" being considered in the coarsening process.
47 \end{itemize}
48 \item The graph of the "strength matrix" for A is a subgraph of the
49 graph of A, but requires nonsymmetric storage even if A is
50 symmetric. This is because of the directional nature of the
51 "strengh of dependence" notion (see below). Since we are using
52 nonsymmetric storage for A right now, this is not a problem. If we
53 ever add the ability to store A symmetrically, then we could store
54 the strength graph as floats instead of doubles to save space.
55 \item This routine currently "compresses" the strength matrix. We
56 should consider the possibility of defining this matrix to have the
57 same "nonzero structure" as A. To do this, we could use the same
58 A\_i and A\_j arrays, and would need only define the S\_data array.
59 There are several pros and cons to discuss.
60 \end{itemize}
61
62 Terminology:
63 \begin{itemize}
64 \item Ruge's terminology: A point is "strongly connected to" $j$, or
65 "strongly depends on" $j$, if $-a_ij >= \theta max_{l != j} \{-a_il\}$.
66 \item Here, we retain some of this terminology, but with a more
67 generalized notion of "strength". We also retain the "natural"
68 graph notation for representing the directed graph of a matrix.
69 That is, the nonzero entry $a_ij$ is represented as: i --> j. In
70 the strength matrix, S, the entry $s_ij$ is also graphically denoted
71 as above, and means both of the following:
72 \begin{itemize}
73 \item $i$ "depends on" $j$ with "strength" $s_ij$
74 \item $j$ "influences" $i$ with "strength" $s_ij$
75 \end{itemize}
76 \end{itemize}
77
78 {\bf Input files:}
79 headers.h
80
81 @return Error code.
82
83 @param A [IN]
84 coefficient matrix
85 @param strength_threshold [IN]
86 threshold parameter used to define strength
87 @param S_ptr [OUT]
88 strength matrix
89 @param CF_marker_ptr [OUT]
90 array indicating C/F points
91
92 @see */
93/*--------------------------------------------------------------------------*/
94
95#define C_PT 1
96#define F_PT -1
97#define SF_PT -3
98#define COMMON_C_PT 2
99#define Z_PT -2
100
101int
102hypre_BoomerAMGCoarsen( hypre_ParCSRMatrix *S,
103 hypre_ParCSRMatrix *A,
104 int CF_init,
105 int debug_flag,
106 int **CF_marker_ptr)
107{
108 MPI_Comm comm = hypre_ParCSRMatrixComm(S);
109 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(S);
110 hypre_ParCSRCommHandle *comm_handle;
111
112 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
113 int *S_diag_i = hypre_CSRMatrixI(S_diag);
114 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
115
116 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
117 int *S_offd_i = hypre_CSRMatrixI(S_offd);
118 int *S_offd_j;
119
120 /*HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(S);*/
121 int num_variables = hypre_CSRMatrixNumRows(S_diag);
122 /*HYPRE_BigInt col_1 = hypre_ParCSRMatrixFirstColDiag(S);
123 HYPRE_BigInt col_n = col_1 + (HYPRE_BigInt)hypre_CSRMatrixNumCols(S_diag);*/
124 int num_cols_offd = 0;
125
126 hypre_CSRMatrix *S_ext;
127 int *S_ext_i;
128 int *S_ext_j;
129
130 int num_sends = 0;
131 int *int_buf_data;
132 double *buf_data;
133
134 int *CF_marker;
135 int *CF_marker_offd;
136
137 double *measure_array;
138 int *graph_array;
139 int *graph_array_offd;
140 int graph_size;
141 HYPRE_BigInt big_graph_size;
142 int graph_offd_size;
143 HYPRE_BigInt global_graph_size;
144
145 int i, j, k, kc, jS, kS, ig, elmt;
146 int index, start, my_id, num_procs, jrow, cnt;
147
148 int ierr = 0;
149 int use_commpkg_A = 0;
150 int break_var = 1;
151
152 double wall_time;
153 int iter = 0;
154
155#if 0 /* debugging */
156 char filename[256];
157 FILE *fp;
158 int iter = 0;
159#endif
160
161 /*--------------------------------------------------------------
162 * Compute a ParCSR strength matrix, S.
163 *
164 * For now, the "strength" of dependence/influence is defined in
165 * the following way: i depends on j if
166 * aij > hypre_max (k != i) aik, aii < 0
167 * or
168 * aij < hypre_min (k != i) aik, aii >= 0
169 * Then S_ij = 1, else S_ij = 0.
170 *
171 * NOTE: the entries are negative initially, corresponding
172 * to "unaccounted-for" dependence.
173 *----------------------------------------------------------------*/
174
175 S_ext = NULL;
176 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
177 MPI_Comm_size(comm,&num_procs);
178 MPI_Comm_rank(comm,&my_id);
179
180 if (!comm_pkg)
181 {
182 use_commpkg_A = 1;
183 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
184 }
185
186 if (!comm_pkg)
187 {
188#ifdef HYPRE_NO_GLOBAL_PARTITION
189 hypre_NewCommPkgCreate(A);
190#else
191 hypre_MatvecCommPkgCreate(A);
192#endif
193 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
194 }
195
196 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
197
198 int_buf_data = hypre_CTAlloc(int, hypre_ParCSRCommPkgSendMapStart(comm_pkg,
199 num_sends));
200 buf_data = hypre_CTAlloc(double, hypre_ParCSRCommPkgSendMapStart(comm_pkg,
201 num_sends));
202
203 num_cols_offd = hypre_CSRMatrixNumCols(S_offd);
204
205 S_diag_j = hypre_CSRMatrixJ(S_diag);
206
207 if (num_cols_offd)
208 {
209 S_offd_j = hypre_CSRMatrixJ(S_offd);
210 }
211 /*----------------------------------------------------------
212 * Compute the measures
213 *
214 * The measures are currently given by the column sums of S.
215 * Hence, measure_array[i] is the number of influences
216 * of variable i.
217 *
218 * The measures are augmented by a random number
219 * between 0 and 1.
220 *----------------------------------------------------------*/
221
222 measure_array = hypre_CTAlloc(double, num_variables+num_cols_offd);
223
224 for (i=0; i < S_offd_i[num_variables]; i++)
225 {
226 measure_array[num_variables + S_offd_j[i]] += 1.0;
227 }
228 if (num_procs > 1)
229 comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg,
230 &measure_array[num_variables], buf_data);
231
232 for (i=0; i < S_diag_i[num_variables]; i++)
233 {
234 measure_array[S_diag_j[i]] += 1.0;
235 }
236
237 if (num_procs > 1)
238 hypre_ParCSRCommHandleDestroy(comm_handle);
239
240 index = 0;
241 for (i=0; i < num_sends; i++)
242 {
243 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
244 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
245 measure_array[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)]
246 += buf_data[index++];
247 }
248
249 for (i=num_variables; i < num_variables+num_cols_offd; i++)
250 {
251 measure_array[i] = 0;
252 }
253
254 /* this augments the measures */
255 if (CF_init == 2)
256 hypre_BoomerAMGIndepSetInit(S, measure_array, 1);
257 else
258 hypre_BoomerAMGIndepSetInit(S, measure_array, 0);
259
260 /*---------------------------------------------------
261 * Initialize the graph array
262 * graph_array contains interior points in elements 0 ... num_variables-1
263 * followed by boundary values
264 *---------------------------------------------------*/
265
266 graph_array = hypre_CTAlloc(int, num_variables);
267 if (num_cols_offd)
268 graph_array_offd = hypre_CTAlloc(int, num_cols_offd);
269 else
270 graph_array_offd = NULL;
271
272 /* initialize measure array and graph array */
273
274 for (ig = 0; ig < num_cols_offd; ig++)
275 graph_array_offd[ig] = ig;
276
277 /*---------------------------------------------------
278 * Initialize the C/F marker array
279 * C/F marker array contains interior points in elements 0 ...
280 * num_variables-1 followed by boundary values
281 *---------------------------------------------------*/
282
283 graph_offd_size = num_cols_offd;
284
285 if (CF_init==1)
286 {
287 CF_marker = *CF_marker_ptr;
288 cnt = 0;
289 for (i=0; i < num_variables; i++)
290 {
291 if ( (S_offd_i[i+1]-S_offd_i[i]) > 0
292 || CF_marker[i] == -1)
293 {
294 CF_marker[i] = 0;
295 }
296 if ( CF_marker[i] == Z_PT)
297 {
298 if (measure_array[i] >= 1.0 ||
299 (S_diag_i[i+1]-S_diag_i[i]) > 0)
300 {
301 CF_marker[i] = 0;
302 graph_array[cnt++] = i;
303 }
304 else
305 {
306 CF_marker[i] = F_PT;
307 }
308 }
309 else if (CF_marker[i] == SF_PT)
310 measure_array[i] = 0;
311 else
312 graph_array[cnt++] = i;
313 }
314 }
315 else
316 {
317 CF_marker = hypre_CTAlloc(int, num_variables);
318 cnt = 0;
319 for (i=0; i < num_variables; i++)
320 {
321 CF_marker[i] = 0;
322 if ( (S_diag_i[i+1]-S_diag_i[i]) == 0
323 && (S_offd_i[i+1]-S_offd_i[i]) == 0)
324 {
325 CF_marker[i] = SF_PT;
326 measure_array[i] = 0;
327 }
328 else
329 graph_array[cnt++] = i;
330 }
331 }
332 graph_size = cnt;
333 if (num_cols_offd)
334 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
335 else
336 CF_marker_offd = NULL;
337 for (i=0; i < num_cols_offd; i++)
338 CF_marker_offd[i] = 0;
339
340 /*---------------------------------------------------
341 * Loop until all points are either fine or coarse.
342 *---------------------------------------------------*/
343
344 if (num_procs > 1)
345 {
346 if (use_commpkg_A)
347 S_ext = hypre_ParCSRMatrixExtractConvBExt(S,A,0);
348 else
349 S_ext = hypre_ParCSRMatrixExtractConvBExt(S,S,0);
350 S_ext_i = hypre_CSRMatrixI(S_ext);
351 S_ext_j = hypre_CSRMatrixJ(S_ext);
352 }
353
354 /* compress S_ext and convert column numbers*/
355
356 /* index = 0;
357 for (i=0; i < num_cols_offd; i++)
358 {
359 for (j=S_ext_i[i]; j < S_ext_i[i+1]; j++)
360 {
361 big_k = S_ext_j[j];
362 if (big_k >= col_1 && big_k < col_n)
363 {
364 S_ext_j[index++] = (int)(big_k - col_1);
365 }
366 else
367 {
368 kc = hypre_BigBinarySearch(col_map_offd,big_k,num_cols_offd);
369 if (kc > -1) S_ext_j[index++] = -kc-1;
370 }
371 }
372 S_ext_i[i] = index;
373 }
374 for (i = num_cols_offd; i > 0; i--)
375 S_ext_i[i] = S_ext_i[i-1];
376 if (num_procs > 1) S_ext_i[0] = 0; */
377
378 if (debug_flag == 3)
379 {
380 wall_time = time_getWallclockSeconds() - wall_time;
381 printf("Proc = %d Initialize CLJP phase = %f\n",
382 my_id, wall_time);
383 }
384
385 while (1)
386 {
387 /*------------------------------------------------
388 * Exchange boundary data, i.i. get measures and S_ext_data
389 *------------------------------------------------*/
390
391 if (num_procs > 1)
392 comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg,
393 &measure_array[num_variables], buf_data);
394
395 if (num_procs > 1)
396 hypre_ParCSRCommHandleDestroy(comm_handle);
397
398 index = 0;
399 for (i=0; i < num_sends; i++)
400 {
401 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
402 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
403 measure_array[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)]
404 += buf_data[index++];
405 }
406
407 /*------------------------------------------------
408 * Set F-pts and update subgraph
409 *------------------------------------------------*/
410
411 if (iter || (CF_init != 1))
412 {
413 for (ig = 0; ig < graph_size; ig++)
414 {
415 i = graph_array[ig];
416
417 if ( (CF_marker[i] != C_PT) && (measure_array[i] < 1) )
418 {
419 /* set to be an F-pt */
420 CF_marker[i] = F_PT;
421
422 /* make sure all dependencies have been accounted for */
423 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
424 {
425 if (S_diag_j[jS] > -1)
426 {
427 CF_marker[i] = 0;
428 }
429 }
430 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
431 {
432 if (S_offd_j[jS] > -1)
433 {
434 CF_marker[i] = 0;
435 }
436 }
437 }
438 if (CF_marker[i])
439 {
440 measure_array[i] = 0;
441
442 /* take point out of the subgraph */
443 graph_size--;
444 graph_array[ig] = graph_array[graph_size];
445 graph_array[graph_size] = i;
446 ig--;
447 }
448 }
449 }
450
451 /*------------------------------------------------
452 * Exchange boundary data, i.i. get measures
453 *------------------------------------------------*/
454
455 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
456
457 index = 0;
458 for (i = 0; i < num_sends; i++)
459 {
460 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
461 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
462 {
463 jrow = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
464 buf_data[index++] = measure_array[jrow];
465 }
466 }
467
468 if (num_procs > 1)
469 {
470 comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg, buf_data,
471 &measure_array[num_variables]);
472
473 hypre_ParCSRCommHandleDestroy(comm_handle);
474
475 }
476 /*------------------------------------------------
477 * Debugging:
478 *
479 * Uncomment the sections of code labeled
480 * "debugging" to generate several files that
481 * can be visualized using the `coarsen.m'
482 * matlab routine.
483 *------------------------------------------------*/
484
485#if 0 /* debugging */
486 /* print out measures */
487 sprintf(filename, "coarsen.out.measures.%04d", iter);
488 fp = fopen(filename, "w");
489 for (i = 0; i < num_variables; i++)
490 {
491 fprintf(fp, "%f\n", measure_array[i]);
492 }
493 fclose(fp);
494
495 /* print out strength matrix */
496 sprintf(filename, "coarsen.out.strength.%04d", iter);
497 hypre_CSRMatrixPrint(S, filename);
498
499 /* print out C/F marker */
500 sprintf(filename, "coarsen.out.CF.%04d", iter);
501 fp = fopen(filename, "w");
502 for (i = 0; i < num_variables; i++)
503 {
504 fprintf(fp, "%d\n", CF_marker[i]);
505 }
506 fclose(fp);
507
508 iter++;
509#endif
510
511 /*------------------------------------------------
512 * Test for convergence
513 *------------------------------------------------*/
514 big_graph_size = (HYPRE_BigInt) graph_size;
515
516 MPI_Allreduce(&big_graph_size,&global_graph_size,1,MPI_HYPRE_BIG_INT,MPI_SUM,comm);
517
518 if (global_graph_size == 0)
519 break;
520
521 /*------------------------------------------------
522 * Pick an independent set of points with
523 * maximal measure.
524 *------------------------------------------------*/
525 if (iter || (CF_init != 1))
526 {
527 hypre_BoomerAMGIndepSet(S, measure_array, graph_array,
528 graph_size,
529 graph_array_offd, graph_offd_size,
530 CF_marker, CF_marker_offd);
531 if (num_procs > 1)
532 {
533 comm_handle = hypre_ParCSRCommHandleCreate(12, comm_pkg,
534 CF_marker_offd, int_buf_data);
535
536 hypre_ParCSRCommHandleDestroy(comm_handle);
537 }
538
539 index = 0;
540 for (i = 0; i < num_sends; i++)
541 {
542 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
543 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg,i+1);j++) {
544 elmt = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
545 if (!int_buf_data[index++] && CF_marker[elmt] > 0)
546 {
547 CF_marker[elmt] = 0;
548 }
549 }
550 }
551 }
552
553 iter++;
554 /*------------------------------------------------
555 * Exchange boundary data for CF_marker
556 *------------------------------------------------*/
557
558
559 index = 0;
560 for (i = 0; i < num_sends; i++)
561 {
562 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
563 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
564 {
565 elmt = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
566 int_buf_data[index++] = CF_marker[elmt];
567 }
568 }
569
570 if (num_procs > 1)
571 {
572 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
573 CF_marker_offd);
574
575 hypre_ParCSRCommHandleDestroy(comm_handle);
576 }
577
578 for (ig = 0; ig < graph_offd_size; ig++)
579 {
580 i = graph_array_offd[ig];
581
582 if (CF_marker_offd[i] < 0)
583 {
584 /* take point out of the subgraph */
585 graph_offd_size--;
586 graph_array_offd[ig] = graph_array_offd[graph_offd_size];
587 graph_array_offd[graph_offd_size] = i;
588 ig--;
589 }
590 }
591 if (debug_flag == 3)
592 {
593 wall_time = time_getWallclockSeconds() - wall_time;
594 printf("Proc = %d iter %d comm. and subgraph update = %f\n",
595 my_id, iter, wall_time);
596 }
597 /*------------------------------------------------
598 * Set C_pts and apply heuristics.
599 *------------------------------------------------*/
600
601 for (i=num_variables; i < num_variables+num_cols_offd; i++)
602 {
603 measure_array[i] = 0;
604 }
605
606 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
607 for (ig = 0; ig < graph_size; ig++)
608 {
609 i = graph_array[ig];
610
611 /*---------------------------------------------
612 * Heuristic: C-pts don't interpolate from
613 * neighbors that influence them.
614 *---------------------------------------------*/
615
616 if (CF_marker[i] > 0)
617 {
618 /* set to be a C-pt */
619 CF_marker[i] = C_PT;
620
621 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
622 {
623 j = S_diag_j[jS];
624 if (j > -1)
625 {
626
627 /* "remove" edge from S */
628 S_diag_j[jS] = -S_diag_j[jS]-1;
629
630 /* decrement measures of unmarked neighbors */
631 if (!CF_marker[j])
632 {
633 measure_array[j]--;
634 }
635 }
636 }
637 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
638 {
639 j = S_offd_j[jS];
640 if (j > -1)
641 {
642
643 /* "remove" edge from S */
644 S_offd_j[jS] = -S_offd_j[jS]-1;
645
646 /* decrement measures of unmarked neighbors */
647 if (!CF_marker_offd[j])
648 {
649 measure_array[j+num_variables]--;
650 }
651 }
652 }
653 }
654 else
655 {
656 /* marked dependencies */
657 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
658 {
659 j = S_diag_j[jS];
660 if (j < 0) j = -j-1;
661
662 if (CF_marker[j] > 0)
663 {
664 if (S_diag_j[jS] > -1)
665 {
666 /* "remove" edge from S */
667 S_diag_j[jS] = -S_diag_j[jS]-1;
668 }
669
670 /* IMPORTANT: consider all dependencies */
671 /* temporarily modify CF_marker */
672 CF_marker[j] = COMMON_C_PT;
673 }
674 else if (CF_marker[j] == SF_PT)
675 {
676 if (S_diag_j[jS] > -1)
677 {
678 /* "remove" edge from S */
679 S_diag_j[jS] = -S_diag_j[jS]-1;
680 }
681 }
682 }
683 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
684 {
685 j = S_offd_j[jS];
686 if (j < 0) j = -j-1;
687
688 if (CF_marker_offd[j] > 0)
689 {
690 if (S_offd_j[jS] > -1)
691 {
692 /* "remove" edge from S */
693 S_offd_j[jS] = -S_offd_j[jS]-1;
694 }
695
696 /* IMPORTANT: consider all dependencies */
697 /* temporarily modify CF_marker */
698 CF_marker_offd[j] = COMMON_C_PT;
699 }
700 else if (CF_marker_offd[j] == SF_PT)
701 {
702 if (S_offd_j[jS] > -1)
703 {
704 /* "remove" edge from S */
705 S_offd_j[jS] = -S_offd_j[jS]-1;
706 }
707 }
708 }
709
710 /* unmarked dependencies */
711 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
712 {
713 if (S_diag_j[jS] > -1)
714 {
715 j = S_diag_j[jS];
716 break_var = 1;
717 /* check for common C-pt */
718 for (kS = S_diag_i[j]; kS < S_diag_i[j+1]; kS++)
719 {
720 k = S_diag_j[kS];
721 if (k < 0) k = -k-1;
722
723 /* IMPORTANT: consider all dependencies */
724 if (CF_marker[k] == COMMON_C_PT)
725 {
726 /* "remove" edge from S and update measure*/
727 S_diag_j[jS] = -S_diag_j[jS]-1;
728 measure_array[j]--;
729 break_var = 0;
730 break;
731 }
732 }
733 if (break_var)
734 {
735 for (kS = S_offd_i[j]; kS < S_offd_i[j+1]; kS++)
736 {
737 k = S_offd_j[kS];
738 if (k < 0) k = -k-1;
739
740 /* IMPORTANT: consider all dependencies */
741 if ( CF_marker_offd[k] == COMMON_C_PT)
742 {
743 /* "remove" edge from S and update measure*/
744 S_diag_j[jS] = -S_diag_j[jS]-1;
745 measure_array[j]--;
746 break;
747 }
748 }
749 }
750 }
751 }
752 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
753 {
754 if (S_offd_j[jS] > -1)
755 {
756 j = S_offd_j[jS];
757
758 /* check for common C-pt */
759 for (kS = S_ext_i[j]; kS < S_ext_i[j+1]; kS++)
760 {
761 k = S_ext_j[kS];
762 if (k >= 0)
763 {
764 /* IMPORTANT: consider all dependencies */
765 if (CF_marker[k] == COMMON_C_PT)
766 {
767 /* "remove" edge from S and update measure*/
768 S_offd_j[jS] = -S_offd_j[jS]-1;
769 measure_array[j+num_variables]--;
770 break;
771 }
772 }
773 else
774 {
775 kc = -k-1;
776 if (kc > -1 && CF_marker_offd[kc] == COMMON_C_PT)
777 {
778 /* "remove" edge from S and update measure*/
779 S_offd_j[jS] = -S_offd_j[jS]-1;
780 measure_array[j+num_variables]--;
781 break;
782 }
783 }
784 }
785 }
786 }
787 }
788
789 /* reset CF_marker */
790 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
791 {
792 j = S_diag_j[jS];
793 if (j < 0) j = -j-1;
794
795 if (CF_marker[j] == COMMON_C_PT)
796 {
797 CF_marker[j] = C_PT;
798 }
799 }
800 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
801 {
802 j = S_offd_j[jS];
803 if (j < 0) j = -j-1;
804
805 if (CF_marker_offd[j] == COMMON_C_PT)
806 {
807 CF_marker_offd[j] = C_PT;
808 }
809 }
810 }
811 if (debug_flag == 3)
812 {
813 wall_time = time_getWallclockSeconds() - wall_time;
814 printf("Proc = %d CLJP phase = %f graph_size = %d nc_offd = %d\n",
815 my_id, wall_time, graph_size, num_cols_offd);
816 }
817 }
818
819 /*---------------------------------------------------
820 * Clean up and return
821 *---------------------------------------------------*/
822
823 /* Reset S_matrix */
824 for (i=0; i < S_diag_i[num_variables]; i++)
825 {
826 if (S_diag_j[i] < 0)
827 S_diag_j[i] = -S_diag_j[i]-1;
828 }
829 for (i=0; i < S_offd_i[num_variables]; i++)
830 {
831 if (S_offd_j[i] < 0)
832 S_offd_j[i] = -S_offd_j[i]-1;
833 }
834 /*for (i=0; i < num_variables; i++)
835 if (CF_marker[i] == SF_PT) CF_marker[i] = F_PT;*/
836
837 hypre_TFree(measure_array);
838 hypre_TFree(graph_array);
839 if (num_cols_offd) hypre_TFree(graph_array_offd);
840 hypre_TFree(buf_data);
841 hypre_TFree(int_buf_data);
842 hypre_TFree(CF_marker_offd);
843 if (num_procs > 1) hypre_CSRMatrixDestroy(S_ext);
844
845 *CF_marker_ptr = CF_marker;
846
847 return (ierr);
848}
849
850/*==========================================================================
851 * Ruge's coarsening algorithm
852 *==========================================================================*/
853
854#define C_PT 1
855#define F_PT -1
856#define Z_PT -2
857#define SF_PT -3 /* special fine points */
858#define SC_PT 3 /* special coarse points */
859#define UNDECIDED 0
860
861
862/**************************************************************
863 *
864 * Ruge Coarsening routine
865 *
866 **************************************************************/
867int
868hypre_BoomerAMGCoarsenRuge( hypre_ParCSRMatrix *S,
869 hypre_ParCSRMatrix *A,
870 int measure_type,
871 int coarsen_type,
872 int debug_flag,
873 int **CF_marker_ptr)
874{
875 MPI_Comm comm = hypre_ParCSRMatrixComm(S);
876 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(S);
877 hypre_ParCSRCommHandle *comm_handle;
878 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
879 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
880 int *S_i = hypre_CSRMatrixI(S_diag);
881 int *S_j = hypre_CSRMatrixJ(S_diag);
882 int *S_offd_i = hypre_CSRMatrixI(S_offd);
883 int *S_offd_j;
884 int num_variables = hypre_CSRMatrixNumRows(S_diag);
885 int num_cols_offd = hypre_CSRMatrixNumCols(S_offd);
886 /*HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(S);*/
887
888 hypre_CSRMatrix *S_ext;
889 int *S_ext_i;
890 int *S_ext_j;
891
892 hypre_CSRMatrix *ST;
893 int *ST_i;
894 int *ST_j;
895
896 int *CF_marker;
897 int *CF_marker_offd;
898 int ci_tilde = -1;
899 int ci_tilde_mark = -1;
900 int ci_tilde_offd = -1;
901 int ci_tilde_offd_mark = -1;
902
903 int *measure_array;
904 int *graph_array;
905 int *int_buf_data = NULL;
906 int *ci_array;
907
908 int i, j, k, jS;
909 int ji, jj, jk, jm, index;
910 int set_empty = 1;
911 int C_i_nonempty = 0;
912 /*int num_nonzeros;*/
913 int num_procs, my_id;
914 int num_sends = 0;
915 int start;
916 /*HYPRE_BigInt first_col;
917 HYPRE_BigInt col_0, col_n;*/
918
919 hypre_LinkList LoL_head;
920 hypre_LinkList LoL_tail;
921
922 int *lists, *where;
923 int measure, new_meas;
924 int meas_type = 0;
925 int agg_2 = 0;
926 int num_left, elmt;
927 int nabor, nabor_two;
928
929 int ierr = 0;
930 int use_commpkg_A = 0;
931 int break_var = 0;
932 int f_pnt = F_PT;
933 double wall_time;
934
935 if (coarsen_type < 0) coarsen_type = -coarsen_type;
936 if (measure_type == 1 || measure_type == 4) meas_type = 1;
937 if (measure_type == 4 || measure_type == 3) agg_2 = 1;
938
939 /*-------------------------------------------------------
940 * Initialize the C/F marker, LoL_head, LoL_tail arrays
941 *-------------------------------------------------------*/
942
943 LoL_head = NULL;
944 LoL_tail = NULL;
945 lists = hypre_CTAlloc(int, num_variables);
946 where = hypre_CTAlloc(int, num_variables);
947
948#if 0 /* debugging */
949 char filename[256];
950 FILE *fp;
951 int iter = 0;
952#endif
953
954 /*--------------------------------------------------------------
955 * Compute a CSR strength matrix, S.
956 *
957 * For now, the "strength" of dependence/influence is defined in
958 * the following way: i depends on j if
959 * aij > hypre_max (k != i) aik, aii < 0
960 * or
961 * aij < hypre_min (k != i) aik, aii >= 0
962 * Then S_ij = 1, else S_ij = 0.
963 *
964 * NOTE: the entries are negative initially, corresponding
965 * to "unaccounted-for" dependence.
966 *----------------------------------------------------------------*/
967
968 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
969
970 MPI_Comm_size(comm,&num_procs);
971 MPI_Comm_rank(comm,&my_id);
972
973 if (!comm_pkg)
974 {
975 use_commpkg_A = 1;
976 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
977 }
978
979 if (!comm_pkg)
980 {
981#ifdef HYPRE_NO_GLOBAL_PARTITION
982 hypre_NewCommPkgCreate(A);
983#else
984 hypre_MatvecCommPkgCreate(A);
985#endif
986 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
987 }
988
989 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
990
991 if (num_cols_offd) S_offd_j = hypre_CSRMatrixJ(S_offd);
992
993 jS = S_i[num_variables];
994
995 ST = hypre_CSRMatrixCreate(num_variables, num_variables, jS);
996 ST_i = hypre_CTAlloc(int,num_variables+1);
997 ST_j = hypre_CTAlloc(int,jS);
998 hypre_CSRMatrixI(ST) = ST_i;
999 hypre_CSRMatrixJ(ST) = ST_j;
1000
1001 /*----------------------------------------------------------
1002 * generate transpose of S, ST
1003 *----------------------------------------------------------*/
1004
1005 for (i=0; i <= num_variables; i++)
1006 ST_i[i] = 0;
1007
1008 for (i=0; i < jS; i++)
1009 {
1010 ST_i[S_j[i]+1]++;
1011 }
1012 for (i=0; i < num_variables; i++)
1013 {
1014 ST_i[i+1] += ST_i[i];
1015 }
1016 for (i=0; i < num_variables; i++)
1017 {
1018 for (j=S_i[i]; j < S_i[i+1]; j++)
1019 {
1020 index = S_j[j];
1021 ST_j[ST_i[index]] = i;
1022 ST_i[index]++;
1023 }
1024 }
1025 for (i = num_variables; i > 0; i--)
1026 {
1027 ST_i[i] = ST_i[i-1];
1028 }
1029 ST_i[0] = 0;
1030
1031 /*----------------------------------------------------------
1032 * Compute the measures
1033 *
1034 * The measures are given by the row sums of ST.
1035 * Hence, measure_array[i] is the number of influences
1036 * of variable i.
1037 * correct actual measures through adding influences from
1038 * neighbor processors
1039 *----------------------------------------------------------*/
1040
1041 measure_array = hypre_CTAlloc(int, num_variables);
1042
1043 for (i = 0; i < num_variables; i++)
1044 {
1045 measure_array[i] = ST_i[i+1]-ST_i[i];
1046 }
1047
1048 /* special case for Falgout coarsening */
1049 if (coarsen_type == 6)
1050 {
1051 f_pnt = Z_PT;
1052 coarsen_type = 1;
1053 }
1054 if (coarsen_type == 10)
1055 {
1056 f_pnt = Z_PT;
1057 coarsen_type = 11;
1058 }
1059
1060 if (meas_type && num_procs > 1)
1061 {
1062 int *measure_offd;
1063 measure_offd = hypre_CTAlloc(int, num_cols_offd);
1064 int_buf_data = hypre_CTAlloc(int,
1065 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1066
1067 for (i=0; i < S_offd_i[num_variables]; i++)
1068 measure_offd[S_offd_j[i]]++;
1069
1070 comm_handle = hypre_ParCSRCommHandleCreate(12, comm_pkg,
1071 measure_offd, int_buf_data);
1072
1073 hypre_ParCSRCommHandleDestroy(comm_handle);
1074
1075 index = 0;
1076 for (i=0; i < num_sends; i++)
1077 {
1078 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1079 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1080 measure_array[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)]
1081 += int_buf_data[index++];
1082 }
1083 hypre_TFree(measure_offd);
1084 }
1085
1086 if ((coarsen_type != 1 && coarsen_type != 11)
1087 && num_procs > 1)
1088 {
1089 if (use_commpkg_A)
1090 S_ext = hypre_ParCSRMatrixExtractConvBExt(S,A,0);
1091 else
1092 S_ext = hypre_ParCSRMatrixExtractConvBExt(S,S,0);
1093 S_ext_i = hypre_CSRMatrixI(S_ext);
1094 S_ext_j = hypre_CSRMatrixJ(S_ext);
1095 /*num_nonzeros = S_ext_i[num_cols_offd];
1096 col_0 = hypre_ParCSRMatrixFirstColDiag(S)-1;
1097 col_n = col_0 + (HYPRE_BigInt) num_variables;
1098 if (measure_type)
1099 {
1100 for (i=0; i < num_nonzeros; i++)
1101 {
1102 index = S_ext_j[i] - first_col;
1103 if (index > -1 && index < num_variables)
1104 measure_array[index]++;
1105 }
1106 } */
1107 }
1108
1109 /*---------------------------------------------------
1110 * Loop until all points are either fine or coarse.
1111 *---------------------------------------------------*/
1112
1113 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
1114
1115 /* first coarsening phase */
1116
1117 /*************************************************************
1118 *
1119 * Initialize the lists
1120 *
1121 *************************************************************/
1122
1123 CF_marker = hypre_CTAlloc(int, num_variables);
1124
1125 num_left = 0;
1126 for (j = 0; j < num_variables; j++)
1127 {
1128 if ((S_i[j+1]-S_i[j])== 0 &&
1129 (S_offd_i[j+1]-S_offd_i[j]) == 0)
1130 {
1131 CF_marker[j] = SF_PT;
1132 if (agg_2) CF_marker[j] = SC_PT;
1133 measure_array[j] = 0;
1134 }
1135 else
1136 {
1137 CF_marker[j] = UNDECIDED;
1138 num_left++;
1139 }
1140 }
1141
1142 for (j = 0; j < num_variables; j++)
1143 {
1144 measure = measure_array[j];
1145 if (CF_marker[j] != SF_PT && CF_marker[j] != SC_PT)
1146 {
1147 if (measure > 0)
1148 {
1149 enter_on_lists(&LoL_head, &LoL_tail, measure, j, lists, where);
1150 }
1151 else
1152 {
1153 if (measure < 0) printf("negative measure!\n");
1154 CF_marker[j] = f_pnt;
1155 for (k = S_i[j]; k < S_i[j+1]; k++)
1156 {
1157 nabor = S_j[k];
1158 if (CF_marker[nabor] != SF_PT && CF_marker[nabor] != SC_PT)
1159 {
1160 if (nabor < j)
1161 {
1162 new_meas = measure_array[nabor];
1163 if (new_meas > 0)
1164 remove_point(&LoL_head, &LoL_tail, new_meas,
1165 nabor, lists, where);
1166
1167 new_meas = ++(measure_array[nabor]);
1168 enter_on_lists(&LoL_head, &LoL_tail, new_meas,
1169 nabor, lists, where);
1170 }
1171 else
1172 {
1173 new_meas = ++(measure_array[nabor]);
1174 }
1175 }
1176 }
1177 --num_left;
1178 }
1179 }
1180 }
1181
1182 /****************************************************************
1183 *
1184 * Main loop of Ruge-Stueben first coloring pass.
1185 *
1186 * WHILE there are still points to classify DO:
1187 * 1) find first point, i, on list with max_measure
1188 * make i a C-point, remove it from the lists
1189 * 2) For each point, j, in S_i^T,
1190 * a) Set j to be an F-point
1191 * b) For each point, k, in S_j
1192 * move k to the list in LoL with measure one
1193 * greater than it occupies (creating new LoL
1194 * entry if necessary)
1195 * 3) For each point, j, in S_i,
1196 * move j to the list in LoL with measure one
1197 * smaller than it occupies (creating new LoL
1198 * entry if necessary)
1199 *
1200 ****************************************************************/
1201
1202 while (num_left > 0)
1203 {
1204 index = LoL_head -> head;
1205
1206 CF_marker[index] = C_PT;
1207 measure = measure_array[index];
1208 measure_array[index] = 0;
1209 --num_left;
1210
1211 remove_point(&LoL_head, &LoL_tail, measure, index, lists, where);
1212
1213 for (j = ST_i[index]; j < ST_i[index+1]; j++)
1214 {
1215 nabor = ST_j[j];
1216 if (CF_marker[nabor] == UNDECIDED)
1217 {
1218 CF_marker[nabor] = F_PT;
1219 measure = measure_array[nabor];
1220
1221 remove_point(&LoL_head, &LoL_tail, measure, nabor, lists, where);
1222 --num_left;
1223
1224 for (k = S_i[nabor]; k < S_i[nabor+1]; k++)
1225 {
1226 nabor_two = S_j[k];
1227 if (CF_marker[nabor_two] == UNDECIDED)
1228 {
1229 measure = measure_array[nabor_two];
1230 remove_point(&LoL_head, &LoL_tail, measure,
1231 nabor_two, lists, where);
1232
1233 new_meas = ++(measure_array[nabor_two]);
1234
1235 enter_on_lists(&LoL_head, &LoL_tail, new_meas,
1236 nabor_two, lists, where);
1237 }
1238 }
1239 }
1240 }
1241 for (j = S_i[index]; j < S_i[index+1]; j++)
1242 {
1243 nabor = S_j[j];
1244 if (CF_marker[nabor] == UNDECIDED)
1245 {
1246 measure = measure_array[nabor];
1247
1248 remove_point(&LoL_head, &LoL_tail, measure, nabor, lists, where);
1249
1250 measure_array[nabor] = --measure;
1251
1252 if (measure > 0)
1253 enter_on_lists(&LoL_head, &LoL_tail, measure, nabor,
1254 lists, where);
1255 else
1256 {
1257 CF_marker[nabor] = F_PT;
1258 --num_left;
1259
1260 for (k = S_i[nabor]; k < S_i[nabor+1]; k++)
1261 {
1262 nabor_two = S_j[k];
1263 if (CF_marker[nabor_two] == UNDECIDED)
1264 {
1265 new_meas = measure_array[nabor_two];
1266 remove_point(&LoL_head, &LoL_tail, new_meas,
1267 nabor_two, lists, where);
1268
1269 new_meas = ++(measure_array[nabor_two]);
1270
1271 enter_on_lists(&LoL_head, &LoL_tail, new_meas,
1272 nabor_two, lists, where);
1273 }
1274 }
1275 }
1276 }
1277 }
1278 }
1279
1280 hypre_TFree(measure_array);
1281 hypre_CSRMatrixDestroy(ST);
1282
1283 if (debug_flag == 3)
1284 {
1285 wall_time = time_getWallclockSeconds() - wall_time;
1286 printf("Proc = %d Coarsen 1st pass = %f\n",
1287 my_id, wall_time);
1288 }
1289
1290 hypre_TFree(lists);
1291 hypre_TFree(where);
1292 hypre_TFree(LoL_head);
1293 hypre_TFree(LoL_tail);
1294
1295 for (i=0; i < num_variables; i++)
1296 if (CF_marker[i] == SC_PT) CF_marker[i] = C_PT;
1297
1298 if (coarsen_type == 11)
1299 {
1300 hypre_TFree(int_buf_data);
1301 int_buf_data = NULL;
1302 *CF_marker_ptr = CF_marker;
1303 return 0;
1304 }
1305
1306 /* second pass, check fine points for coarse neighbors
1307 for coarsen_type = 2, the second pass includes
1308 off-processore boundary points */
1309
1310 /*---------------------------------------------------
1311 * Initialize the graph array
1312 *---------------------------------------------------*/
1313
1314 graph_array = hypre_CTAlloc(int, num_variables);
1315
1316 for (i = 0; i < num_variables; i++)
1317 {
1318 graph_array[i] = -1;
1319 }
1320
1321 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
1322
1323 if (coarsen_type == 2)
1324 {
1325 /*------------------------------------------------
1326 * Exchange boundary data for CF_marker
1327 *------------------------------------------------*/
1328
1329 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
1330 if (int_buf_data == NULL)
1331 int_buf_data = hypre_CTAlloc(int,
1332 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1333
1334 index = 0;
1335 for (i = 0; i < num_sends; i++)
1336 {
1337 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1338 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1339 int_buf_data[index++]
1340 = CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
1341 }
1342
1343 if (num_procs > 1)
1344 {
1345 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
1346 CF_marker_offd);
1347
1348 hypre_ParCSRCommHandleDestroy(comm_handle);
1349 }
1350
1351 ci_array = hypre_CTAlloc(int,num_cols_offd);
1352 for (i=0; i < num_cols_offd; i++)
1353 ci_array[i] = -1;
1354
1355 for (i=0; i < num_variables; i++)
1356 {
1357 if (ci_tilde_mark |= i) ci_tilde = -1;
1358 if (ci_tilde_offd_mark |= i) ci_tilde_offd = -1;
1359 if (CF_marker[i] == -1)
1360 {
1361 break_var = 1;
1362 for (ji = S_i[i]; ji < S_i[i+1]; ji++)
1363 {
1364 j = S_j[ji];
1365 if (CF_marker[j] > 0)
1366 graph_array[j] = i;
1367 }
1368 for (ji = S_offd_i[i]; ji < S_offd_i[i+1]; ji++)
1369 {
1370 j = S_offd_j[ji];
1371 if (CF_marker_offd[j] > 0)
1372 ci_array[j] = i;
1373 }
1374 for (ji = S_i[i]; ji < S_i[i+1]; ji++)
1375 {
1376 j = S_j[ji];
1377 if (CF_marker[j] == -1)
1378 {
1379 set_empty = 1;
1380 for (jj = S_i[j]; jj < S_i[j+1]; jj++)
1381 {
1382 index = S_j[jj];
1383 if (graph_array[index] == i)
1384 {
1385 set_empty = 0;
1386 break;
1387 }
1388 }
1389 if (set_empty)
1390 {
1391 for (jj = S_offd_i[j]; jj < S_offd_i[j+1]; jj++)
1392 {
1393 index = S_offd_j[jj];
1394 if (ci_array[index] == i)
1395 {
1396 set_empty = 0;
1397 break;
1398 }
1399 }
1400 }
1401 if (set_empty)
1402 {
1403 if (C_i_nonempty)
1404 {
1405 CF_marker[i] = 1;
1406 if (ci_tilde > -1)
1407 {
1408 CF_marker[ci_tilde] = -1;
1409 ci_tilde = -1;
1410 }
1411 if (ci_tilde_offd > -1)
1412 {
1413 CF_marker_offd[ci_tilde_offd] = -1;
1414 ci_tilde_offd = -1;
1415 }
1416 C_i_nonempty = 0;
1417 break_var = 0;
1418 break;
1419 }
1420 else
1421 {
1422 ci_tilde = j;
1423 ci_tilde_mark = i;
1424 CF_marker[j] = 1;
1425 C_i_nonempty = 1;
1426 i--;
1427 break_var = 0;
1428 break;
1429 }
1430 }
1431 }
1432 }
1433 if (break_var)
1434 {
1435 for (ji = S_offd_i[i]; ji < S_offd_i[i+1]; ji++)
1436 {
1437 j = S_offd_j[ji];
1438 if (CF_marker_offd[j] == -1)
1439 {
1440 set_empty = 1;
1441 for (jj = S_ext_i[j]; jj < S_ext_i[j+1]; jj++)
1442 {
1443 index = S_ext_j[jj];
1444 if (index >= 0) /* index interior */
1445 {
1446 if (graph_array[index] == i)
1447 {
1448 set_empty = 0;
1449 break;
1450 }
1451 }
1452 else
1453 {
1454 jk = -index-1;
1455 if (ci_array[jk] == i)
1456 {
1457 set_empty = 0;
1458 break;
1459 }
1460 }
1461 }
1462 if (set_empty)
1463 {
1464 if (C_i_nonempty)
1465 {
1466 CF_marker[i] = 1;
1467 if (ci_tilde > -1)
1468 {
1469 CF_marker[ci_tilde] = -1;
1470 ci_tilde = -1;
1471 }
1472 if (ci_tilde_offd > -1)
1473 {
1474 CF_marker_offd[ci_tilde_offd] = -1;
1475 ci_tilde_offd = -1;
1476 }
1477 C_i_nonempty = 0;
1478 break;
1479 }
1480 else
1481 {
1482 ci_tilde_offd = j;
1483 ci_tilde_offd_mark = i;
1484 CF_marker_offd[j] = 1;
1485 C_i_nonempty = 1;
1486 i--;
1487 break;
1488 }
1489 }
1490 }
1491 }
1492 }
1493 }
1494 }
1495 }
1496 else
1497 {
1498 for (i=0; i < num_variables; i++)
1499 {
1500 if (ci_tilde_mark |= i) ci_tilde = -1;
1501 if (CF_marker[i] == -1)
1502 {
1503 for (ji = S_i[i]; ji < S_i[i+1]; ji++)
1504 {
1505 j = S_j[ji];
1506 if (CF_marker[j] > 0)
1507 graph_array[j] = i;
1508 }
1509 for (ji = S_i[i]; ji < S_i[i+1]; ji++)
1510 {
1511 j = S_j[ji];
1512 if (CF_marker[j] == -1)
1513 {
1514 set_empty = 1;
1515 for (jj = S_i[j]; jj < S_i[j+1]; jj++)
1516 {
1517 index = S_j[jj];
1518 if (graph_array[index] == i)
1519 {
1520 set_empty = 0;
1521 break;
1522 }
1523 }
1524 if (set_empty)
1525 {
1526 if (C_i_nonempty)
1527 {
1528 CF_marker[i] = 1;
1529 if (ci_tilde > -1)
1530 {
1531 CF_marker[ci_tilde] = -1;
1532 ci_tilde = -1;
1533 }
1534 C_i_nonempty = 0;
1535 break;
1536 }
1537 else
1538 {
1539 ci_tilde = j;
1540 ci_tilde_mark = i;
1541 CF_marker[j] = 1;
1542 C_i_nonempty = 1;
1543 i--;
1544 break;
1545 }
1546 }
1547 }
1548 }
1549 }
1550 }
1551 }
1552
1553 if (debug_flag == 3 && coarsen_type != 2)
1554 {
1555 wall_time = time_getWallclockSeconds() - wall_time;
1556 printf("Proc = %d Coarsen 2nd pass = %f\n",
1557 my_id, wall_time);
1558 }
1559
1560 /* third pass, check boundary fine points for coarse neighbors */
1561
1562 if (coarsen_type == 3 || coarsen_type == 4)
1563 {
1564 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
1565
1566 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
1567 if (int_buf_data == NULL)
1568 int_buf_data = hypre_CTAlloc(int,
1569 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1570
1571 /*------------------------------------------------
1572 * Exchange boundary data for CF_marker
1573 *------------------------------------------------*/
1574
1575 index = 0;
1576 for (i = 0; i < num_sends; i++)
1577 {
1578 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1579 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1580 int_buf_data[index++]
1581 = CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
1582 }
1583
1584 if (num_procs > 1)
1585 {
1586 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
1587 int_buf_data, CF_marker_offd);
1588
1589 hypre_ParCSRCommHandleDestroy(comm_handle);
1590 }
1591
1592 ci_array = hypre_CTAlloc(int,num_cols_offd);
1593 for (i=0; i < num_cols_offd; i++)
1594 ci_array[i] = -1;
1595 }
1596
1597 if (coarsen_type > 1 && coarsen_type < 5)
1598 {
1599 for (i=0; i < num_variables; i++)
1600 graph_array[i] = -1;
1601 for (i=0; i < num_cols_offd; i++)
1602 {
1603 if (ci_tilde_mark |= i) ci_tilde = -1;
1604 if (ci_tilde_offd_mark |= i) ci_tilde_offd = -1;
1605 if (CF_marker_offd[i] == -1)
1606 {
1607 for (ji = S_ext_i[i]; ji < S_ext_i[i+1]; ji++)
1608 {
1609 j = S_ext_j[ji];
1610 if (j >= 0)
1611 {
1612 if (CF_marker[j] > 0)
1613 graph_array[j] = i;
1614 }
1615 else
1616 {
1617 jj = -j-1;
1618 if (CF_marker_offd[jj] > 0)
1619 ci_array[jj] = i;
1620 }
1621 }
1622 for (ji = S_ext_i[i]; ji < S_ext_i[i+1]; ji++)
1623 {
1624 j = S_ext_j[ji];
1625 if (j >= 0)
1626 {
1627 if ( CF_marker[j] == -1)
1628 {
1629 set_empty = 1;
1630 for (jj = S_i[j]; jj < S_i[j+1]; jj++)
1631 {
1632 index = S_j[jj];
1633 if (graph_array[index] == i)
1634 {
1635 set_empty = 0;
1636 break;
1637 }
1638 }
1639 for (jj = S_offd_i[j]; jj < S_offd_i[j+1]; jj++)
1640 {
1641 index = S_offd_j[jj];
1642 if (ci_array[index] == i)
1643 {
1644 set_empty = 0;
1645 break;
1646 }
1647 }
1648 if (set_empty)
1649 {
1650 if (C_i_nonempty)
1651 {
1652 CF_marker_offd[i] = 1;
1653 if (ci_tilde > -1)
1654 {
1655 CF_marker[ci_tilde] = -1;
1656 ci_tilde = -1;
1657 }
1658 if (ci_tilde_offd > -1)
1659 {
1660 CF_marker_offd[ci_tilde_offd] = -1;
1661 ci_tilde_offd = -1;
1662 }
1663 C_i_nonempty = 0;
1664 break;
1665 }
1666 else
1667 {
1668 ci_tilde = j;
1669 ci_tilde_mark = i;
1670 CF_marker[j] = 1;
1671 C_i_nonempty = 1;
1672 i--;
1673 break;
1674 }
1675 }
1676 }
1677 }
1678 else
1679 {
1680 jm = -j-1;
1681 if (CF_marker_offd[jm] == -1)
1682 {
1683 set_empty = 1;
1684 for (jj = S_ext_i[jm]; jj < S_ext_i[jm+1]; jj++)
1685 {
1686 index = S_ext_j[jj];
1687 if (index >= 0)
1688 {
1689 if (graph_array[index] == i)
1690 {
1691 set_empty = 0;
1692 break;
1693 }
1694 }
1695 else
1696 {
1697 jk = -index-1;
1698 if (ci_array[jk] == i)
1699 {
1700 set_empty = 0;
1701 break;
1702 }
1703 }
1704 }
1705 if (set_empty)
1706 {
1707 if (C_i_nonempty)
1708 {
1709 CF_marker_offd[i] = 1;
1710 if (ci_tilde > -1)
1711 {
1712 CF_marker[ci_tilde] = -1;
1713 ci_tilde = -1;
1714 }
1715 if (ci_tilde_offd > -1)
1716 {
1717 CF_marker_offd[ci_tilde_offd] = -1;
1718 ci_tilde_offd = -1;
1719 }
1720 C_i_nonempty = 0;
1721 break;
1722 }
1723 else
1724 {
1725 ci_tilde_offd = jm;
1726 ci_tilde_offd_mark = i;
1727 CF_marker_offd[jm] = 1;
1728 C_i_nonempty = 1;
1729 i--;
1730 break;
1731 }
1732 }
1733 }
1734 }
1735 }
1736 }
1737 }
1738 /*------------------------------------------------
1739 * Send boundary data for CF_marker back
1740 *------------------------------------------------*/
1741 if (num_procs > 1)
1742 {
1743 comm_handle = hypre_ParCSRCommHandleCreate(12, comm_pkg,
1744 CF_marker_offd, int_buf_data);
1745
1746 hypre_ParCSRCommHandleDestroy(comm_handle);
1747 }
1748
1749 /* only CF_marker entries from larger procs are accepted
1750 if coarsen_type = 4 coarse points are not overwritten */
1751
1752 index = 0;
1753 if (coarsen_type != 4)
1754 {
1755 for (i = 0; i < num_sends; i++)
1756 {
1757 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1758 if (hypre_ParCSRCommPkgSendProc(comm_pkg,i) > my_id)
1759 {
1760 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1761 CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)] =
1762 int_buf_data[index++];
1763 }
1764 else
1765 {
1766 index += hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1) - start;
1767 }
1768 }
1769 }
1770 else
1771 {
1772 for (i = 0; i < num_sends; i++)
1773 {
1774 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1775 if (hypre_ParCSRCommPkgSendProc(comm_pkg,i) > my_id)
1776 {
1777 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1778 {
1779 elmt = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
1780 if (CF_marker[elmt] != 1)
1781 CF_marker[elmt] = int_buf_data[index];
1782 index++;
1783 }
1784 }
1785 else
1786 {
1787 index += hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1) - start;
1788 }
1789 }
1790 }
1791 if (debug_flag == 3)
1792 {
1793 wall_time = time_getWallclockSeconds() - wall_time;
1794 if (coarsen_type == 4)
1795 printf("Proc = %d Coarsen 3rd pass = %f\n",
1796 my_id, wall_time);
1797 if (coarsen_type == 3)
1798 printf("Proc = %d Coarsen 3rd pass = %f\n",
1799 my_id, wall_time);
1800 if (coarsen_type == 2)
1801 printf("Proc = %d Coarsen 2nd pass = %f\n",
1802 my_id, wall_time);
1803 }
1804 }
1805 if (coarsen_type == 5)
1806 {
1807 /*------------------------------------------------
1808 * Exchange boundary data for CF_marker
1809 *------------------------------------------------*/
1810
1811 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
1812
1813 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
1814 if (int_buf_data == NULL)
1815 int_buf_data = hypre_CTAlloc(int,
1816 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1817
1818 index = 0;
1819 for (i = 0; i < num_sends; i++)
1820 {
1821 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1822 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1823 int_buf_data[index++]
1824 = CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
1825 }
1826
1827 if (num_procs > 1)
1828 {
1829 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
1830 int_buf_data, CF_marker_offd);
1831
1832 hypre_ParCSRCommHandleDestroy(comm_handle);
1833 }
1834
1835 ci_array = hypre_CTAlloc(int,num_cols_offd);
1836 for (i=0; i < num_cols_offd; i++)
1837 ci_array[i] = -1;
1838 for (i=0; i < num_variables; i++)
1839 graph_array[i] = -1;
1840
1841 for (i=0; i < num_variables; i++)
1842 {
1843 if (CF_marker[i] == -1 && (S_offd_i[i+1]-S_offd_i[i]) > 0)
1844 {
1845 break_var = 1;
1846 for (ji = S_i[i]; ji < S_i[i+1]; ji++)
1847 {
1848 j = S_j[ji];
1849 if (CF_marker[j] > 0)
1850 graph_array[j] = i;
1851 }
1852 for (ji = S_offd_i[i]; ji < S_offd_i[i+1]; ji++)
1853 {
1854 j = S_offd_j[ji];
1855 if (CF_marker_offd[j] > 0)
1856 ci_array[j] = i;
1857 }
1858 for (ji = S_offd_i[i]; ji < S_offd_i[i+1]; ji++)
1859 {
1860 j = S_offd_j[ji];
1861 if (CF_marker_offd[j] == -1)
1862 {
1863 set_empty = 1;
1864 for (jj = S_ext_i[j]; jj < S_ext_i[j+1]; jj++)
1865 {
1866 index = S_ext_j[jj];
1867 if (index >= 0) /* index interior */
1868 {
1869 if (graph_array[index] == i)
1870 {
1871 set_empty = 0;
1872 break;
1873 }
1874 }
1875 else
1876 {
1877 jk = -index-1;
1878 if (ci_array[jk] == i)
1879 {
1880 set_empty = 0;
1881 break;
1882 }
1883 }
1884 }
1885 if (set_empty)
1886 {
1887 if (C_i_nonempty)
1888 {
1889 CF_marker[i] = -2;
1890 C_i_nonempty = 0;
1891 break;
1892 }
1893 else
1894 {
1895 C_i_nonempty = 1;
1896 i--;
1897 break;
1898 }
1899 }
1900 }
1901 }
1902 }
1903 }
1904 if (debug_flag == 3)
1905 {
1906 wall_time = time_getWallclockSeconds() - wall_time;
1907 printf("Proc = %d Coarsen special points = %f\n",
1908 my_id, wall_time);
1909 }
1910
1911 }
1912 /*---------------------------------------------------
1913 * Clean up and return
1914 *---------------------------------------------------*/
1915
1916 hypre_TFree(int_buf_data);
1917 if (coarsen_type != 1)
1918 {
1919 hypre_TFree(CF_marker_offd);
1920 hypre_TFree(ci_array);
1921 }
1922 hypre_TFree(graph_array);
1923 if ((meas_type || coarsen_type != 1) && num_procs > 1)
1924 hypre_CSRMatrixDestroy(S_ext);
1925
1926 *CF_marker_ptr = CF_marker;
1927
1928 return (ierr);
1929}
1930
1931
1932int
1933hypre_BoomerAMGCoarsenFalgout( hypre_ParCSRMatrix *S,
1934 hypre_ParCSRMatrix *A,
1935 int measure_type,
1936 int debug_flag,
1937 int **CF_marker_ptr)
1938{
1939 int ierr = 0;
1940
1941 /*-------------------------------------------------------
1942 * Perform Ruge coarsening followed by CLJP coarsening
1943 *-------------------------------------------------------*/
1944
1945 ierr += hypre_BoomerAMGCoarsenRuge (S, A, measure_type, 6, debug_flag,
1946 CF_marker_ptr);
1947
1948 ierr += hypre_BoomerAMGCoarsen (S, A, 1, debug_flag,
1949 CF_marker_ptr);
1950
1951 return (ierr);
1952}
1953
1954int
1955hypre_BoomerAMGCoarsenHMIS( hypre_ParCSRMatrix *S,
1956 hypre_ParCSRMatrix *A,
1957 int measure_type,
1958 int debug_flag,
1959 int **CF_marker_ptr)
1960{
1961 int ierr = 0;
1962
1963 /*-------------------------------------------------------
1964 * Perform Ruge coarsening followed by CLJP coarsening
1965 *-------------------------------------------------------*/
1966
1967 ierr += hypre_BoomerAMGCoarsenRuge (S, A, measure_type, 10, debug_flag,
1968 CF_marker_ptr);
1969
1970 ierr += hypre_BoomerAMGCoarsenPMIS (S, A, 1, debug_flag,
1971 CF_marker_ptr);
1972
1973 return (ierr);
1974}
1975
1976/*--------------------------------------------------------------------------*/
1977
1978#define C_PT 1
1979#define F_PT -1
1980#define SF_PT -3
1981#define COMMON_C_PT 2
1982#define Z_PT -2
1983
1984 /* begin HANS added */
1985/**************************************************************
1986 *
1987 * Modified Independent Set Coarsening routine
1988 * (don't worry about strong F-F connections
1989 * without a common C point)
1990 *
1991 **************************************************************/
1992int
1993hypre_BoomerAMGCoarsenPMIS( hypre_ParCSRMatrix *S,
1994 hypre_ParCSRMatrix *A,
1995 int CF_init,
1996 int debug_flag,
1997 int **CF_marker_ptr)
1998{
1999 MPI_Comm comm = hypre_ParCSRMatrixComm(S);
2000 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(S);
2001 hypre_ParCSRCommHandle *comm_handle;
2002
2003 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
2004 int *S_diag_i = hypre_CSRMatrixI(S_diag);
2005 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
2006
2007 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
2008 int *S_offd_i = hypre_CSRMatrixI(S_offd);
2009 int *S_offd_j;
2010
2011 int num_variables = hypre_CSRMatrixNumRows(S_diag);
2012 int num_cols_offd = 0;
2013
2014 /* hypre_CSRMatrix *S_ext;
2015 int *S_ext_i;
2016 int *S_ext_j; */
2017
2018 int num_sends = 0;
2019 int *int_buf_data;
2020 double *buf_data;
2021
2022 int *CF_marker;
2023 int *CF_marker_offd;
2024
2025 double *measure_array;
2026 int *graph_array;
2027 int *graph_array_offd;
2028 int graph_size;
2029 HYPRE_BigInt big_graph_size;
2030 int graph_offd_size;
2031 HYPRE_BigInt global_graph_size;
2032
2033 int i, j, jS, ig;
2034 int index, start, my_id, num_procs, jrow, cnt, elmt;
2035
2036 int ierr = 0;
2037
2038 double wall_time;
2039 int iter = 0;
2040
2041
2042
2043#if 0 /* debugging */
2044 char filename[256];
2045 FILE *fp;
2046 int iter = 0;
2047#endif
2048
2049 /*******************************************************************************
2050 BEFORE THE INDEPENDENT SET COARSENING LOOP:
2051 measure_array: calculate the measures, and communicate them
2052 (this array contains measures for both local and external nodes)
2053 CF_marker, CF_marker_offd: initialize CF_marker
2054 (separate arrays for local and external; 0=unassigned, negative=F point, positive=C point)
2055 ******************************************************************************/
2056
2057 /*--------------------------------------------------------------
2058 * Use the ParCSR strength matrix, S.
2059 *
2060 * For now, the "strength" of dependence/influence is defined in
2061 * the following way: i depends on j if
2062 * aij > hypre_max (k != i) aik, aii < 0
2063 * or
2064 * aij < hypre_min (k != i) aik, aii >= 0
2065 * Then S_ij = 1, else S_ij = 0.
2066 *
2067 * NOTE: S_data is not used; in stead, only strong columns are retained
2068 * in S_j, which can then be used like S_data
2069 *----------------------------------------------------------------*/
2070
2071 /*S_ext = NULL; */
2072 if (debug_flag == 3) wall_time = time_getWallclockSeconds();
2073 MPI_Comm_size(comm,&num_procs);
2074 MPI_Comm_rank(comm,&my_id);
2075
2076 if (!comm_pkg)
2077 {
2078 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
2079 }
2080
2081 if (!comm_pkg)
2082 {
2083#ifdef HYPRE_NO_GLOBAL_PARTITION
2084 hypre_NewCommPkgCreate(A);
2085#else
2086 hypre_MatvecCommPkgCreate(A);
2087#endif
2088 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
2089 }
2090
2091 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
2092
2093 int_buf_data = hypre_CTAlloc(int, hypre_ParCSRCommPkgSendMapStart(comm_pkg,
2094 num_sends));
2095 buf_data = hypre_CTAlloc(double, hypre_ParCSRCommPkgSendMapStart(comm_pkg,
2096 num_sends));
2097
2098 num_cols_offd = hypre_CSRMatrixNumCols(S_offd);
2099
2100 S_diag_j = hypre_CSRMatrixJ(S_diag);
2101
2102 if (num_cols_offd)
2103 {
2104 S_offd_j = hypre_CSRMatrixJ(S_offd);
2105 }
2106
2107 /*----------------------------------------------------------
2108 * Compute the measures
2109 *
2110 * The measures are currently given by the column sums of S.
2111 * Hence, measure_array[i] is the number of influences
2112 * of variable i.
2113 *
2114 * The measures are augmented by a random number
2115 * between 0 and 1.
2116 *----------------------------------------------------------*/
2117
2118 measure_array = hypre_CTAlloc(double, num_variables+num_cols_offd);
2119
2120 /* first calculate the local part of the sums for the external nodes */
2121 for (i=0; i < S_offd_i[num_variables]; i++)
2122 {
2123 measure_array[num_variables + S_offd_j[i]] += 1.0;
2124 }
2125
2126 /* now send those locally calculated values for the external nodes to the neighboring processors */
2127 if (num_procs > 1)
2128 comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg,
2129 &measure_array[num_variables], buf_data);
2130
2131 /* calculate the local part for the local nodes */
2132 for (i=0; i < S_diag_i[num_variables]; i++)
2133 {
2134 measure_array[S_diag_j[i]] += 1.0;
2135 }
2136
2137 /* finish the communication */
2138 if (num_procs > 1)
2139 hypre_ParCSRCommHandleDestroy(comm_handle);
2140
2141 /* now add the externally calculated part of the local nodes to the local nodes */
2142 index = 0;
2143 for (i=0; i < num_sends; i++)
2144 {
2145 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
2146 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
2147 measure_array[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)]
2148 += buf_data[index++];
2149 }
2150
2151 /* set the measures of the external nodes to zero */
2152 for (i=num_variables; i < num_variables+num_cols_offd; i++)
2153 {
2154 measure_array[i] = 0;
2155 }
2156
2157 /* this augments the measures with a random number between 0 and 1 */
2158 /* (only for the local part) */
2159 /* this augments the measures */
2160 if (CF_init == 2 || CF_init == 4)
2161 hypre_BoomerAMGIndepSetInit(S, measure_array, 1);
2162 else
2163 hypre_BoomerAMGIndepSetInit(S, measure_array, 0);
2164
2165 /*---------------------------------------------------
2166 * Initialize the graph arrays, and CF_marker arrays
2167 *---------------------------------------------------*/
2168
2169 /* first the off-diagonal part of the graph array */
2170 if (num_cols_offd)
2171 graph_array_offd = hypre_CTAlloc(int, num_cols_offd);
2172 else
2173 graph_array_offd = NULL;
2174
2175 for (ig = 0; ig < num_cols_offd; ig++)
2176 graph_array_offd[ig] = ig;
2177
2178 graph_offd_size = num_cols_offd;
2179
2180 /* now the local part of the graph array, and the local CF_marker array */
2181 graph_array = hypre_CTAlloc(int, num_variables);
2182
2183 if (CF_init==1)
2184 {
2185 CF_marker = *CF_marker_ptr;
2186 cnt = 0;
2187 for (i=0; i < num_variables; i++)
2188 {
2189 if ( (S_offd_i[i+1]-S_offd_i[i]) > 0 || CF_marker[i] == -1)
2190 {
2191 CF_marker[i] = 0;
2192 }
2193 if ( CF_marker[i] == Z_PT)
2194 {
2195 if (measure_array[i] >= 1.0 ||
2196 (S_diag_i[i+1]-S_diag_i[i]) > 0)
2197 {
2198 CF_marker[i] = 0;
2199 graph_array[cnt++] = i;
2200 }
2201 else
2202 {
2203 CF_marker[i] = F_PT;
2204 }
2205 }
2206 else if (CF_marker[i] == SF_PT)
2207 measure_array[i] = 0;
2208 else
2209 graph_array[cnt++] = i;
2210 }
2211 }
2212 else
2213 {
2214 CF_marker = hypre_CTAlloc(int, num_variables);
2215 cnt = 0;
2216 for (i=0; i < num_variables; i++)
2217 {
2218 CF_marker[i] = 0;
2219 if ( (S_diag_i[i+1]-S_diag_i[i]) == 0
2220 && (S_offd_i[i+1]-S_offd_i[i]) == 0)
2221 {
2222 CF_marker[i] = SF_PT;
2223 if (CF_init == 3 || CF_init == 4) CF_marker[i] = C_PT;
2224 measure_array[i] = 0;
2225 }
2226 else
2227 graph_array[cnt++] = i;
2228 }
2229 }
2230 graph_size = cnt;
2231
2232 /* now the off-diagonal part of CF_marker */
2233 if (num_cols_offd)
2234 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
2235 else
2236 CF_marker_offd = NULL;
2237
2238 for (i=0; i < num_cols_offd; i++)
2239 CF_marker_offd[i] = 0;
2240
2241 /*------------------------------------------------
2242 * Communicate the local measures, which are complete,
2243 to the external nodes
2244 *------------------------------------------------*/
2245 index = 0;
2246 for (i = 0; i < num_sends; i++)
2247 {
2248 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
2249 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
2250 {
2251 jrow = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
2252 buf_data[index++] = measure_array[jrow];
2253 }
2254 }
2255
2256 if (num_procs > 1)
2257 {
2258 comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg, buf_data,
2259 &measure_array[num_variables]);
2260
2261 hypre_ParCSRCommHandleDestroy(comm_handle);
2262
2263 }
2264
2265 /* we need S_ext: the columns of the S matrix for the local nodes */
2266 /* we need this because the independent set routine can only decide
2267 which local nodes are in it when it knows both the rows and columns
2268 of S */
2269
2270 /* if (num_procs > 1)
2271 {
2272 S_ext = hypre_ParCSRMatrixExtractBExt(S,A,0);
2273 S_ext_i = hypre_CSRMatrixI(S_ext);
2274 S_ext_j = hypre_CSRMatrixJ(S_ext);
2275 } */
2276
2277 /* compress S_ext and convert column numbers*/
2278
2279 /* index = 0;
2280 for (i=0; i < num_cols_offd; i++)
2281 {
2282 for (j=S_ext_i[i]; j < S_ext_i[i+1]; j++)
2283 {
2284 k = S_ext_j[j];
2285 if (k >= col_1 && k < col_n)
2286 {
2287 S_ext_j[index++] = k - col_1;
2288 }
2289 else
2290 {
2291 kc = hypre_BinarySearch(col_map_offd,k,num_cols_offd);
2292 if (kc > -1) S_ext_j[index++] = -kc-1;
2293 }
2294 }
2295 S_ext_i[i] = index;
2296 }
2297 for (i = num_cols_offd; i > 0; i--)
2298 S_ext_i[i] = S_ext_i[i-1];
2299 if (num_procs > 1) S_ext_i[0] = 0; */
2300
2301 if (debug_flag == 3)
2302 {
2303 wall_time = time_getWallclockSeconds() - wall_time;
2304 printf("Proc = %d Initialize CLJP phase = %f\n",
2305 my_id, wall_time);
2306 }
2307
2308 /*******************************************************************************
2309 THE INDEPENDENT SET COARSENING LOOP:
2310 ******************************************************************************/
2311
2312 /*---------------------------------------------------
2313 * Loop until all points are either fine or coarse.
2314 *---------------------------------------------------*/
2315
2316 while (1)
2317 {
2318
2319 big_graph_size = (HYPRE_BigInt) graph_size;
2320 /* stop the coarsening if nothing left to be coarsened */
2321 MPI_Allreduce(&big_graph_size,&global_graph_size,1,MPI_HYPRE_BIG_INT,MPI_SUM,comm);
2322
2323 if (global_graph_size == 0)
2324 break;
2325
2326 /* printf("\n");
2327 printf("*** MIS iteration %d\n",iter);
2328 printf("graph_size remaining %d\n",graph_size);*/
2329
2330 /*------------------------------------------------
2331 * Pick an independent set of points with
2332 * maximal measure.
2333 At the end, CF_marker is complete, but still needs to be
2334 communicated to CF_marker_offd
2335 *------------------------------------------------*/
2336 if (!CF_init || iter)
2337 {
2338 hypre_BoomerAMGIndepSet(S, measure_array, graph_array,
2339 graph_size,
2340 graph_array_offd, graph_offd_size,
2341 CF_marker, CF_marker_offd);
2342
2343 /*------------------------------------------------
2344 * Exchange boundary data for CF_marker: send internal
2345 points to external points
2346 *------------------------------------------------*/
2347
2348 if (num_procs > 1)
2349 {
2350 comm_handle = hypre_ParCSRCommHandleCreate(12, comm_pkg,
2351 CF_marker_offd, int_buf_data);
2352
2353 hypre_ParCSRCommHandleDestroy(comm_handle);
2354 }
2355
2356 index = 0;
2357 for (i = 0; i < num_sends; i++)
2358 {
2359 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
2360 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
2361 {
2362 elmt = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
2363 if (!int_buf_data[index] && CF_marker[elmt] > 0)
2364 {
2365 CF_marker[elmt] = 0;
2366 index++;
2367 }
2368 else
2369 {
2370 int_buf_data[index++] = CF_marker[elmt];
2371 }
2372 }
2373 }
2374
2375 if (num_procs > 1)
2376 {
2377 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
2378 CF_marker_offd);
2379
2380 hypre_ParCSRCommHandleDestroy(comm_handle);
2381 }
2382 }
2383
2384 iter++;
2385 /*------------------------------------------------
2386 * Set C-pts and F-pts.
2387 *------------------------------------------------*/
2388
2389 for (ig = 0; ig < graph_size; ig++)
2390 {
2391 i = graph_array[ig];
2392
2393 /*---------------------------------------------
2394 * If the measure of i is smaller than 1, then
2395 * make i and F point (because it does not influence
2396 * any other point), and remove all edges of
2397 * equation i.
2398 *---------------------------------------------*/
2399
2400 if(measure_array[i]<1.)
2401 {
2402 /* make point i an F point*/
2403 CF_marker[i]= F_PT;
2404
2405 /* remove the edges in equation i */
2406 /* first the local part */
2407 /*for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++) {
2408 j = S_diag_j[jS];
2409 if (j > -1){
2410 S_diag_j[jS] = -S_diag_j[jS]-1;
2411 }
2412 }*/
2413 /* now the external part */
2414 /*for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++) {
2415 j = S_offd_j[jS];
2416 if (j > -1){
2417 S_offd_j[jS] = -S_offd_j[jS]-1;
2418 }
2419 }*/
2420 }
2421
2422 /*---------------------------------------------
2423 * First treat the case where point i is in the
2424 * independent set: make i a C point,
2425 * take out all the graph edges for
2426 * equation i.
2427 *---------------------------------------------*/
2428
2429 if (CF_marker[i] > 0)
2430 {
2431 /* set to be a C-pt */
2432 CF_marker[i] = C_PT;
2433
2434 /* remove the edges in equation i */
2435 /* first the local part */
2436 /*for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++) {
2437 j = S_diag_j[jS];
2438 if (j > -1){
2439 S_diag_j[jS] = -S_diag_j[jS]-1;
2440 }
2441 }*/
2442 /* now the external part */
2443 /*for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++) {
2444 j = S_offd_j[jS];
2445 if (j > -1){
2446 S_offd_j[jS] = -S_offd_j[jS]-1;
2447 }
2448 }*/
2449 }
2450
2451 /*---------------------------------------------
2452 * Now treat the case where point i is not in the
2453 * independent set: loop over
2454 * all the points j that influence equation i; if
2455 * j is a C point, then make i an F point.
2456 * If i is a new F point, then remove all the edges
2457 * from the graph for equation i.
2458 *---------------------------------------------*/
2459
2460 else
2461 {
2462
2463 /* first the local part */
2464 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++)
2465 {
2466 /* j is the column number, or the local number of the point influencing i */
2467 j = S_diag_j[jS];
2468 /*if(j<0) j=-j-1;*/
2469
2470 if (CF_marker[j] > 0)
2471 { /* j is a C-point */
2472 CF_marker[i] = F_PT;
2473 }
2474 }
2475 /* now the external part */
2476 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++)
2477 {
2478 j = S_offd_j[jS];
2479 /*if(j<0) j=-j-1;*/
2480 if (CF_marker_offd[j] > 0)
2481 { /* j is a C-point */
2482 CF_marker[i] = F_PT;
2483 }
2484 }
2485
2486 /* remove all the edges for equation i if i is a new F point */
2487 /*if (CF_marker[i] == F_PT){
2488 for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++) {
2489 j = S_diag_j[jS];
2490 if (j > -1){
2491 S_diag_j[jS] = -S_diag_j[jS]-1;
2492 }
2493 }
2494 for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++) {
2495 j = S_offd_j[jS];
2496 if (j > -1){
2497 S_offd_j[jS] = -S_offd_j[jS]-1;
2498 }
2499 }
2500 } */
2501 } /* end else */
2502 } /* end first loop over graph */
2503
2504 /* now communicate CF_marker to CF_marker_offd, to make
2505 sure that new external F points are known on this processor */
2506
2507 /*------------------------------------------------
2508 * Exchange boundary data for CF_marker: send internal
2509 points to external points
2510 *------------------------------------------------*/
2511
2512 index = 0;
2513 for (i = 0; i < num_sends; i++)
2514 {
2515 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
2516 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
2517 int_buf_data[index++]
2518 = CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
2519 }
2520
2521 if (num_procs > 1)
2522 {
2523 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
2524 CF_marker_offd);
2525
2526 hypre_ParCSRCommHandleDestroy(comm_handle);
2527 }
2528
2529 /*---------------------------------------------
2530 * Now loop over the points i in the unassigned
2531 * graph again. For all points i that are no new C or
2532 * F points, remove the edges in equation i that
2533 * connect to C or F points.
2534 * (We have removed the rows for the new C and F
2535 * points above; now remove the columns.)
2536 *---------------------------------------------*/
2537
2538 /*for (ig = 0; ig < graph_size; ig++) {
2539 i = graph_array[ig];
2540
2541 if(CF_marker[i]==0) {*/
2542
2543 /* first the local part */
2544 /*for (jS = S_diag_i[i]; jS < S_diag_i[i+1]; jS++) {
2545 j = S_diag_j[jS];
2546 if(j<0) j=-j-1;
2547
2548 if (!CF_marker[j]==0 && S_diag_j[jS] > -1){ */ /* connection to C or F point, and
2549 column number is still positive; not accounted for yet */
2550 /*S_diag_j[jS] = -S_diag_j[jS]-1;
2551 }
2552 }*/
2553 /* now the external part */
2554 /*for (jS = S_offd_i[i]; jS < S_offd_i[i+1]; jS++) {
2555 j = S_offd_j[jS];
2556 if(j<0) j=-j-1;
2557
2558 if (!CF_marker_offd[j]==0 && S_offd_j[jS] > -1){*/ /* connection to C or F point, and
2559 column number is still positive; not accounted for yet */
2560 /*S_offd_j[jS] = -S_offd_j[jS]-1;
2561 }
2562 }
2563 }
2564 } *//* end second loop over graph */
2565
2566 /*------------------------------------------------
2567 * Update subgraph
2568 *------------------------------------------------*/
2569
2570 for (ig = 0; ig < graph_size; ig++)
2571 {
2572 i = graph_array[ig];
2573
2574 if (!CF_marker[i]==0) /* C or F point */
2575 {
2576 /* the independent set subroutine needs measure 0 for
2577 removed nodes */
2578 measure_array[i] = 0;
2579 /* take point out of the subgraph */
2580 graph_size--;
2581 graph_array[ig] = graph_array[graph_size];
2582 graph_array[graph_size] = i;
2583 ig--;
2584 }
2585 }
2586 for (ig = 0; ig < graph_offd_size; ig++)
2587 {
2588 i = graph_array_offd[ig];
2589
2590 if (!CF_marker_offd[i]==0) /* C or F point */
2591 {
2592 /* the independent set subroutine needs measure 0 for
2593 removed nodes */
2594 measure_array[i+num_variables] = 0;
2595 /* take point out of the subgraph */
2596 graph_offd_size--;
2597 graph_array_offd[ig] = graph_array_offd[graph_offd_size];
2598 graph_array_offd[graph_offd_size] = i;
2599 ig--;
2600 }
2601 }
2602
2603 } /* end while */
2604
2605 /* printf("*** MIS iteration %d\n",iter);
2606 printf("graph_size remaining %d\n",graph_size);
2607
2608 printf("num_cols_offd %d\n",num_cols_offd);
2609 for (i=0;i<num_variables;i++)
2610 {
2611 if(CF_marker[i]==1)
2612 printf("node %d CF %d\n",i,CF_marker[i]);
2613 }*/
2614
2615
2616 /*---------------------------------------------------
2617 * Clean up and return
2618 *---------------------------------------------------*/
2619
2620 /* Reset S_matrix */
2621 /*for (i=0; i < S_diag_i[num_variables]; i++)
2622 {
2623 if (S_diag_j[i] < 0)
2624 S_diag_j[i] = -S_diag_j[i]-1;
2625 }
2626 for (i=0; i < S_offd_i[num_variables]; i++)
2627 {
2628 if (S_offd_j[i] < 0)
2629 S_offd_j[i] = -S_offd_j[i]-1;
2630 }*/
2631 /*for (i=0; i < num_variables; i++)
2632 if (CF_marker[i] == SF_PT) CF_marker[i] = F_PT;*/
2633
2634 hypre_TFree(measure_array);
2635 hypre_TFree(graph_array);
2636 if (num_cols_offd) hypre_TFree(graph_array_offd);
2637 hypre_TFree(buf_data);
2638 hypre_TFree(int_buf_data);
2639 hypre_TFree(CF_marker_offd);
2640 /*if (num_procs > 1) hypre_CSRMatrixDestroy(S_ext);*/
2641
2642 *CF_marker_ptr = CF_marker;
2643
2644 return (ierr);
2645}
Note: See TracBrowser for help on using the repository browser.