source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_strength.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: 60.0 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
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 Generates strength 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 coarsening and interpolation routines.
38 \end{itemize}
39 \item The graph of the "strength matrix" for A is a subgraph of the
40 graph of A, but requires nonsymmetric storage even if A is
41 symmetric. This is because of the directional nature of the
42 "strengh of dependence" notion (see below). Since we are using
43 nonsymmetric storage for A right now, this is not a problem. If we
44 ever add the ability to store A symmetrically, then we could store
45 the strength graph as floats instead of doubles to save space.
46 \item This routine currently "compresses" the strength matrix. We
47 should consider the possibility of defining this matrix to have the
48 same "nonzero structure" as A. To do this, we could use the same
49 A\_i and A\_j arrays, and would need only define the S\_data array.
50 There are several pros and cons to discuss.
51 \end{itemize}
52
53 Terminology:
54 \begin{itemize}
55 \item Ruge's terminology: A point is "strongly connected to" $j$, or
56 "strongly depends on" $j$, if $-a_ij >= \theta max_{l != j} \{-a_il\}$.
57 \item Here, we retain some of this terminology, but with a more
58 generalized notion of "strength". We also retain the "natural"
59 graph notation for representing the directed graph of a matrix.
60 That is, the nonzero entry $a_ij$ is represented as: i --> j. In
61 the strength matrix, S, the entry $s_ij$ is also graphically denoted
62 as above, and means both of the following:
63 \begin{itemize}
64 \item $i$ "depends on" $j$ with "strength" $s_ij$
65 \item $j$ "influences" $i$ with "strength" $s_ij$
66 \end{itemize}
67 \end{itemize}
68
69 {\bf Input files:}
70 headers.h
71
72 @return Error code.
73
74 @param A [IN]
75 coefficient matrix
76 @param strength_threshold [IN]
77 threshold parameter used to define strength
78 @param max_row_sum [IN]
79 parameter used to modify definition of strength for diagonal dominant matrices
80 @param S_ptr [OUT]
81 strength matrix
82
83 @see */
84/*--------------------------------------------------------------------------*/
85
86int
87hypre_BoomerAMGCreateS(hypre_ParCSRMatrix *A,
88 double strength_threshold,
89 double max_row_sum,
90 int num_functions,
91 int *dof_func,
92 hypre_ParCSRMatrix **S_ptr)
93{
94 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
95 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
96 hypre_ParCSRCommHandle *comm_handle;
97 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
98 int *A_diag_i = hypre_CSRMatrixI(A_diag);
99 double *A_diag_data = hypre_CSRMatrixData(A_diag);
100
101
102 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
103 int *A_offd_i = hypre_CSRMatrixI(A_offd);
104 double *A_offd_data;
105 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
106 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
107
108 HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(A);
109 int num_variables = hypre_CSRMatrixNumRows(A_diag);
110 HYPRE_BigInt global_num_vars = hypre_ParCSRMatrixGlobalNumRows(A);
111 int num_nonzeros_diag;
112 int num_nonzeros_offd = 0;
113 int num_cols_offd = 0;
114
115 hypre_ParCSRMatrix *S;
116 hypre_CSRMatrix *S_diag;
117 int *S_diag_i;
118 int *S_diag_j;
119 /* double *S_diag_data; */
120 hypre_CSRMatrix *S_offd;
121 int *S_offd_i;
122 int *S_offd_j;
123 /* double *S_offd_data; */
124
125 double diag, row_scale, row_sum;
126 int i, jA, jS;
127
128 int ierr = 0;
129
130 int *dof_func_offd;
131 int num_sends;
132 int *int_buf_data;
133 int index, start, j;
134
135 /*--------------------------------------------------------------
136 * Compute a ParCSR strength matrix, S.
137 *
138 * For now, the "strength" of dependence/influence is defined in
139 * the following way: i depends on j if
140 * aij > hypre_max (k != i) aik, aii < 0
141 * or
142 * aij < hypre_min (k != i) aik, aii >= 0
143 * Then S_ij = 1, else S_ij = 0.
144 *
145 * NOTE: the entries are negative initially, corresponding
146 * to "unaccounted-for" dependence.
147 *----------------------------------------------------------------*/
148
149 num_nonzeros_diag = A_diag_i[num_variables];
150 num_cols_offd = hypre_CSRMatrixNumCols(A_offd);
151
152 A_offd_i = hypre_CSRMatrixI(A_offd);
153 num_nonzeros_offd = A_offd_i[num_variables];
154
155 S = hypre_ParCSRMatrixCreate(comm, global_num_vars, global_num_vars,
156 row_starts, row_starts,
157 num_cols_offd, num_nonzeros_diag, num_nonzeros_offd);
158/* row_starts is owned by A, col_starts = row_starts */
159 hypre_ParCSRMatrixSetRowStartsOwner(S,0);
160 S_diag = hypre_ParCSRMatrixDiag(S);
161
162 hypre_CSRMatrixI(S_diag) = hypre_CTAlloc(int, num_variables+1);
163 hypre_CSRMatrixJ(S_diag) = hypre_CTAlloc(int, num_nonzeros_diag);
164 S_offd = hypre_ParCSRMatrixOffd(S);
165 hypre_CSRMatrixI(S_offd) = hypre_CTAlloc(int, num_variables+1);
166
167 S_diag_i = hypre_CSRMatrixI(S_diag);
168 S_diag_j = hypre_CSRMatrixJ(S_diag);
169 S_offd_i = hypre_CSRMatrixI(S_offd);
170
171 dof_func_offd = NULL;
172
173 if (num_cols_offd)
174 {
175 A_offd_data = hypre_CSRMatrixData(A_offd);
176 hypre_CSRMatrixJ(S_offd) = hypre_CTAlloc(int, num_nonzeros_offd);
177
178 S_offd_j = hypre_CSRMatrixJ(S_offd);
179 hypre_ParCSRMatrixColMapOffd(S) = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd);
180
181 if (num_functions > 1)
182 dof_func_offd = hypre_CTAlloc(int, num_cols_offd);
183 }
184
185
186 /*-------------------------------------------------------------------
187 * Get the dof_func data for the off-processor columns
188 *-------------------------------------------------------------------*/
189
190 if (!comm_pkg)
191 {
192#ifdef HYPRE_NO_GLOBAL_PARTITION
193 hypre_NewCommPkgCreate(A);
194#else
195 hypre_MatvecCommPkgCreate(A);
196#endif
197 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
198 }
199
200 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
201 if (num_functions > 1)
202 {
203 int_buf_data = hypre_CTAlloc(int,hypre_ParCSRCommPkgSendMapStart(comm_pkg,
204 num_sends));
205 index = 0;
206 for (i = 0; i < num_sends; i++)
207 {
208 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
209 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
210 int_buf_data[index++]
211 = dof_func[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
212 }
213
214 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data,
215 dof_func_offd);
216
217 hypre_ParCSRCommHandleDestroy(comm_handle);
218 hypre_TFree(int_buf_data);
219 }
220
221 /* give S same nonzero structure as A */
222 hypre_ParCSRMatrixCopy(A,S,0);
223
224#define HYPRE_SMP_PRIVATE i,diag,row_scale,row_sum,jA
225#include "../utilities/hypre_smp_forloop.h"
226 for (i = 0; i < num_variables; i++)
227 {
228 diag = A_diag_data[A_diag_i[i]];
229
230 /* compute scaling factor and row sum */
231 row_scale = 0.0;
232 row_sum = diag;
233 if (num_functions > 1)
234 {
235 if (diag < 0)
236 {
237 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
238 {
239 if (dof_func[i] == dof_func[A_diag_j[jA]])
240 {
241 row_scale = hypre_max(row_scale, A_diag_data[jA]);
242 row_sum += A_diag_data[jA];
243 }
244 }
245 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
246 {
247 if (dof_func[i] == dof_func_offd[A_offd_j[jA]])
248 {
249 row_scale = hypre_max(row_scale, A_offd_data[jA]);
250 row_sum += A_offd_data[jA];
251 }
252 }
253 }
254 else
255 {
256 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
257 {
258 if (dof_func[i] == dof_func[A_diag_j[jA]])
259 {
260 row_scale = hypre_min(row_scale, A_diag_data[jA]);
261 row_sum += A_diag_data[jA];
262 }
263 }
264 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
265 {
266 if (dof_func[i] == dof_func_offd[A_offd_j[jA]])
267 {
268 row_scale = hypre_min(row_scale, A_offd_data[jA]);
269 row_sum += A_offd_data[jA];
270 }
271 }
272 }
273 }
274 else
275 {
276 if (diag < 0)
277 {
278 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
279 {
280 row_scale = hypre_max(row_scale, A_diag_data[jA]);
281 row_sum += A_diag_data[jA];
282 }
283 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
284 {
285 row_scale = hypre_max(row_scale, A_offd_data[jA]);
286 row_sum += A_offd_data[jA];
287 }
288 }
289 else
290 {
291 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
292 {
293 row_scale = hypre_min(row_scale, A_diag_data[jA]);
294 row_sum += A_diag_data[jA];
295 }
296 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
297 {
298 row_scale = hypre_min(row_scale, A_offd_data[jA]);
299 row_sum += A_offd_data[jA];
300 }
301 }
302 }
303 row_sum = fabs( row_sum / diag );
304
305 /* compute row entries of S */
306 S_diag_j[A_diag_i[i]] = -1;
307 if ((row_sum > max_row_sum) && (max_row_sum < 1.0))
308 {
309 /* make all dependencies weak */
310 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
311 {
312 S_diag_j[jA] = -1;
313 }
314 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
315 {
316 S_offd_j[jA] = -1;
317 }
318 }
319 else
320 {
321 if (num_functions > 1)
322 {
323 if (diag < 0)
324 {
325 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
326 {
327 if (A_diag_data[jA] <= strength_threshold * row_scale
328 || dof_func[i] != dof_func[A_diag_j[jA]])
329 {
330 S_diag_j[jA] = -1;
331 }
332 }
333 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
334 {
335 if (A_offd_data[jA] <= strength_threshold * row_scale
336 || dof_func[i] != dof_func_offd[A_offd_j[jA]])
337 {
338 S_offd_j[jA] = -1;
339 }
340 }
341 }
342 else
343 {
344 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
345 {
346 if (A_diag_data[jA] >= strength_threshold * row_scale
347 || dof_func[i] != dof_func[A_diag_j[jA]])
348 {
349 S_diag_j[jA] = -1;
350 }
351 }
352 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
353 {
354 if (A_offd_data[jA] >= strength_threshold * row_scale
355 || dof_func[i] != dof_func_offd[A_offd_j[jA]])
356 {
357 S_offd_j[jA] = -1;
358 }
359 }
360 }
361 }
362 else
363 {
364 if (diag < 0)
365 {
366 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
367 {
368 if (A_diag_data[jA] <= strength_threshold * row_scale)
369 {
370 S_diag_j[jA] = -1;
371 }
372 }
373 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
374 {
375 if (A_offd_data[jA] <= strength_threshold * row_scale)
376 {
377 S_offd_j[jA] = -1;
378 }
379 }
380 }
381 else
382 {
383 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
384 {
385 if (A_diag_data[jA] >= strength_threshold * row_scale)
386 {
387 S_diag_j[jA] = -1;
388 }
389 }
390 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
391 {
392 if (A_offd_data[jA] >= strength_threshold * row_scale)
393 {
394 S_offd_j[jA] = -1;
395 }
396 }
397 }
398 }
399 }
400 }
401
402 /*--------------------------------------------------------------
403 * "Compress" the strength matrix.
404 *
405 * NOTE: S has *NO DIAGONAL ELEMENT* on any row. Caveat Emptor!
406 *
407 * NOTE: This "compression" section of code may be removed, and
408 * coarsening will still be done correctly. However, the routine
409 * that builds interpolation would have to be modified first.
410 *----------------------------------------------------------------*/
411
412/* RDF: not sure if able to thread this loop */
413 jS = 0;
414 for (i = 0; i < num_variables; i++)
415 {
416 S_diag_i[i] = jS;
417 for (jA = A_diag_i[i]; jA < A_diag_i[i+1]; jA++)
418 {
419 if (S_diag_j[jA] > -1)
420 {
421 S_diag_j[jS] = S_diag_j[jA];
422 jS++;
423 }
424 }
425 }
426 S_diag_i[num_variables] = jS;
427 hypre_CSRMatrixNumNonzeros(S_diag) = jS;
428
429/* RDF: not sure if able to thread this loop */
430 jS = 0;
431 for (i = 0; i < num_variables; i++)
432 {
433 S_offd_i[i] = jS;
434 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
435 {
436 if (S_offd_j[jA] > -1)
437 {
438 S_offd_j[jS] = S_offd_j[jA];
439 jS++;
440 }
441 }
442 }
443 S_offd_i[num_variables] = jS;
444 hypre_CSRMatrixNumNonzeros(S_offd) = jS;
445 hypre_ParCSRMatrixCommPkg(S) = NULL;
446
447 *S_ptr = S;
448
449 hypre_TFree(dof_func_offd);
450
451 return (ierr);
452}
453
454
455/*==========================================================================*/
456/*==========================================================================*/
457/**
458 Generates strength matrix
459
460 Notes:
461 \begin{itemize}
462 \item The underlying matrix storage scheme is a hypre_ParCSR matrix.
463 \item The routine returns the following:
464 \begin{itemize}
465 \item S - a ParCSR matrix representing the "strength matrix". This is
466 used in the coarsening and interpolation routines.
467 \end{itemize}
468 \item The graph of the "strength matrix" for A is a subgraph of the
469 graph of A, but requires nonsymmetric storage even if A is
470 symmetric. This is because of the directional nature of the
471 "strengh of dependence" notion (see below). Since we are using
472 nonsymmetric storage for A right now, this is not a problem. If we
473 ever add the ability to store A symmetrically, then we could store
474 the strength graph as floats instead of doubles to save space.
475 \item This routine currently "compresses" the strength matrix. We
476 should consider the possibility of defining this matrix to have the
477 same "nonzero structure" as A. To do this, we could use the same
478 A\_i and A\_j arrays, and would need only define the S\_data array.
479 There are several pros and cons to discuss.
480 \end{itemize}
481
482 Terminology:
483 \begin{itemize}
484 \item Ruge's terminology: A point is "strongly connected to" $j$, or
485 "strongly depends on" $j$, if $|a_ij| >= \theta max_{l != j} |a_il|}$.
486 \item Here, we retain some of this terminology, but with a more
487 generalized notion of "strength". We also retain the "natural"
488 graph notation for representing the directed graph of a matrix.
489 That is, the nonzero entry $a_ij$ is represented as: i --> j. In
490 the strength matrix, S, the entry $s_ij$ is also graphically denoted
491 as above, and means both of the following:
492 \begin{itemize}
493 \item $i$ "depends on" $j$ with "strength" $s_ij$
494 \item $j$ "influences" $i$ with "strength" $s_ij$
495 \end{itemize}
496 \end{itemize}
497
498 {\bf Input files:}
499 headers.h
500
501 @return Error code.
502
503 @param A [IN]
504 coefficient matrix
505 @param strength_threshold [IN]
506 threshold parameter used to define strength
507 @param max_row_sum [IN]
508 parameter used to modify definition of strength for diagonal dominant matrices
509 @param S_ptr [OUT]
510 strength matrix
511
512 @see */
513/*--------------------------------------------------------------------------*/
514
515int
516hypre_BoomerAMGCreateSabs(hypre_ParCSRMatrix *A,
517 double strength_threshold,
518 double max_row_sum,
519 int num_functions,
520 int *dof_func,
521 hypre_ParCSRMatrix **S_ptr)
522{
523 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
524 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
525 hypre_ParCSRCommHandle *comm_handle;
526 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
527 int *A_diag_i = hypre_CSRMatrixI(A_diag);
528 double *A_diag_data = hypre_CSRMatrixData(A_diag);
529
530
531 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
532 int *A_offd_i = hypre_CSRMatrixI(A_offd);
533 double *A_offd_data;
534 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
535 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
536
537 HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(A);
538 int num_variables = hypre_CSRMatrixNumRows(A_diag);
539 HYPRE_BigInt global_num_vars = hypre_ParCSRMatrixGlobalNumRows(A);
540 int num_nonzeros_diag;
541 int num_nonzeros_offd = 0;
542 int num_cols_offd = 0;
543
544 hypre_ParCSRMatrix *S;
545 hypre_CSRMatrix *S_diag;
546 int *S_diag_i;
547 int *S_diag_j;
548 /* double *S_diag_data; */
549 hypre_CSRMatrix *S_offd;
550 int *S_offd_i;
551 int *S_offd_j;
552 /* double *S_offd_data; */
553
554 double diag, row_scale, row_sum;
555 int i, jA, jS;
556
557 int ierr = 0;
558
559 int *dof_func_offd;
560 int num_sends;
561 int *int_buf_data;
562 int index, start, j;
563
564 /*--------------------------------------------------------------
565 * Compute a ParCSR strength matrix, S.
566 *
567 * For now, the "strength" of dependence/influence is defined in
568 * the following way: i depends on j if
569 * aij > hypre_max (k != i) aik, aii < 0
570 * or
571 * aij < hypre_min (k != i) aik, aii >= 0
572 * Then S_ij = 1, else S_ij = 0.
573 *
574 * NOTE: the entries are negative initially, corresponding
575 * to "unaccounted-for" dependence.
576 *----------------------------------------------------------------*/
577
578 num_nonzeros_diag = A_diag_i[num_variables];
579 num_cols_offd = hypre_CSRMatrixNumCols(A_offd);
580
581 A_offd_i = hypre_CSRMatrixI(A_offd);
582 num_nonzeros_offd = A_offd_i[num_variables];
583
584 S = hypre_ParCSRMatrixCreate(comm, global_num_vars, global_num_vars,
585 row_starts, row_starts,
586 num_cols_offd, num_nonzeros_diag, num_nonzeros_offd);
587/* row_starts is owned by A, col_starts = row_starts */
588 hypre_ParCSRMatrixSetRowStartsOwner(S,0);
589 S_diag = hypre_ParCSRMatrixDiag(S);
590
591 hypre_CSRMatrixI(S_diag) = hypre_CTAlloc(int, num_variables+1);
592 hypre_CSRMatrixJ(S_diag) = hypre_CTAlloc(int, num_nonzeros_diag);
593 S_offd = hypre_ParCSRMatrixOffd(S);
594 hypre_CSRMatrixI(S_offd) = hypre_CTAlloc(int, num_variables+1);
595 S_diag_i = hypre_CSRMatrixI(S_diag);
596 S_diag_j = hypre_CSRMatrixJ(S_diag);
597 S_offd_i = hypre_CSRMatrixI(S_offd);
598
599 dof_func_offd = NULL;
600
601 if (num_cols_offd)
602 {
603 A_offd_data = hypre_CSRMatrixData(A_offd);
604 hypre_CSRMatrixJ(S_offd) = hypre_CTAlloc(int, num_nonzeros_offd);
605 S_offd_j = hypre_CSRMatrixJ(S_offd);
606 hypre_ParCSRMatrixColMapOffd(S) = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd);
607 if (num_functions > 1)
608 dof_func_offd = hypre_CTAlloc(int, num_cols_offd);
609 }
610
611
612 /*-------------------------------------------------------------------
613 * Get the dof_func data for the off-processor columns
614 *-------------------------------------------------------------------*/
615
616 if (!comm_pkg)
617 {
618#ifdef HYPRE_NO_GLOBAL_PARTITION
619 hypre_NewCommPkgCreate(A);
620#else
621 hypre_MatvecCommPkgCreate(A);
622#endif
623 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
624 }
625
626 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
627 if (num_functions > 1)
628 {
629 int_buf_data = hypre_CTAlloc(int,hypre_ParCSRCommPkgSendMapStart(comm_pkg,
630 num_sends));
631 index = 0;
632 for (i = 0; i < num_sends; i++)
633 {
634 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
635 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
636 int_buf_data[index++]
637 = dof_func[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
638 }
639
640 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data,
641 dof_func_offd);
642
643 hypre_ParCSRCommHandleDestroy(comm_handle);
644 hypre_TFree(int_buf_data);
645 }
646
647 /* give S same nonzero structure as A */
648 hypre_ParCSRMatrixCopy(A,S,0);
649
650#define HYPRE_SMP_PRIVATE i,diag,row_scale,row_sum,jA
651#include "../utilities/hypre_smp_forloop.h"
652 for (i = 0; i < num_variables; i++)
653 {
654 diag = A_diag_data[A_diag_i[i]];
655
656 /* compute scaling factor and row sum */
657 row_scale = 0.0;
658 row_sum = diag;
659 if (num_functions > 1)
660 {
661 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
662 {
663 if (dof_func[i] == dof_func[A_diag_j[jA]])
664 {
665 row_scale = hypre_max(row_scale, fabs(A_diag_data[jA]));
666 row_sum += fabs(A_diag_data[jA]);
667 }
668 }
669 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
670 {
671 if (dof_func[i] == dof_func_offd[A_offd_j[jA]])
672 {
673 row_scale = hypre_max(row_scale, fabs(A_offd_data[jA]));
674 row_sum += fabs(A_offd_data[jA]);
675 }
676 }
677 }
678 else
679 {
680 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
681 {
682 row_scale = hypre_max(row_scale, fabs(A_diag_data[jA]));
683 row_sum += fabs(A_diag_data[jA]);
684 }
685 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
686 {
687 row_scale = hypre_max(row_scale, fabs(A_offd_data[jA]));
688 row_sum += fabs(A_offd_data[jA]);
689 }
690 }
691 row_sum = fabs( row_sum / diag );
692
693 /* compute row entries of S */
694 S_diag_j[A_diag_i[i]] = -1;
695 if ((row_sum > max_row_sum) && (max_row_sum < 1.0))
696 {
697 /* make all dependencies weak */
698 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
699 {
700 S_diag_j[jA] = -1;
701 }
702 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
703 {
704 S_offd_j[jA] = -1;
705 }
706 }
707 else
708 {
709 if (num_functions > 1)
710 {
711 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
712 {
713 if (fabs(A_diag_data[jA]) <= strength_threshold * row_scale
714 || dof_func[i] != dof_func[A_diag_j[jA]])
715 {
716 S_diag_j[jA] = -1;
717 }
718 }
719 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
720 {
721 if (fabs(A_offd_data[jA]) <= strength_threshold * row_scale
722 || dof_func[i] != dof_func_offd[A_offd_j[jA]])
723 {
724 S_offd_j[jA] = -1;
725 }
726 }
727 }
728 else
729 {
730 for (jA = A_diag_i[i]+1; jA < A_diag_i[i+1]; jA++)
731 {
732 if (fabs(A_diag_data[jA]) <= strength_threshold * row_scale)
733 {
734 S_diag_j[jA] = -1;
735 }
736 }
737 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
738 {
739 if (fabs(A_offd_data[jA]) <= strength_threshold * row_scale)
740 {
741 S_offd_j[jA] = -1;
742 }
743 }
744 }
745 }
746 }
747
748 /*--------------------------------------------------------------
749 * "Compress" the strength matrix.
750 *
751 * NOTE: S has *NO DIAGONAL ELEMENT* on any row. Caveat Emptor!
752 *
753 * NOTE: This "compression" section of code may be removed, and
754 * coarsening will still be done correctly. However, the routine
755 * that builds interpolation would have to be modified first.
756 *----------------------------------------------------------------*/
757
758/* RDF: not sure if able to thread this loop */
759 jS = 0;
760 for (i = 0; i < num_variables; i++)
761 {
762 S_diag_i[i] = jS;
763 for (jA = A_diag_i[i]; jA < A_diag_i[i+1]; jA++)
764 {
765 if (S_diag_j[jA] > -1)
766 {
767 S_diag_j[jS] = S_diag_j[jA];
768 jS++;
769 }
770 }
771 }
772 S_diag_i[num_variables] = jS;
773 hypre_CSRMatrixNumNonzeros(S_diag) = jS;
774
775/* RDF: not sure if able to thread this loop */
776 jS = 0;
777 for (i = 0; i < num_variables; i++)
778 {
779 S_offd_i[i] = jS;
780 for (jA = A_offd_i[i]; jA < A_offd_i[i+1]; jA++)
781 {
782 if (S_offd_j[jA] > -1)
783 {
784 S_offd_j[jS] = S_offd_j[jA];
785 jS++;
786 }
787 }
788 }
789 S_offd_i[num_variables] = jS;
790 hypre_CSRMatrixNumNonzeros(S_offd) = jS;
791 hypre_ParCSRMatrixCommPkg(S) = NULL;
792
793 *S_ptr = S;
794
795 hypre_TFree(dof_func_offd);
796
797 return (ierr);
798}
799
800/*--------------------------------------------------------------------------*/
801
802int
803hypre_BoomerAMGCreateSCommPkg(hypre_ParCSRMatrix *A,
804 hypre_ParCSRMatrix *S,
805 int **col_offd_S_to_A_ptr)
806{
807 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
808 MPI_Status *status;
809 MPI_Request *requests;
810 hypre_ParCSRCommPkg *comm_pkg_A = hypre_ParCSRMatrixCommPkg(A);
811 hypre_ParCSRCommPkg *comm_pkg_S;
812 hypre_ParCSRCommHandle *comm_handle;
813 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
814 HYPRE_BigInt *col_map_offd_A = hypre_ParCSRMatrixColMapOffd(A);
815
816 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
817 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
818 int *S_offd_i = hypre_CSRMatrixI(S_offd);
819 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
820 HYPRE_BigInt *col_map_offd_S = hypre_ParCSRMatrixColMapOffd(S);
821
822 int *recv_procs_A = hypre_ParCSRCommPkgRecvProcs(comm_pkg_A);
823 int *recv_vec_starts_A =
824 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_A);
825 int *send_procs_A =
826 hypre_ParCSRCommPkgSendProcs(comm_pkg_A);
827 int *send_map_starts_A =
828 hypre_ParCSRCommPkgSendMapStarts(comm_pkg_A);
829 int *recv_procs_S;
830 int *recv_vec_starts_S;
831 int *send_procs_S;
832 int *send_map_starts_S;
833 int *send_map_elmts_S = NULL;
834 HYPRE_BigInt *big_send_map_elmts_S = NULL;
835 int *col_offd_S_to_A;
836
837 int *S_marker;
838 int *send_change;
839 int *recv_change;
840
841 int num_variables = hypre_CSRMatrixNumRows(S_diag);
842 int num_cols_offd_A = hypre_CSRMatrixNumCols(A_offd);
843 int num_cols_offd_S;
844 int i, j, jcol;
845 int proc, cnt, proc_cnt, total_nz;
846 HYPRE_BigInt first_row;
847
848 int ierr = 0;
849
850 int num_sends_A = hypre_ParCSRCommPkgNumSends(comm_pkg_A);
851 int num_recvs_A = hypre_ParCSRCommPkgNumRecvs(comm_pkg_A);
852 int num_sends_S;
853 int num_recvs_S;
854 int num_nonzeros;
855
856 num_nonzeros = S_offd_i[num_variables];
857
858 S_marker = NULL;
859 if (num_cols_offd_A)
860 S_marker = hypre_CTAlloc(int,num_cols_offd_A);
861
862 for (i=0; i < num_cols_offd_A; i++)
863 S_marker[i] = -1;
864
865 for (i=0; i < num_nonzeros; i++)
866 {
867 jcol = S_offd_j[i];
868 S_marker[jcol] = 0;
869 }
870
871 proc = 0;
872 proc_cnt = 0;
873 cnt = 0;
874 num_recvs_S = 0;
875 for (i=0; i < num_recvs_A; i++)
876 {
877 for (j=recv_vec_starts_A[i]; j < recv_vec_starts_A[i+1]; j++)
878 {
879 if (!S_marker[j])
880 {
881 S_marker[j] = cnt;
882 cnt++;
883 proc = 1;
884 }
885 }
886 if (proc) {num_recvs_S++; proc = 0;}
887 }
888
889
890 num_cols_offd_S = cnt;
891 recv_change = NULL;
892 recv_procs_S = NULL;
893 send_change = NULL;
894 if (col_map_offd_S) hypre_TFree(col_map_offd_S);
895 col_map_offd_S = NULL;
896 col_offd_S_to_A = NULL;
897 if (num_recvs_A) recv_change = hypre_CTAlloc(int, num_recvs_A);
898 if (num_sends_A) send_change = hypre_CTAlloc(int, num_sends_A);
899 if (num_recvs_S) recv_procs_S = hypre_CTAlloc(int, num_recvs_S);
900 recv_vec_starts_S = hypre_CTAlloc(int, num_recvs_S+1);
901 if (num_cols_offd_S)
902 {
903 col_map_offd_S = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd_S);
904 col_offd_S_to_A = hypre_CTAlloc(int,num_cols_offd_S);
905 }
906 if (num_cols_offd_S < num_cols_offd_A)
907 {
908 for (i=0; i < num_nonzeros; i++)
909 {
910 jcol = S_offd_j[i];
911 S_offd_j[i] = S_marker[jcol];
912 }
913
914 proc = 0;
915 proc_cnt = 0;
916 cnt = 0;
917 recv_vec_starts_S[0] = 0;
918 for (i=0; i < num_recvs_A; i++)
919 {
920 for (j=recv_vec_starts_A[i]; j < recv_vec_starts_A[i+1]; j++)
921 {
922 if (S_marker[j] != -1)
923 {
924 col_map_offd_S[cnt] = col_map_offd_A[j];
925 col_offd_S_to_A[cnt++] = j;
926 proc = 1;
927 }
928 }
929 recv_change[i] = j-cnt-recv_vec_starts_A[i]
930 +recv_vec_starts_S[proc_cnt];
931 if (proc)
932 {
933 recv_procs_S[proc_cnt++] = recv_procs_A[i];
934 recv_vec_starts_S[proc_cnt] = cnt;
935 proc = 0;
936 }
937 }
938 }
939 else
940 {
941 for (i=0; i < num_recvs_A; i++)
942 {
943 for (j=recv_vec_starts_A[i]; j < recv_vec_starts_A[i+1]; j++)
944 {
945 col_map_offd_S[j] = col_map_offd_A[j];
946 col_offd_S_to_A[j] = j;
947 }
948 recv_procs_S[i] = recv_procs_A[i];
949 recv_vec_starts_S[i] = recv_vec_starts_A[i];
950 }
951 recv_vec_starts_S[num_recvs_A] = recv_vec_starts_A[num_recvs_A];
952 }
953
954 requests = hypre_CTAlloc(MPI_Request,num_sends_A+num_recvs_A);
955 j=0;
956 for (i=0; i < num_sends_A; i++)
957 MPI_Irecv(&send_change[i],1,MPI_INT,send_procs_A[i],
958 0,comm,&requests[j++]);
959
960 for (i=0; i < num_recvs_A; i++)
961 MPI_Isend(&recv_change[i],1,MPI_INT,recv_procs_A[i],
962 0,comm,&requests[j++]);
963
964 status = hypre_CTAlloc(MPI_Status,j);
965 MPI_Waitall(j,requests,status);
966 hypre_TFree(status);
967 hypre_TFree(requests);
968
969 num_sends_S = 0;
970 total_nz = send_map_starts_A[num_sends_A];
971 for (i=0; i < num_sends_A; i++)
972 {
973 if (send_change[i])
974 {
975 if ((send_map_starts_A[i+1]-send_map_starts_A[i]) > send_change[i])
976 num_sends_S++;
977 }
978 else
979 num_sends_S++;
980 total_nz -= send_change[i];
981 }
982
983 send_procs_S = NULL;
984 if (num_sends_S)
985 send_procs_S = hypre_CTAlloc(int,num_sends_S);
986 send_map_starts_S = hypre_CTAlloc(int,num_sends_S+1);
987 send_map_elmts_S = NULL;
988 if (total_nz)
989 {
990 send_map_elmts_S = hypre_CTAlloc(int,total_nz);
991 big_send_map_elmts_S = hypre_CTAlloc(HYPRE_BigInt,total_nz);
992 }
993
994
995 proc = 0;
996 proc_cnt = 0;
997 for (i=0; i < num_sends_A; i++)
998 {
999 cnt = send_map_starts_A[i+1]-send_map_starts_A[i]-send_change[i];
1000 if (cnt)
1001 {
1002 send_procs_S[proc_cnt++] = send_procs_A[i];
1003 send_map_starts_S[proc_cnt] = send_map_starts_S[proc_cnt-1]+cnt;
1004 }
1005 }
1006
1007 comm_pkg_S = hypre_CTAlloc(hypre_ParCSRCommPkg,1);
1008 hypre_ParCSRCommPkgComm(comm_pkg_S) = comm;
1009 hypre_ParCSRCommPkgNumRecvs(comm_pkg_S) = num_recvs_S;
1010 hypre_ParCSRCommPkgRecvProcs(comm_pkg_S) = recv_procs_S;
1011 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_S) = recv_vec_starts_S;
1012 hypre_ParCSRCommPkgNumSends(comm_pkg_S) = num_sends_S;
1013 hypre_ParCSRCommPkgSendProcs(comm_pkg_S) = send_procs_S;
1014 hypre_ParCSRCommPkgSendMapStarts(comm_pkg_S) = send_map_starts_S;
1015
1016 comm_handle = hypre_ParCSRCommHandleCreate(22, comm_pkg_S, col_map_offd_S,
1017 big_send_map_elmts_S);
1018 hypre_ParCSRCommHandleDestroy(comm_handle);
1019
1020 first_row = hypre_ParCSRMatrixFirstRowIndex(A);
1021 if (first_row)
1022 for (i=0; i < send_map_starts_S[num_sends_S]; i++)
1023 send_map_elmts_S[i] = (int)(big_send_map_elmts_S[i]-first_row);
1024
1025 hypre_ParCSRCommPkgSendMapElmts(comm_pkg_S) = send_map_elmts_S;
1026
1027 hypre_ParCSRMatrixCommPkg(S) = comm_pkg_S;
1028 hypre_ParCSRMatrixColMapOffd(S) = col_map_offd_S;
1029 hypre_CSRMatrixNumCols(S_offd) = num_cols_offd_S;
1030
1031 hypre_TFree(S_marker);
1032 hypre_TFree(send_change);
1033 hypre_TFree(recv_change);
1034 hypre_TFree(big_send_map_elmts_S);
1035
1036 *col_offd_S_to_A_ptr = col_offd_S_to_A;
1037
1038 return ierr;
1039}
1040
1041/*--------------------------------------------------------------------------
1042 * hypre_BoomerAMGCreate2ndS : creates strength matrix on coarse points
1043 * for second coarsening pass in aggressive coarsening (S*S+2S)
1044 *--------------------------------------------------------------------------*/
1045
1046int hypre_BoomerAMGCreate2ndS( hypre_ParCSRMatrix *S, int *CF_marker,
1047 int num_paths, HYPRE_BigInt *coarse_row_starts, hypre_ParCSRMatrix **C_ptr)
1048{
1049 MPI_Comm comm = hypre_ParCSRMatrixComm(S);
1050 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(S);
1051 hypre_ParCSRCommPkg *tmp_comm_pkg;
1052 hypre_ParCSRCommHandle *comm_handle;
1053
1054 hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S);
1055
1056 int *S_diag_i = hypre_CSRMatrixI(S_diag);
1057 int *S_diag_j = hypre_CSRMatrixJ(S_diag);
1058
1059 hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S);
1060
1061 int *S_offd_i = hypre_CSRMatrixI(S_offd);
1062 int *S_offd_j = hypre_CSRMatrixJ(S_offd);
1063
1064 int num_cols_diag_S = hypre_CSRMatrixNumCols(S_diag);
1065 int num_cols_offd_S = hypre_CSRMatrixNumCols(S_offd);
1066
1067 hypre_ParCSRMatrix *S2;
1068 HYPRE_BigInt *col_map_offd_C = NULL;
1069
1070 hypre_CSRMatrix *C_diag;
1071
1072 int *C_diag_data = NULL;
1073 int *C_diag_i;
1074 int *C_diag_j = NULL;
1075
1076 hypre_CSRMatrix *C_offd;
1077
1078 int *C_offd_data=NULL;
1079 int *C_offd_i;
1080 int *C_offd_j=NULL;
1081
1082 int C_diag_size;
1083 int C_offd_size;
1084 int num_cols_offd_C = 0;
1085
1086 int *S_ext_diag_i = NULL;
1087 int *S_ext_diag_j = NULL;
1088 int S_ext_diag_size = 0;
1089
1090 int *S_ext_offd_i = NULL;
1091 int *S_ext_offd_j = NULL;
1092 int S_ext_offd_size = 0;
1093
1094 int *CF_marker_offd = NULL;
1095
1096 int *S_marker = NULL;
1097 int *S_marker_offd = NULL;
1098
1099 int *fine_to_coarse = NULL;
1100 HYPRE_BigInt *big_fine_to_coarse_offd = NULL;
1101 int *map_S_to_C = NULL;
1102
1103 int num_sends = 0;
1104 int num_recvs = 0;
1105 int *send_map_starts;
1106 int *tmp_send_map_starts = NULL;
1107 int *send_map_elmts;
1108 int *recv_vec_starts;
1109 int *tmp_recv_vec_starts = NULL;
1110 int *int_buf_data = NULL;
1111 HYPRE_BigInt *big_int_buf_data = NULL;
1112 HYPRE_BigInt *temp = NULL;
1113
1114 int i, j, k;
1115 int i1, i2, i3;
1116 HYPRE_BigInt big_i1;
1117 int jj1, jj2, jcol, jrow, j_cnt;
1118
1119 int jj_count_diag, jj_count_offd;
1120 int jj_row_begin_diag, jj_row_begin_offd;
1121 int cnt, cnt_offd, cnt_diag;
1122 int num_procs, my_id;
1123 int index;
1124 HYPRE_BigInt value;
1125 int num_coarse;
1126 int num_coarse_offd;
1127 int num_nonzeros;
1128 int num_nonzeros_diag;
1129 int num_nonzeros_offd;
1130 int num_nnz_diag;
1131 int num_nnz_offd;
1132 int num_threads;
1133 int ns, ne, rest, size;
1134 int *ns_thr, *ne_thr;
1135 int *nnz, *nnz_offd, *cnt_coarse;
1136 HYPRE_BigInt global_num_coarse;
1137 HYPRE_BigInt my_first_cpt, my_last_cpt;
1138
1139 int *S_int_i = NULL;
1140 HYPRE_BigInt *S_int_j = NULL;
1141 int *S_ext_i = NULL;
1142 HYPRE_BigInt *S_ext_j = NULL;
1143
1144 /*-----------------------------------------------------------------------
1145 * Extract S_ext, i.e. portion of B that is stored on neighbor procs
1146 * and needed locally for matrix matrix product
1147 *-----------------------------------------------------------------------*/
1148
1149 MPI_Comm_size(comm, &num_procs);
1150 MPI_Comm_rank(comm, &my_id);
1151 num_threads = hypre_NumThreads();
1152
1153#ifdef HYPRE_NO_GLOBAL_PARTITION
1154 my_first_cpt = coarse_row_starts[0];
1155 my_last_cpt = coarse_row_starts[1]-1;
1156 if (my_id == (num_procs -1)) global_num_coarse = coarse_row_starts[1];
1157 MPI_Bcast(&global_num_coarse, 1, MPI_HYPRE_BIG_INT, num_procs-1, comm);
1158#else
1159 my_first_cpt = coarse_row_starts[my_id];
1160 my_last_cpt = coarse_row_starts[my_id+1]-1;
1161 global_num_coarse = coarse_row_starts[num_procs];
1162#endif
1163
1164 if (num_cols_offd_S)
1165 {
1166 CF_marker_offd = hypre_CTAlloc(int, num_cols_offd_S);
1167 big_fine_to_coarse_offd = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd_S);
1168 }
1169
1170 if (num_cols_diag_S) fine_to_coarse = hypre_CTAlloc(int, num_cols_diag_S);
1171
1172 num_coarse = 0;
1173 for (i=0; i < num_cols_diag_S; i++)
1174 {
1175 if (CF_marker[i] > 0)
1176 {
1177 fine_to_coarse[i] = num_coarse;
1178 num_coarse++;
1179 }
1180 else
1181 {
1182 fine_to_coarse[i] = -1;
1183 }
1184 }
1185
1186 if (num_procs > 1)
1187 {
1188 if (!comm_pkg)
1189 {
1190#ifdef HYPRE_NO_GLOBAL_PARTITION
1191 hypre_NewCommPkgCreate(S);
1192#else
1193 hypre_MatvecCommPkgCreate(S);
1194#endif
1195 comm_pkg = hypre_ParCSRMatrixCommPkg(S);
1196 }
1197 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1198 send_map_starts = hypre_ParCSRCommPkgSendMapStarts(comm_pkg);
1199 send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
1200 num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
1201 recv_vec_starts = hypre_ParCSRCommPkgRecvVecStarts(comm_pkg);
1202 big_int_buf_data = hypre_CTAlloc(HYPRE_BigInt, send_map_starts[num_sends]);
1203
1204 index = 0;
1205 for (i = 0; i < num_sends; i++)
1206 {
1207 for (j = send_map_starts[i]; j < send_map_starts[i+1]; j++)
1208 big_int_buf_data[index++]
1209 = my_first_cpt+ (HYPRE_BigInt)fine_to_coarse[send_map_elmts[j]];
1210 }
1211
1212 comm_handle = hypre_ParCSRCommHandleCreate( 21, comm_pkg, big_int_buf_data,
1213 big_fine_to_coarse_offd);
1214
1215 hypre_ParCSRCommHandleDestroy(comm_handle);
1216 hypre_TFree(big_int_buf_data);
1217 int_buf_data = hypre_CTAlloc(int, send_map_starts[num_sends]);
1218
1219 index = 0;
1220 for (i = 0; i < num_sends; i++)
1221 {
1222 for (j = send_map_starts[i]; j < send_map_starts[i+1]; j++)
1223 {
1224 int_buf_data[index++] = CF_marker[send_map_elmts[j]];
1225 }
1226 }
1227
1228 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
1229 CF_marker_offd);
1230
1231 hypre_ParCSRCommHandleDestroy(comm_handle);
1232 hypre_TFree(int_buf_data);
1233
1234 S_int_i = hypre_CTAlloc(int, send_map_starts[num_sends]+1);
1235 S_ext_i = hypre_CTAlloc(int, recv_vec_starts[num_recvs]+1);
1236
1237/*--------------------------------------------------------------------------
1238 * generate S_int_i through adding number of coarse row-elements of offd and diag
1239 * for corresponding rows. S_int_i[j+1] contains the number of coarse elements of
1240 * a row j (which is determined through send_map_elmts)
1241 *--------------------------------------------------------------------------*/
1242 S_int_i[0] = 0;
1243 j_cnt = 0;
1244 num_nonzeros = 0;
1245 for (i=0; i < num_sends; i++)
1246 {
1247 for (j = send_map_starts[i]; j < send_map_starts[i+1]; j++)
1248 {
1249 jrow = send_map_elmts[j];
1250 index = 0;
1251 for (k = S_diag_i[jrow]; k < S_diag_i[jrow+1]; k++)
1252 {
1253 if (CF_marker[S_diag_j[k]] > 0) index++;
1254 }
1255 for (k = S_offd_i[jrow]; k < S_offd_i[jrow+1]; k++)
1256 {
1257 if (CF_marker_offd[S_offd_j[k]] > 0) index++;
1258 }
1259 S_int_i[++j_cnt] = index;
1260 num_nonzeros += S_int_i[j_cnt];
1261 }
1262 }
1263
1264/*--------------------------------------------------------------------------
1265 * initialize communication
1266 *--------------------------------------------------------------------------*/
1267 if (num_procs > 1)
1268 comm_handle =
1269 hypre_ParCSRCommHandleCreate(11,comm_pkg,&S_int_i[1],&S_ext_i[1]);
1270
1271 if (num_nonzeros) S_int_j = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros);
1272
1273 tmp_send_map_starts = hypre_CTAlloc(int, num_sends+1);
1274 tmp_recv_vec_starts = hypre_CTAlloc(int, num_recvs+1);
1275
1276 tmp_send_map_starts[0] = 0;
1277 j_cnt = 0;
1278 for (i=0; i < num_sends; i++)
1279 {
1280 for (j = send_map_starts[i]; j < send_map_starts[i+1]; j++)
1281 {
1282 jrow = send_map_elmts[j];
1283 for (k=S_diag_i[jrow]; k < S_diag_i[jrow+1]; k++)
1284 {
1285 if (CF_marker[S_diag_j[k]] > 0)
1286 S_int_j[j_cnt++] = (HYPRE_BigInt)fine_to_coarse[S_diag_j[k]]+my_first_cpt;
1287 }
1288 for (k=S_offd_i[jrow]; k < S_offd_i[jrow+1]; k++)
1289 {
1290 if (CF_marker_offd[S_offd_j[k]] > 0)
1291 S_int_j[j_cnt++] = big_fine_to_coarse_offd[S_offd_j[k]];
1292 }
1293 }
1294 tmp_send_map_starts[i+1] = j_cnt;
1295 }
1296
1297 tmp_comm_pkg = hypre_CTAlloc(hypre_ParCSRCommPkg,1);
1298 hypre_ParCSRCommPkgComm(tmp_comm_pkg) = comm;
1299 hypre_ParCSRCommPkgNumSends(tmp_comm_pkg) = num_sends;
1300 hypre_ParCSRCommPkgNumRecvs(tmp_comm_pkg) = num_recvs;
1301 hypre_ParCSRCommPkgSendProcs(tmp_comm_pkg) =
1302 hypre_ParCSRCommPkgSendProcs(comm_pkg);
1303 hypre_ParCSRCommPkgRecvProcs(tmp_comm_pkg) =
1304 hypre_ParCSRCommPkgRecvProcs(comm_pkg);
1305 hypre_ParCSRCommPkgSendMapStarts(tmp_comm_pkg) = tmp_send_map_starts;
1306
1307 hypre_ParCSRCommHandleDestroy(comm_handle);
1308 comm_handle = NULL;
1309/*--------------------------------------------------------------------------
1310 * after communication exchange S_ext_i[j+1] contains the number of coarse elements
1311 * of a row j !
1312 * evaluate S_ext_i and compute num_nonzeros for S_ext
1313 *--------------------------------------------------------------------------*/
1314
1315 for (i=0; i < recv_vec_starts[num_recvs]; i++)
1316 S_ext_i[i+1] += S_ext_i[i];
1317
1318 num_nonzeros = S_ext_i[recv_vec_starts[num_recvs]];
1319
1320 if (num_nonzeros) S_ext_j = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros);
1321
1322 tmp_recv_vec_starts[0] = 0;
1323 for (i=0; i < num_recvs; i++)
1324 tmp_recv_vec_starts[i+1] = S_ext_i[recv_vec_starts[i+1]];
1325
1326 hypre_ParCSRCommPkgRecvVecStarts(tmp_comm_pkg) = tmp_recv_vec_starts;
1327
1328 comm_handle = hypre_ParCSRCommHandleCreate(21,tmp_comm_pkg,S_int_j,S_ext_j);
1329 hypre_ParCSRCommHandleDestroy(comm_handle);
1330 comm_handle = NULL;
1331
1332 hypre_TFree(tmp_send_map_starts);
1333 hypre_TFree(tmp_recv_vec_starts);
1334 hypre_TFree(tmp_comm_pkg);
1335
1336 hypre_TFree(S_int_i);
1337 hypre_TFree(S_int_j);
1338
1339 S_ext_diag_size = 0;
1340 S_ext_offd_size = 0;
1341
1342 for (i=0; i < num_cols_offd_S; i++)
1343 {
1344 for (j=S_ext_i[i]; j < S_ext_i[i+1]; j++)
1345 {
1346 if (S_ext_j[j] < my_first_cpt || S_ext_j[j] > my_last_cpt)
1347 S_ext_offd_size++;
1348 else
1349 S_ext_diag_size++;
1350 }
1351 }
1352 S_ext_diag_i = hypre_CTAlloc(int, num_cols_offd_S+1);
1353 S_ext_offd_i = hypre_CTAlloc(int, num_cols_offd_S+1);
1354
1355 if (S_ext_diag_size)
1356 {
1357 S_ext_diag_j = hypre_CTAlloc(int, S_ext_diag_size);
1358 }
1359 if (S_ext_offd_size)
1360 {
1361 S_ext_offd_j = hypre_CTAlloc(int, S_ext_offd_size);
1362 }
1363
1364 cnt_offd = 0;
1365 cnt_diag = 0;
1366 cnt = 0;
1367 num_coarse_offd = 0;
1368 for (i=0; i < num_cols_offd_S; i++)
1369 {
1370 if (CF_marker_offd[i] > 0) num_coarse_offd++;
1371
1372 for (j=S_ext_i[i]; j < S_ext_i[i+1]; j++)
1373 {
1374 big_i1 = S_ext_j[j];
1375 if (big_i1 < my_first_cpt || big_i1 > my_last_cpt)
1376 S_ext_j[cnt_offd++] = big_i1;
1377 else
1378 S_ext_diag_j[cnt_diag++] = (int)(big_i1 - my_first_cpt);
1379 }
1380 S_ext_diag_i[++cnt] = cnt_diag;
1381 S_ext_offd_i[cnt] = cnt_offd;
1382 }
1383
1384 hypre_TFree(S_ext_i);
1385
1386 if (S_ext_offd_size || num_coarse_offd)
1387 temp = hypre_CTAlloc(HYPRE_BigInt, S_ext_offd_size+num_coarse_offd);
1388 for (i=0; i < S_ext_offd_size; i++)
1389 temp[i] = S_ext_j[i];
1390 cnt = S_ext_offd_size;
1391 for (i=0; i < num_cols_offd_S; i++)
1392 if (CF_marker_offd[i] > 0) temp[cnt++] = big_fine_to_coarse_offd[i];
1393 if (cnt)
1394 {
1395 hypre_BigQsort0(temp, 0, cnt-1);
1396
1397 num_cols_offd_C = 1;
1398 value = temp[0];
1399 for (i=1; i < cnt; i++)
1400 {
1401 if (temp[i] > value)
1402 {
1403 value = temp[i];
1404 temp[num_cols_offd_C++] = value;
1405 }
1406 }
1407 }
1408
1409 if (num_cols_offd_C)
1410 col_map_offd_C = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd_C);
1411
1412 for (i=0; i < num_cols_offd_C; i++)
1413 col_map_offd_C[i] = temp[i];
1414
1415 hypre_TFree(temp);
1416
1417#define HYPRE_SMP_PRIVATE i
1418#include "../utilities/hypre_smp_forloop.h"
1419 for (i=0 ; i < S_ext_offd_size; i++)
1420 S_ext_offd_j[i] = hypre_BigBinarySearch(col_map_offd_C,
1421 S_ext_j[i],
1422 num_cols_offd_C);
1423 hypre_TFree(S_ext_j);
1424 if (num_cols_offd_S)
1425 {
1426 map_S_to_C = hypre_CTAlloc(int,num_cols_offd_S);
1427
1428 cnt = 0;
1429 for (i=0; i < num_cols_offd_S; i++)
1430 {
1431 if (CF_marker_offd[i] > 0)
1432 {
1433 while (big_fine_to_coarse_offd[i] > col_map_offd_C[cnt])
1434 {
1435 cnt++;
1436 }
1437 map_S_to_C[i] = cnt++;
1438 }
1439 else
1440 {
1441 map_S_to_C[i] = -1;
1442 }
1443 }
1444 }
1445 }
1446
1447 /*-----------------------------------------------------------------------
1448 * Allocate and initialize some stuff.
1449 *-----------------------------------------------------------------------*/
1450
1451 /*if (num_coarse) S_marker = hypre_CTAlloc(int, num_coarse);
1452
1453 for (i1 = 0; i1 < num_coarse; i1++)
1454 S_marker[i1] = -1;
1455
1456 if (num_cols_offd_C) S_marker_offd = hypre_CTAlloc(int, num_cols_offd_C);
1457
1458 for (i1 = 0; i1 < num_cols_offd_C; i1++)
1459 S_marker_offd[i1] = -1;*/
1460
1461 C_diag_i = hypre_CTAlloc(int, num_coarse+1);
1462 C_offd_i = hypre_CTAlloc(int, num_coarse+1);
1463
1464
1465 /*-----------------------------------------------------------------------
1466 * Loop over rows of S
1467 *-----------------------------------------------------------------------*/
1468
1469 cnt = 0;
1470 num_nonzeros_diag = 0;
1471 num_nonzeros_offd = 0;
1472
1473 ns_thr = hypre_CTAlloc(int, num_threads);
1474 ne_thr = hypre_CTAlloc(int, num_threads);
1475 nnz = hypre_CTAlloc(int, num_threads);
1476 nnz_offd = hypre_CTAlloc(int, num_threads);
1477 cnt_coarse = hypre_CTAlloc(int, num_threads);
1478
1479 size = num_cols_diag_S/num_threads;
1480 rest = num_cols_diag_S - size*num_threads;
1481
1482 for (k = 0; k < num_threads; k++)
1483 {
1484 if (k < rest)
1485 {
1486 ns_thr[k] = k*size+k;
1487 ne_thr[k] = (k+1)*size+k+1;
1488 }
1489 else
1490 {
1491 ns_thr[k] = k*size+rest;
1492 ne_thr[k] = (k+1)*size+rest;
1493 }
1494 }
1495
1496#define HYPRE_SMP_PRIVATE k,i1,i2,i3,jj1,jj2,jcol,ns,ne,S_marker,S_marker_offd,num_nnz_diag,num_nnz_offd,cnt
1497#include "../utilities/hypre_smp_forloop.h"
1498 for (k = 0; k < num_threads; k++)
1499 {
1500 ns = ns_thr[k];
1501 ne = ne_thr[k];
1502 num_nnz_diag = 0;
1503 num_nnz_offd = 0;
1504 cnt_coarse[k] = 0;
1505 /*for (i1 = 0; i1 < num_cols_diag_S; i1++)
1506 {*/
1507 S_marker = NULL;
1508 if (num_coarse) S_marker = hypre_CTAlloc(int, num_coarse);
1509
1510 for (i1 = 0; i1 < num_coarse; i1++)
1511 S_marker[i1] = -1;
1512
1513 S_marker_offd = NULL;
1514 if (num_cols_offd_C) S_marker_offd = hypre_CTAlloc(int, num_cols_offd_C);
1515
1516 for (i1 = 0; i1 < num_cols_offd_C; i1++)
1517 S_marker_offd[i1] = -1;
1518
1519 for (i1 = ns; i1 < ne; i1++)
1520 {
1521 if (CF_marker[i1] > 0)
1522 {
1523 cnt_coarse[k]++;
1524 for (jj1 = S_diag_i[i1]; jj1 < S_diag_i[i1+1]; jj1++)
1525 {
1526 jcol = S_diag_j[jj1];
1527 if (CF_marker[jcol] > 0)
1528 {
1529 S_marker[fine_to_coarse[jcol]] = i1;
1530 num_nnz_diag++;
1531 }
1532 }
1533 for (jj1 = S_offd_i[i1]; jj1 < S_offd_i[i1+1]; jj1++)
1534 {
1535 jcol = S_offd_j[jj1];
1536 if (CF_marker_offd[jcol] > 0)
1537 {
1538 S_marker_offd[map_S_to_C[jcol]] = i1;
1539 num_nnz_offd++;
1540 }
1541 }
1542 for (jj1 = S_diag_i[i1]; jj1 < S_diag_i[i1+1]; jj1++)
1543 {
1544 i2 = S_diag_j[jj1];
1545 for (jj2 = S_diag_i[i2]; jj2 < S_diag_i[i2+1]; jj2++)
1546 {
1547 i3 = S_diag_j[jj2];
1548 if (CF_marker[i3] > 0 && S_marker[fine_to_coarse[i3]] != i1)
1549 {
1550 S_marker[fine_to_coarse[i3]] = i1;
1551 num_nnz_diag++;
1552 }
1553 }
1554 for (jj2 = S_offd_i[i2]; jj2 < S_offd_i[i2+1]; jj2++)
1555 {
1556 i3 = S_offd_j[jj2];
1557 if (CF_marker_offd[i3] > 0 &&
1558 S_marker_offd[map_S_to_C[i3]] != i1)
1559 {
1560 S_marker_offd[map_S_to_C[i3]] = i1;
1561 num_nnz_offd++;
1562 }
1563 }
1564 }
1565 for (jj1 = S_offd_i[i1]; jj1 < S_offd_i[i1+1]; jj1++)
1566 {
1567 i2 = S_offd_j[jj1];
1568 for (jj2 = S_ext_diag_i[i2]; jj2 < S_ext_diag_i[i2+1]; jj2++)
1569 {
1570 i3 = S_ext_diag_j[jj2];
1571 if (S_marker[i3] != i1)
1572 {
1573 S_marker[i3] = i1;
1574 num_nnz_diag++;
1575 }
1576 }
1577 for (jj2 = S_ext_offd_i[i2]; jj2 < S_ext_offd_i[i2+1]; jj2++)
1578 {
1579 i3 = S_ext_offd_j[jj2];
1580 if (S_marker_offd[i3] != i1)
1581 {
1582 S_marker_offd[i3] = i1;
1583 num_nnz_offd++;
1584 }
1585 }
1586 }
1587 }
1588 }
1589 nnz[k] = num_nnz_diag;
1590 nnz_offd[k] = num_nnz_offd;
1591 if (num_coarse) hypre_TFree(S_marker);
1592 if (num_cols_offd_C) hypre_TFree(S_marker_offd);
1593 }
1594
1595 num_nonzeros_diag = 0;
1596 num_nonzeros_offd = 0;
1597 for (k = 0; k < num_threads; k++)
1598 {
1599 num_nonzeros_diag += nnz[k];
1600 num_nonzeros_offd += nnz_offd[k];
1601 }
1602 C_diag_i[num_coarse] = num_nonzeros_diag;
1603 C_offd_i[num_coarse] = num_nonzeros_offd;
1604
1605 if (num_nonzeros_diag)
1606 {
1607 C_diag_j = hypre_CTAlloc(int,num_nonzeros_diag);
1608 C_diag_data = hypre_CTAlloc(int,num_nonzeros_diag);
1609 }
1610 if (num_nonzeros_offd)
1611 {
1612 C_offd_j = hypre_CTAlloc(int,num_nonzeros_offd);
1613 C_offd_data = hypre_CTAlloc(int,num_nonzeros_offd);
1614 }
1615
1616#define HYPRE_SMP_PRIVATE k,i1,i2,i3,jj1,jj2,jcol,ns,ne,S_marker,S_marker_offd,jj_count_diag,jj_count_offd,jj_row_begin_diag,jj_row_begin_offd,index,cnt
1617#include "../utilities/hypre_smp_forloop.h"
1618 for (k = 0; k < num_threads; k++)
1619 {
1620 ns = ns_thr[k];
1621 ne = ne_thr[k];
1622 jj_count_diag = 0;
1623 jj_count_offd = 0;
1624 for (i1 = 0; i1 < k; i1++)
1625 {
1626 jj_count_diag += nnz[i1];
1627 jj_count_offd += nnz_offd[i1];
1628 }
1629
1630 S_marker = NULL;
1631 if (num_coarse) S_marker = hypre_CTAlloc(int, num_coarse);
1632
1633 for (i1 = 0; i1 < num_coarse; i1++)
1634 S_marker[i1] = -1;
1635
1636 S_marker_offd = NULL;
1637 if (num_cols_offd_C) S_marker_offd = hypre_CTAlloc(int, num_cols_offd_C);
1638
1639 for (i1 = 0; i1 < num_cols_offd_C; i1++)
1640 S_marker_offd[i1] = -1;
1641
1642 cnt = 0;
1643 for (i1 = 0; i1 < k; i1++)
1644 cnt += cnt_coarse[i1];
1645 for (i1 = ns; i1 < ne; i1++)
1646 {
1647 /*for (i1 = 0; i1 < num_cols_diag_S; i1++)
1648 {*/
1649
1650 /*--------------------------------------------------------------------
1651 * Set marker for diagonal entry, C_{i1,i1} (for square matrices).
1652 *--------------------------------------------------------------------*/
1653
1654 jj_row_begin_diag = jj_count_diag;
1655 jj_row_begin_offd = jj_count_offd;
1656
1657 if (CF_marker[i1] > 0)
1658 {
1659 C_diag_i[cnt] = jj_row_begin_diag;
1660 C_offd_i[cnt++] = jj_row_begin_offd;
1661 for (jj1 = S_diag_i[i1]; jj1 < S_diag_i[i1+1]; jj1++)
1662 {
1663 jcol = S_diag_j[jj1];
1664 if (CF_marker[jcol] > 0)
1665 {
1666 S_marker[fine_to_coarse[jcol]] = jj_count_diag;
1667 C_diag_j[jj_count_diag] = fine_to_coarse[jcol];
1668 C_diag_data[jj_count_diag] = 2;
1669 jj_count_diag++;
1670 }
1671 }
1672 for (jj1 = S_offd_i[i1]; jj1 < S_offd_i[i1+1]; jj1++)
1673 {
1674 jcol = S_offd_j[jj1];
1675 if (CF_marker_offd[jcol] > 0)
1676 {
1677 index = map_S_to_C[jcol];
1678 S_marker_offd[index] = jj_count_offd;
1679 C_offd_j[jj_count_offd] = index;
1680 C_offd_data[jj_count_offd] = 2;
1681 jj_count_offd++;
1682 }
1683 }
1684 for (jj1 = S_diag_i[i1]; jj1 < S_diag_i[i1+1]; jj1++)
1685 {
1686 i2 = S_diag_j[jj1];
1687 for (jj2 = S_diag_i[i2]; jj2 < S_diag_i[i2+1]; jj2++)
1688 {
1689 i3 = S_diag_j[jj2];
1690 if (CF_marker[i3] > 0)
1691 {
1692 if (S_marker[fine_to_coarse[i3]] < jj_row_begin_diag)
1693 {
1694 S_marker[fine_to_coarse[i3]] = jj_count_diag;
1695 C_diag_j[jj_count_diag] = fine_to_coarse[i3];
1696 C_diag_data[jj_count_diag]++;
1697 jj_count_diag++;
1698 }
1699 else
1700 {
1701 C_diag_data[S_marker[fine_to_coarse[i3]]]++;
1702 }
1703 }
1704 }
1705 for (jj2 = S_offd_i[i2]; jj2 < S_offd_i[i2+1]; jj2++)
1706 {
1707 i3 = S_offd_j[jj2];
1708 if (CF_marker_offd[i3] > 0)
1709 {
1710 index = map_S_to_C[i3];
1711 if (S_marker_offd[index] < jj_row_begin_offd)
1712 {
1713 S_marker_offd[index] = jj_count_offd;
1714 C_offd_j[jj_count_offd] = index;
1715 C_offd_data[jj_count_offd]++;
1716 jj_count_offd++;
1717 }
1718 else
1719 {
1720 C_offd_data[S_marker_offd[index]]++;
1721 }
1722 }
1723 }
1724 }
1725 for (jj1 = S_offd_i[i1]; jj1 < S_offd_i[i1+1]; jj1++)
1726 {
1727 i2 = S_offd_j[jj1];
1728 for (jj2 = S_ext_diag_i[i2]; jj2 < S_ext_diag_i[i2+1]; jj2++)
1729 {
1730 i3 = S_ext_diag_j[jj2];
1731 if (S_marker[i3] < jj_row_begin_diag)
1732 {
1733 S_marker[i3] = jj_count_diag;
1734 C_diag_j[jj_count_diag] = i3;
1735 C_diag_data[jj_count_diag]++;
1736 jj_count_diag++;
1737 }
1738 else
1739 {
1740 C_diag_data[S_marker[i3]]++;
1741 }
1742 }
1743 for (jj2 = S_ext_offd_i[i2]; jj2 < S_ext_offd_i[i2+1]; jj2++)
1744 {
1745 i3 = S_ext_offd_j[jj2];
1746 if (S_marker_offd[i3] < jj_row_begin_offd)
1747 {
1748 S_marker_offd[i3] = jj_count_offd;
1749 C_offd_j[jj_count_offd] = i3;
1750 C_offd_data[jj_count_offd]++;
1751 jj_count_offd++;
1752 }
1753 else
1754 {
1755 C_offd_data[S_marker_offd[i3]]++;
1756 }
1757 }
1758 }
1759 }
1760 }
1761 if (num_coarse) hypre_TFree(S_marker);
1762 if (num_cols_offd_C) hypre_TFree(S_marker_offd);
1763 }
1764
1765
1766 cnt = 0;
1767
1768 for (i=0; i < num_coarse; i++)
1769 {
1770 for (j=C_diag_i[i]; j < C_diag_i[i+1]; j++)
1771 {
1772 jcol = C_diag_j[j];
1773 if (C_diag_data[j] >= num_paths && jcol != i)
1774 C_diag_j[cnt++] = jcol;
1775 }
1776 C_diag_i[i] = cnt;
1777 }
1778
1779 if (num_nonzeros_diag) hypre_TFree(C_diag_data);
1780
1781 for (i=num_coarse; i > 0; i--)
1782 C_diag_i[i] = C_diag_i[i-1];
1783
1784 C_diag_i[0] = 0;
1785
1786 cnt = 0;
1787 for (i=0; i < num_coarse; i++)
1788 {
1789 for (j=C_offd_i[i]; j < C_offd_i[i+1]; j++)
1790 {
1791 jcol = C_offd_j[j];
1792 if (C_offd_data[j] >= num_paths)
1793 C_offd_j[cnt++] = jcol;
1794 }
1795 C_offd_i[i] = cnt;
1796 }
1797 if (num_nonzeros_offd) hypre_TFree(C_offd_data);
1798 for (i=num_coarse; i > 0; i--)
1799 C_offd_i[i] = C_offd_i[i-1];
1800
1801 C_offd_i[0] = 0;
1802
1803 cnt = 0;
1804 for (i=0; i < num_cols_diag_S; i++)
1805 {
1806 if (CF_marker[i] > 0)
1807 {
1808 if (!(C_diag_i[cnt+1]-C_diag_i[cnt]) &&
1809 !(C_offd_i[cnt+1]-C_offd_i[cnt]))
1810 CF_marker[i] = 2;
1811 cnt++;
1812 }
1813 }
1814
1815 C_diag_size = C_diag_i[num_coarse];
1816 C_offd_size = C_offd_i[num_coarse];
1817
1818 S2 = hypre_ParCSRMatrixCreate(comm, global_num_coarse,
1819 global_num_coarse, coarse_row_starts,
1820 coarse_row_starts, num_cols_offd_C, C_diag_size, C_offd_size);
1821
1822 hypre_ParCSRMatrixOwnsRowStarts(S2) = 0;
1823
1824 C_diag = hypre_ParCSRMatrixDiag(S2);
1825 hypre_CSRMatrixI(C_diag) = C_diag_i;
1826 if (num_nonzeros_diag) hypre_CSRMatrixJ(C_diag) = C_diag_j;
1827
1828 C_offd = hypre_ParCSRMatrixOffd(S2);
1829 hypre_CSRMatrixI(C_offd) = C_offd_i;
1830 hypre_ParCSRMatrixOffd(S2) = C_offd;
1831
1832 if (num_cols_offd_C)
1833 {
1834 if (num_nonzeros_offd) hypre_CSRMatrixJ(C_offd) = C_offd_j;
1835 hypre_ParCSRMatrixColMapOffd(S2) = col_map_offd_C;
1836
1837 }
1838
1839 /*-----------------------------------------------------------------------
1840 * Free various arrays
1841 *-----------------------------------------------------------------------*/
1842
1843 hypre_TFree(ns_thr);
1844 hypre_TFree(ne_thr);
1845 hypre_TFree(nnz);
1846 hypre_TFree(nnz_offd);
1847 hypre_TFree(cnt_coarse);
1848 hypre_TFree(S_ext_diag_i);
1849 hypre_TFree(fine_to_coarse);
1850 if (S_ext_diag_size)
1851 {
1852 hypre_TFree(S_ext_diag_j);
1853 }
1854 hypre_TFree(S_ext_offd_i);
1855 if (S_ext_offd_size)
1856 {
1857 hypre_TFree(S_ext_offd_j);
1858 }
1859 if (num_cols_offd_S)
1860 {
1861 hypre_TFree(map_S_to_C);
1862 hypre_TFree(CF_marker_offd);
1863 hypre_TFree(big_fine_to_coarse_offd);
1864 }
1865
1866 *C_ptr = S2;
1867
1868 return 0;
1869
1870}
1871
1872/*--------------------------------------------------------------------------
1873 * hypre_BoomerAMGCorrectCFMarker : corrects CF_marker after aggr. coarsening
1874 *--------------------------------------------------------------------------*/
1875int
1876hypre_BoomerAMGCorrectCFMarker(int *CF_marker, int num_var, int *new_CF_marker)
1877{
1878 int i, cnt;
1879
1880 cnt = 0;
1881 for (i=0; i < num_var; i++)
1882 {
1883 if (CF_marker[i] > 0 )
1884 {
1885 if (CF_marker[i] == 1) CF_marker[i] = new_CF_marker[cnt++];
1886 else { CF_marker[i] = 1; cnt++;}
1887 }
1888 }
1889
1890 return 0;
1891}
1892
1893/*--------------------------------------------------------------------------
1894 * hypre_BoomerAMGCorrectCFMarker2 : corrects CF_marker after aggr. coarsening,
1895 * but marks new F-points (previous C-points) as -2
1896 *--------------------------------------------------------------------------*/
1897int
1898hypre_BoomerAMGCorrectCFMarker2(int *CF_marker, int num_var, int *new_CF_marker){
1899 int i, cnt;
1900
1901 cnt = 0;
1902 for (i=0; i < num_var; i++)
1903 {
1904 if (CF_marker[i] > 0 )
1905 {
1906 if (new_CF_marker[cnt] == -1) CF_marker[i] = -2;
1907 else CF_marker[i] = 1;
1908 cnt++;
1909 }
1910 }
1911
1912 return 0;
1913}
1914
Note: See TracBrowser for help on using the repository browser.