MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
hypre_parcsr.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "../config/config.hpp"
13#include "../general/error.hpp"
14#include "../general/forall.hpp"
15
16#ifdef MFEM_USE_MPI
17
18#include "hypre_parcsr.hpp"
19#include "hypre.hpp"
20#include <limits>
21#include <cmath>
22
23namespace mfem
24{
25namespace internal
26{
27
28/*--------------------------------------------------------------------------
29 * A*X = B style elimination
30 *--------------------------------------------------------------------------*/
31
32/*
33 Function: hypre_CSRMatrixEliminateAXB
34
35 Eliminate the rows and columns of A corresponding to the
36 given sorted (!) list of rows. Put I on the eliminated diagonal,
37 subtract columns times X from B.
38*/
39void hypre_CSRMatrixEliminateAXB(hypre_CSRMatrix *A,
40 HYPRE_Int nrows_to_eliminate,
41 HYPRE_Int *rows_to_eliminate,
42 hypre_Vector *X,
43 hypre_Vector *B)
44{
45 HYPRE_Int i, j;
46 HYPRE_Int irow, jcol, ibeg, iend, pos;
47 HYPRE_Real a;
48
49 HYPRE_Int *Ai = hypre_CSRMatrixI(A);
50 HYPRE_Int *Aj = hypre_CSRMatrixJ(A);
51 HYPRE_Real *Adata = hypre_CSRMatrixData(A);
52 HYPRE_Int nrows = hypre_CSRMatrixNumRows(A);
53
54 HYPRE_Real *Xdata = hypre_VectorData(X);
55 HYPRE_Real *Bdata = hypre_VectorData(B);
56
57 /* eliminate the columns */
58 for (i = 0; i < nrows; i++)
59 {
60 ibeg = Ai[i];
61 iend = Ai[i+1];
62 for (j = ibeg; j < iend; j++)
63 {
64 jcol = Aj[j];
65 pos = hypre_BinarySearch(rows_to_eliminate, jcol, nrows_to_eliminate);
66 if (pos != -1)
67 {
68 a = Adata[j];
69 Adata[j] = 0.0;
70 Bdata[i] -= a * Xdata[jcol];
71 }
72 }
73 }
74
75 /* remove the rows and set the diagonal equal to 1 */
76 for (i = 0; i < nrows_to_eliminate; i++)
77 {
78 irow = rows_to_eliminate[i];
79 ibeg = Ai[irow];
80 iend = Ai[irow+1];
81 for (j = ibeg; j < iend; j++)
82 {
83 if (Aj[j] == irow)
84 {
85 Adata[j] = 1.0;
86 }
87 else
88 {
89 Adata[j] = 0.0;
90 }
91 }
92 }
93}
94
95/*
96 Function: hypre_CSRMatrixEliminateOffdColsAXB
97
98 Eliminate the given sorted (!) list of columns of A, subtract them from B.
99*/
100void hypre_CSRMatrixEliminateOffdColsAXB(hypre_CSRMatrix *A,
101 HYPRE_Int ncols_to_eliminate,
102 HYPRE_Int *eliminate_cols,
103 HYPRE_Real *eliminate_coefs,
104 hypre_Vector *B)
105{
106 HYPRE_Int i, j;
107 HYPRE_Int ibeg, iend, pos;
108 HYPRE_Real a;
109
110 HYPRE_Int *Ai = hypre_CSRMatrixI(A);
111 HYPRE_Int *Aj = hypre_CSRMatrixJ(A);
112 HYPRE_Real *Adata = hypre_CSRMatrixData(A);
113 HYPRE_Int nrows = hypre_CSRMatrixNumRows(A);
114
115 HYPRE_Real *Bdata = hypre_VectorData(B);
116
117 for (i = 0; i < nrows; i++)
118 {
119 ibeg = Ai[i];
120 iend = Ai[i+1];
121 for (j = ibeg; j < iend; j++)
122 {
123 pos = hypre_BinarySearch(eliminate_cols, Aj[j], ncols_to_eliminate);
124 if (pos != -1)
125 {
126 a = Adata[j];
127 Adata[j] = 0.0;
128 Bdata[i] -= a * eliminate_coefs[pos];
129 }
130 }
131 }
132}
133
134/*
135 Function: hypre_CSRMatrixEliminateOffdRowsAXB
136
137 Eliminate (zero) the given list of rows of A.
138*/
139void hypre_CSRMatrixEliminateOffdRowsAXB(hypre_CSRMatrix *A,
140 HYPRE_Int nrows_to_eliminate,
141 HYPRE_Int *rows_to_eliminate)
142{
143 HYPRE_Int *Ai = hypre_CSRMatrixI(A);
144 HYPRE_Real *Adata = hypre_CSRMatrixData(A);
145
146 HYPRE_Int i, j;
147 HYPRE_Int irow, ibeg, iend;
148
149 for (i = 0; i < nrows_to_eliminate; i++)
150 {
151 irow = rows_to_eliminate[i];
152 ibeg = Ai[irow];
153 iend = Ai[irow+1];
154 for (j = ibeg; j < iend; j++)
155 {
156 Adata[j] = 0.0;
157 }
158 }
159}
160
161/*
162 Function: hypre_ParCSRMatrixEliminateAXB
163
164 This function eliminates the global rows and columns of a matrix
165 A corresponding to given lists of sorted (!) local row numbers,
166 so that the solution to the system A*X = B is X_b for the given rows.
167
168 The elimination is done as follows:
169
170 (input) (output)
171
172 / A_ii | A_ib \ / A_ii | 0 \
173 A = | -----+----- | ---> | -----+----- |
174 \ A_bi | A_bb / \ 0 | I /
175
176 / X_i \ / X_i \
177 X = | --- | ---> | --- | (no change)
178 \ X_b / \ X_b /
179
180 / B_i \ / B_i - A_ib * X_b \
181 B = | --- | ---> | ---------------- |
182 \ B_b / \ X_b /
183
184*/
185void hypre_ParCSRMatrixEliminateAXB(hypre_ParCSRMatrix *A,
186 HYPRE_Int num_rowscols_to_elim,
187 HYPRE_Int *rowscols_to_elim,
188 hypre_ParVector *X,
189 hypre_ParVector *B)
190{
191 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
192 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
193 HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows(diag);
194 HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
195
196 hypre_Vector *Xlocal = hypre_ParVectorLocalVector(X);
197 hypre_Vector *Blocal = hypre_ParVectorLocalVector(B);
198
199 HYPRE_Real *Bdata = hypre_VectorData(Blocal);
200 HYPRE_Real *Xdata = hypre_VectorData(Xlocal);
201
202 HYPRE_Int num_offd_cols_to_elim;
203 HYPRE_Int *offd_cols_to_elim;
204 HYPRE_Real *eliminate_coefs;
205
206 /* figure out which offd cols should be eliminated and with what coef */
207 hypre_ParCSRCommHandle *comm_handle;
208 hypre_ParCSRCommPkg *comm_pkg;
209 HYPRE_Int num_sends;
210 HYPRE_Int index, start;
211 HYPRE_Int i, j, k, irow;
212
213 HYPRE_Real *eliminate_row = mfem_hypre_CTAlloc_host(HYPRE_Real, diag_nrows);
214 HYPRE_Real *eliminate_col = mfem_hypre_CTAlloc_host(HYPRE_Real, offd_ncols);
215 HYPRE_Real *buf_data, coef;
216
217 /* make sure A has a communication package */
218 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
219 if (!comm_pkg)
220 {
221 hypre_MatvecCommPkgCreate(A);
222 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
223 }
224
225 /* HACK: rows that shouldn't be eliminated are marked with quiet NaN;
226 those that should are set to the boundary value from X; this is to
227 avoid sending complex type (int+double) or communicating twice. */
228 for (i = 0; i < diag_nrows; i++)
229 {
230 eliminate_row[i] = std::numeric_limits<HYPRE_Real>::quiet_NaN();
231 }
232 for (i = 0; i < num_rowscols_to_elim; i++)
233 {
234 irow = rowscols_to_elim[i];
235 eliminate_row[irow] = Xdata[irow];
236 }
237
238 /* use a Matvec communication pattern to find (in eliminate_col)
239 which of the local offd columns are to be eliminated */
240 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
241 buf_data = mfem_hypre_CTAlloc_host(
242 HYPRE_Real,
243 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
244 index = 0;
245 for (i = 0; i < num_sends; i++)
246 {
247 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
248 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
249 {
250 k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
251 buf_data[index++] = eliminate_row[k];
252 }
253 }
254 comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg,
255 buf_data, eliminate_col);
256
257 /* do sequential part of the elimination while stuff is getting sent */
258 hypre_CSRMatrixEliminateAXB(diag, num_rowscols_to_elim, rowscols_to_elim,
259 Xlocal, Blocal);
260
261 /* finish the communication */
262 hypre_ParCSRCommHandleDestroy(comm_handle);
263
264 /* received eliminate_col[], count offd columns to eliminate */
265 num_offd_cols_to_elim = 0;
266 for (i = 0; i < offd_ncols; i++)
267 {
268 coef = eliminate_col[i];
269 if (coef == coef) // test for NaN
270 {
271 num_offd_cols_to_elim++;
272 }
273 }
274
275 offd_cols_to_elim = mfem_hypre_CTAlloc_host(HYPRE_Int,
276 num_offd_cols_to_elim);
277 eliminate_coefs = mfem_hypre_CTAlloc_host(HYPRE_Real, num_offd_cols_to_elim);
278
279 /* get a list of offd column indices and coefs */
280 num_offd_cols_to_elim = 0;
281 for (i = 0; i < offd_ncols; i++)
282 {
283 coef = eliminate_col[i];
284 if (coef == coef) // test for NaN
285 {
286 offd_cols_to_elim[num_offd_cols_to_elim] = i;
287 eliminate_coefs[num_offd_cols_to_elim] = coef;
288 num_offd_cols_to_elim++;
289 }
290 }
291
292 mfem_hypre_TFree_host(buf_data);
293 mfem_hypre_TFree_host(eliminate_col);
294 mfem_hypre_TFree_host(eliminate_row);
295
296 /* eliminate the off-diagonal part */
297 hypre_CSRMatrixEliminateOffdColsAXB(offd, num_offd_cols_to_elim,
298 offd_cols_to_elim,
299 eliminate_coefs, Blocal);
300
301 hypre_CSRMatrixEliminateOffdRowsAXB(offd, num_rowscols_to_elim,
302 rowscols_to_elim);
303
304 /* set boundary values in the rhs */
305 for (int ii = 0; ii < num_rowscols_to_elim; ii++)
306 {
307 irow = rowscols_to_elim[ii];
308 Bdata[irow] = Xdata[irow];
309 }
310
311 mfem_hypre_TFree_host(offd_cols_to_elim);
312 mfem_hypre_TFree_host(eliminate_coefs);
313}
314
315
316/*--------------------------------------------------------------------------
317 * (A + Ae) style elimination
318 *--------------------------------------------------------------------------*/
319
320/*
321 Function: hypre_CSRMatrixElimCreate
322
323 Prepare the Ae matrix: count nnz, initialize I, allocate J and data.
324*/
325void hypre_CSRMatrixElimCreate(hypre_CSRMatrix *A,
326 hypre_CSRMatrix *Ae,
327 HYPRE_Int nrows, HYPRE_Int *rows,
328 HYPRE_Int ncols, HYPRE_Int *cols,
329 HYPRE_Int *col_mark)
330{
331 HYPRE_Int i, j, col;
332 HYPRE_Int A_beg, A_end;
333
334 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
335 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
336 HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
337
338 hypre_CSRMatrixI(Ae) = mfem_hypre_TAlloc_host(HYPRE_Int, A_rows+1);
339
340 HYPRE_Int *Ae_i = hypre_CSRMatrixI(Ae);
341 HYPRE_Int nnz = 0;
342
343 for (i = 0; i < A_rows; i++)
344 {
345 Ae_i[i] = nnz;
346
347 A_beg = A_i[i];
348 A_end = A_i[i+1];
349
350 if (hypre_BinarySearch(rows, i, nrows) >= 0)
351 {
352 /* full row */
353 nnz += A_end - A_beg;
354
355 if (col_mark)
356 {
357 for (j = A_beg; j < A_end; j++)
358 {
359 col_mark[A_j[j]] = 1;
360 }
361 }
362 }
363 else
364 {
365 /* count columns */
366 for (j = A_beg; j < A_end; j++)
367 {
368 col = A_j[j];
369 if (hypre_BinarySearch(cols, col, ncols) >= 0)
370 {
371 nnz++;
372 if (col_mark) { col_mark[col] = 1; }
373 }
374 }
375 }
376 }
377 Ae_i[A_rows] = nnz;
378
379 hypre_CSRMatrixJ(Ae) = mfem_hypre_TAlloc_host(HYPRE_Int, nnz);
380 hypre_CSRMatrixData(Ae) = mfem_hypre_TAlloc_host(HYPRE_Real, nnz);
381 hypre_CSRMatrixNumNonzeros(Ae) = nnz;
382#if MFEM_HYPRE_VERSION >= 21800
383 hypre_CSRMatrixMemoryLocation(Ae) = HYPRE_MEMORY_HOST;
384#endif
385}
386
387/*
388 Function: hypre_CSRMatrixEliminateRowsCols
389
390 Eliminate rows and columns of A, store eliminated values in Ae.
391 If 'diag' is nonzero, the eliminated diagonal of A is set to identity.
392 If 'col_remap' is not NULL it specifies renumbering of columns of Ae.
393*/
394void hypre_CSRMatrixEliminateRowsCols(hypre_CSRMatrix *A,
395 hypre_CSRMatrix *Ae,
396 HYPRE_Int nrows, HYPRE_Int *rows,
397 HYPRE_Int ncols, HYPRE_Int *cols,
398 int diag, HYPRE_Int* col_remap)
399{
400 HYPRE_Int i, j, k, col;
401 HYPRE_Int A_beg, Ae_beg, A_end;
402 HYPRE_Real a;
403
404 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
405 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
406 HYPRE_Real *A_data = hypre_CSRMatrixData(A);
407 HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
408
409 HYPRE_Int *Ae_i = hypre_CSRMatrixI(Ae);
410 HYPRE_Int *Ae_j = hypre_CSRMatrixJ(Ae);
411 HYPRE_Real *Ae_data = hypre_CSRMatrixData(Ae);
412
413 for (i = 0; i < A_rows; i++)
414 {
415 A_beg = A_i[i];
416 A_end = A_i[i+1];
417 Ae_beg = Ae_i[i];
418
419 if (hypre_BinarySearch(rows, i, nrows) >= 0)
420 {
421 /* eliminate row */
422 for (j = A_beg, k = Ae_beg; j < A_end; j++, k++)
423 {
424 col = A_j[j];
425 Ae_j[k] = col_remap ? col_remap[col] : col;
426 a = (diag && col == i) ? 1.0 : 0.0;
427 Ae_data[k] = A_data[j] - a;
428 A_data[j] = a;
429 }
430 }
431 else
432 {
433 /* eliminate columns */
434 for (j = A_beg, k = Ae_beg; j < A_end; j++)
435 {
436 col = A_j[j];
437 if (hypre_BinarySearch(cols, col, ncols) >= 0)
438 {
439 Ae_j[k] = col_remap ? col_remap[col] : col;
440 Ae_data[k] = A_data[j];
441 A_data[j] = 0.0;
442 k++;
443 }
444 }
445 }
446 }
447}
448
449/*
450 Eliminate rows of A, setting all entries in the eliminated rows to zero.
451*/
452void hypre_CSRMatrixEliminateRows(hypre_CSRMatrix *A,
453 HYPRE_Int nrows, const HYPRE_Int *rows)
454{
455 HYPRE_Int irow, i, j;
456 HYPRE_Int A_beg, A_end;
457
458 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
459 HYPRE_Real *A_data = hypre_CSRMatrixData(A);
460
461 for (i = 0; i < nrows; i++)
462 {
463 irow = rows[i];
464 A_beg = A_i[irow];
465 A_end = A_i[irow+1];
466 /* eliminate row */
467 for (j = A_beg; j < A_end; j++)
468 {
469 A_data[j] = 0.0;
470 }
471 }
472}
473
474
475/*
476 Function: hypre_ParCSRMatrixEliminateAAe
477
478 (input) (output)
479
480 / A_ii | A_ib \ / A_ii | 0 \
481 A = | -----+----- | ---> | -----+----- |
482 \ A_bi | A_bb / \ 0 | I /
483
484
485 / 0 | A_ib \
486 Ae = | -----+--------- |
487 \ A_bi | A_bb - I /
488
489*/
490
491void hypre_ParCSRMatrixEliminateAAe(hypre_ParCSRMatrix *A,
492 hypre_ParCSRMatrix **Ae,
493 HYPRE_Int num_rowscols_to_elim,
494 HYPRE_Int *rowscols_to_elim,
495 int ignore_rows)
496{
497 HYPRE_Int i, j, k;
498
499 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
500 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
501 HYPRE_Int A_diag_ncols = hypre_CSRMatrixNumCols(A_diag);
502 HYPRE_Int A_offd_ncols = hypre_CSRMatrixNumCols(A_offd);
503
504 *Ae = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
505 hypre_ParCSRMatrixGlobalNumRows(A),
506 hypre_ParCSRMatrixGlobalNumCols(A),
507 hypre_ParCSRMatrixRowStarts(A),
508 hypre_ParCSRMatrixColStarts(A),
509 0, 0, 0);
510
511#if MFEM_HYPRE_VERSION <= 22200
512 hypre_ParCSRMatrixSetRowStartsOwner(*Ae, 0);
513 hypre_ParCSRMatrixSetColStartsOwner(*Ae, 0);
514#endif
515
516 hypre_CSRMatrix *Ae_diag = hypre_ParCSRMatrixDiag(*Ae);
517 hypre_CSRMatrix *Ae_offd = hypre_ParCSRMatrixOffd(*Ae);
518 HYPRE_Int Ae_offd_ncols;
519
520 HYPRE_Int num_offd_cols_to_elim;
521 HYPRE_Int *offd_cols_to_elim;
522
523 HYPRE_BigInt *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
524 HYPRE_BigInt *Ae_col_map_offd;
525
526 HYPRE_Int *col_mark;
527 HYPRE_Int *col_remap;
528
529 /* figure out which offd cols should be eliminated */
530 {
531 hypre_ParCSRCommHandle *comm_handle;
532 hypre_ParCSRCommPkg *comm_pkg;
533 HYPRE_Int num_sends, *int_buf_data;
534 HYPRE_Int index, start;
535
536 HYPRE_Int *eliminate_diag_col = mfem_hypre_CTAlloc_host(HYPRE_Int,
537 A_diag_ncols);
538 HYPRE_Int *eliminate_offd_col = mfem_hypre_CTAlloc_host(HYPRE_Int,
539 A_offd_ncols);
540
541 /* make sure A has a communication package */
542 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
543 if (!comm_pkg)
544 {
545 hypre_MatvecCommPkgCreate(A);
546 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
547 }
548
549 /* which of the local rows are to be eliminated */
550 for (i = 0; i < A_diag_ncols; i++)
551 {
552 eliminate_diag_col[i] = 0;
553 }
554 for (i = 0; i < num_rowscols_to_elim; i++)
555 {
556 eliminate_diag_col[rowscols_to_elim[i]] = 1;
557 }
558
559 /* use a Matvec communication pattern to find (in eliminate_col)
560 which of the local offd columns are to be eliminated */
561 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
562 int_buf_data = mfem_hypre_CTAlloc_host(
563 HYPRE_Int,
564 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
565 index = 0;
566 for (i = 0; i < num_sends; i++)
567 {
568 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
569 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
570 {
571 k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
572 int_buf_data[index++] = eliminate_diag_col[k];
573 }
574 }
575 comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
576 int_buf_data,
577 eliminate_offd_col);
578
579 /* eliminate diagonal part, overlapping it with communication */
580 if (ignore_rows)
581 {
582 hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
583 0, nullptr,
584 num_rowscols_to_elim, rowscols_to_elim,
585 NULL);
586
587 hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
588 0, nullptr,
589 num_rowscols_to_elim, rowscols_to_elim,
590 1, NULL);
591 }
592 else
593 {
594 hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
595 num_rowscols_to_elim, rowscols_to_elim,
596 num_rowscols_to_elim, rowscols_to_elim,
597 NULL);
598
599 hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
600 num_rowscols_to_elim, rowscols_to_elim,
601 num_rowscols_to_elim, rowscols_to_elim,
602 1, NULL);
603 }
604
605 hypre_CSRMatrixReorder(Ae_diag);
606
607 /* finish the communication */
608 hypre_ParCSRCommHandleDestroy(comm_handle);
609
610 /* received eliminate_col[], count offd columns to eliminate */
611 num_offd_cols_to_elim = 0;
612 for (i = 0; i < A_offd_ncols; i++)
613 {
614 if (eliminate_offd_col[i]) { num_offd_cols_to_elim++; }
615 }
616
617 offd_cols_to_elim = mfem_hypre_CTAlloc_host(HYPRE_Int,
618 num_offd_cols_to_elim);
619
620 /* get a list of offd column indices and coefs */
621 num_offd_cols_to_elim = 0;
622 for (i = 0; i < A_offd_ncols; i++)
623 {
624 if (eliminate_offd_col[i])
625 {
626 offd_cols_to_elim[num_offd_cols_to_elim++] = i;
627 }
628 }
629
630 mfem_hypre_TFree_host(int_buf_data);
631 mfem_hypre_TFree_host(eliminate_offd_col);
632 mfem_hypre_TFree_host(eliminate_diag_col);
633 }
634
635 /* eliminate the off-diagonal part */
636 col_mark = mfem_hypre_CTAlloc_host(HYPRE_Int, A_offd_ncols);
637 col_remap = mfem_hypre_CTAlloc_host(HYPRE_Int, A_offd_ncols);
638
639 if (ignore_rows)
640 {
641 hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
642 0, nullptr,
643 num_offd_cols_to_elim, offd_cols_to_elim,
644 col_mark);
645
646 for (i = k = 0; i < A_offd_ncols; i++)
647 {
648 if (col_mark[i]) { col_remap[i] = k++; }
649 }
650
651 hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
652 0, nullptr,
653 num_offd_cols_to_elim, offd_cols_to_elim,
654 0, col_remap);
655 }
656 else
657 {
658 hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
659 num_rowscols_to_elim, rowscols_to_elim,
660 num_offd_cols_to_elim, offd_cols_to_elim,
661 col_mark);
662
663 for (i = k = 0; i < A_offd_ncols; i++)
664 {
665 if (col_mark[i]) { col_remap[i] = k++; }
666 }
667
668 hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
669 num_rowscols_to_elim, rowscols_to_elim,
670 num_offd_cols_to_elim, offd_cols_to_elim,
671 0, col_remap);
672 }
673
674 /* create col_map_offd for Ae */
675 Ae_offd_ncols = 0;
676 for (i = 0; i < A_offd_ncols; i++)
677 {
678 if (col_mark[i]) { Ae_offd_ncols++; }
679 }
680
681 Ae_col_map_offd = mfem_hypre_CTAlloc_host(HYPRE_BigInt, Ae_offd_ncols);
682
683 Ae_offd_ncols = 0;
684 for (i = 0; i < A_offd_ncols; i++)
685 {
686 if (col_mark[i])
687 {
688 Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
689 }
690 }
691
692 hypre_ParCSRMatrixColMapOffd(*Ae) = Ae_col_map_offd;
693 hypre_CSRMatrixNumCols(Ae_offd) = Ae_offd_ncols;
694
695 mfem_hypre_TFree_host(col_remap);
696 mfem_hypre_TFree_host(col_mark);
697 mfem_hypre_TFree_host(offd_cols_to_elim);
698
699 hypre_ParCSRMatrixSetNumNonzeros(*Ae);
700 hypre_MatvecCommPkgCreate(*Ae);
701}
702
703
704// Eliminate rows from the diagonal and off-diagonal blocks of the matrix
705void hypre_ParCSRMatrixEliminateRows(hypre_ParCSRMatrix *A,
706 HYPRE_Int num_rows_to_elim,
707 const HYPRE_Int *rows_to_elim)
708{
709 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
710 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
711 hypre_CSRMatrixEliminateRows(A_diag, num_rows_to_elim, rows_to_elim);
712 hypre_CSRMatrixEliminateRows(A_offd, num_rows_to_elim, rows_to_elim);
713}
714
715
716/*--------------------------------------------------------------------------
717 * Split
718 *--------------------------------------------------------------------------*/
719
720void hypre_CSRMatrixSplit(hypre_CSRMatrix *A,
721 HYPRE_Int nr, HYPRE_Int nc,
722 HYPRE_Int *row_block_num, HYPRE_Int *col_block_num,
723 hypre_CSRMatrix **blocks)
724{
725 HYPRE_Int i, j, k, bi, bj;
726
727 HYPRE_Int* A_i = hypre_CSRMatrixI(A);
728 HYPRE_Int* A_j = hypre_CSRMatrixJ(A);
729 HYPRE_Complex* A_data = hypre_CSRMatrixData(A);
730
731 HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
732 HYPRE_Int A_cols = hypre_CSRMatrixNumCols(A);
733
734 HYPRE_Int *num_rows = mfem_hypre_CTAlloc_host(HYPRE_Int, nr);
735 HYPRE_Int *num_cols = mfem_hypre_CTAlloc_host(HYPRE_Int, nc);
736
737 HYPRE_Int *block_row = mfem_hypre_TAlloc_host(HYPRE_Int, A_rows);
738 HYPRE_Int *block_col = mfem_hypre_TAlloc_host(HYPRE_Int, A_cols);
739
740 for (i = 0; i < A_rows; i++)
741 {
742 block_row[i] = num_rows[row_block_num[i]]++;
743 }
744 for (j = 0; j < A_cols; j++)
745 {
746 block_col[j] = num_cols[col_block_num[j]]++;
747 }
748
749 /* allocate the blocks */
750 for (i = 0; i < nr; i++)
751 {
752 for (j = 0; j < nc; j++)
753 {
754 hypre_CSRMatrix *B = hypre_CSRMatrixCreate(num_rows[i],
755 num_cols[j], 0);
756 hypre_CSRMatrixI(B) = mfem_hypre_CTAlloc_host(HYPRE_Int,
757 num_rows[i] + 1);
758#if MFEM_HYPRE_VERSION >= 21800
759 hypre_CSRMatrixMemoryLocation(B) = HYPRE_MEMORY_HOST;
760#endif
761 blocks[i*nc + j] = B;
762 }
763 }
764
765 /* count block row nnz */
766 for (i = 0; i < A_rows; i++)
767 {
768 bi = row_block_num[i];
769 for (j = A_i[i]; j < A_i[i+1]; j++)
770 {
771 bj = col_block_num[A_j[j]];
772 hypre_CSRMatrix *B = blocks[bi*nc + bj];
773 hypre_CSRMatrixI(B)[block_row[i] + 1]++;
774 }
775 }
776
777 /* count block nnz */
778 for (k = 0; k < nr*nc; k++)
779 {
780 hypre_CSRMatrix *B = blocks[k];
781 HYPRE_Int* B_i = hypre_CSRMatrixI(B);
782
783 HYPRE_Int nnz = 0, rs;
784 for (int kk = 1; kk <= hypre_CSRMatrixNumRows(B); kk++)
785 {
786 rs = B_i[kk], B_i[kk] = nnz, nnz += rs;
787 }
788
789 hypre_CSRMatrixJ(B) = mfem_hypre_TAlloc_host(HYPRE_Int, nnz);
790 hypre_CSRMatrixData(B) = mfem_hypre_TAlloc_host(HYPRE_Complex, nnz);
791 hypre_CSRMatrixNumNonzeros(B) = nnz;
792 }
793
794 /* populate blocks */
795 for (i = 0; i < A_rows; i++)
796 {
797 bi = row_block_num[i];
798 for (j = A_i[i]; j < A_i[i+1]; j++)
799 {
800 k = A_j[j];
801 bj = col_block_num[k];
802 hypre_CSRMatrix *B = blocks[bi*nc + bj];
803 HYPRE_Int *bii = hypre_CSRMatrixI(B) + block_row[i] + 1;
804 hypre_CSRMatrixJ(B)[*bii] = block_col[k];
805 hypre_CSRMatrixData(B)[*bii] = A_data[j];
806 (*bii)++;
807 }
808 }
809
810 mfem_hypre_TFree_host(block_col);
811 mfem_hypre_TFree_host(block_row);
812
813 mfem_hypre_TFree_host(num_cols);
814 mfem_hypre_TFree_host(num_rows);
815}
816
817
818void hypre_ParCSRMatrixSplit(hypre_ParCSRMatrix *A,
819 HYPRE_Int nr, HYPRE_Int nc,
820 hypre_ParCSRMatrix **blocks,
821 int interleaved_rows, int interleaved_cols)
822{
823 HYPRE_Int i, j, k;
824
825 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
826
827 hypre_CSRMatrix *Adiag = hypre_ParCSRMatrixDiag(A);
828 hypre_CSRMatrix *Aoffd = hypre_ParCSRMatrixOffd(A);
829
830 HYPRE_BigInt global_rows = hypre_ParCSRMatrixGlobalNumRows(A);
831 HYPRE_BigInt global_cols = hypre_ParCSRMatrixGlobalNumCols(A);
832
833 HYPRE_Int local_rows = hypre_CSRMatrixNumRows(Adiag);
834 HYPRE_Int local_cols = hypre_CSRMatrixNumCols(Adiag);
835 HYPRE_Int offd_cols = hypre_CSRMatrixNumCols(Aoffd);
836
837 hypre_assert(local_rows % nr == 0 && local_cols % nc == 0);
838 hypre_assert(global_rows % nr == 0 && global_cols % nc == 0);
839
840 HYPRE_Int block_rows = local_rows / nr;
841 HYPRE_Int block_cols = local_cols / nc;
842 HYPRE_Int num_blocks = nr * nc;
843
844 /* mark local rows and columns with block number */
845 HYPRE_Int *row_block_num = mfem_hypre_TAlloc_host(HYPRE_Int, local_rows);
846 HYPRE_Int *col_block_num = mfem_hypre_TAlloc_host(HYPRE_Int, local_cols);
847
848 for (i = 0; i < local_rows; i++)
849 {
850 row_block_num[i] = interleaved_rows ? (i % nr) : (i / block_rows);
851 }
852 for (i = 0; i < local_cols; i++)
853 {
854 col_block_num[i] = interleaved_cols ? (i % nc) : (i / block_cols);
855 }
856
857 /* determine the block numbers for offd columns */
858 HYPRE_BigInt *offd_col_block_num = mfem_hypre_TAlloc_host(HYPRE_BigInt,
859 offd_cols);
860 hypre_ParCSRCommHandle *comm_handle;
861 HYPRE_BigInt *int_buf_data;
862 {
863 /* make sure A has a communication package */
864 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
865 if (!comm_pkg)
866 {
867 hypre_MatvecCommPkgCreate(A);
868 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
869 }
870
871 /* calculate the final global column numbers for each block */
872 HYPRE_Int *count = mfem_hypre_CTAlloc_host(HYPRE_Int, nc);
873 HYPRE_BigInt *block_global_col = mfem_hypre_TAlloc_host(HYPRE_BigInt,
874 local_cols);
875 HYPRE_BigInt first_col = hypre_ParCSRMatrixFirstColDiag(A) / nc;
876 for (i = 0; i < local_cols; i++)
877 {
878 block_global_col[i] = first_col + count[col_block_num[i]]++;
879 }
880 mfem_hypre_TFree_host(count);
881
882 /* use a Matvec communication pattern to determine offd_col_block_num */
883 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
884 int_buf_data = mfem_hypre_CTAlloc_host(
885 HYPRE_BigInt,
886 hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
887 HYPRE_Int start, index = 0;
888 for (i = 0; i < num_sends; i++)
889 {
890 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
891 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
892 {
893 k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
894 int_buf_data[index++] = col_block_num[k] + nc*block_global_col[k];
895 }
896 }
897
898 mfem_hypre_TFree_host(block_global_col);
899
900#if MFEM_HYPRE_VERSION < 21600
901 const int job = 11;
902#else
903 const int job = 21;
904#endif
905 comm_handle = hypre_ParCSRCommHandleCreate(job, comm_pkg, int_buf_data,
906 offd_col_block_num);
907 }
908
909 /* create the block matrices */
910 HYPRE_Int num_procs = 1;
911 if (!hypre_ParCSRMatrixAssumedPartition(A))
912 {
913 hypre_MPI_Comm_size(comm, &num_procs);
914 }
915
916 HYPRE_BigInt *row_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, num_procs+1);
917 HYPRE_BigInt *col_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, num_procs+1);
918 for (i = 0; i <= num_procs; i++)
919 {
920 row_starts[i] = hypre_ParCSRMatrixRowStarts(A)[i] / nr;
921 col_starts[i] = hypre_ParCSRMatrixColStarts(A)[i] / nc;
922 }
923
924 for (i = 0; i < num_blocks; i++)
925 {
926 blocks[i] = hypre_ParCSRMatrixCreate(comm,
927 global_rows/nr, global_cols/nc,
928 row_starts, col_starts, 0, 0, 0);
929 }
930
931 /* split diag part */
932 hypre_CSRMatrix **csr_blocks = mfem_hypre_TAlloc_host(hypre_CSRMatrix*,
933 nr*nc);
934 hypre_CSRMatrixSplit(Adiag, nr, nc, row_block_num, col_block_num,
935 csr_blocks);
936
937 for (i = 0; i < num_blocks; i++)
938 {
939 mfem_hypre_TFree_host(hypre_ParCSRMatrixDiag(blocks[i]));
940 hypre_ParCSRMatrixDiag(blocks[i]) = csr_blocks[i];
941 }
942
943 /* finish communication, receive offd_col_block_num */
944 hypre_ParCSRCommHandleDestroy(comm_handle);
945 mfem_hypre_TFree_host(int_buf_data);
946
947 /* decode global offd column numbers */
948 HYPRE_Int *offd_col_block_num_nc = mfem_hypre_TAlloc_host(HYPRE_Int,
949 offd_cols);
950 HYPRE_BigInt* offd_global_col = mfem_hypre_TAlloc_host(HYPRE_BigInt,
951 offd_cols);
952 for (i = 0; i < offd_cols; i++)
953 {
954 offd_global_col[i] = offd_col_block_num[i] / nc;
955 offd_col_block_num_nc[i] = offd_col_block_num[i] % nc;
956 }
957
958 mfem_hypre_TFree_host(offd_col_block_num);
959
960 /* split offd part */
961 hypre_CSRMatrixSplit(Aoffd, nr, nc, row_block_num, offd_col_block_num_nc,
962 csr_blocks);
963
964 for (i = 0; i < num_blocks; i++)
965 {
966 mfem_hypre_TFree_host(hypre_ParCSRMatrixOffd(blocks[i]));
967 hypre_ParCSRMatrixOffd(blocks[i]) = csr_blocks[i];
968 }
969
970 mfem_hypre_TFree_host(csr_blocks);
971 mfem_hypre_TFree_host(col_block_num);
972 mfem_hypre_TFree_host(row_block_num);
973
974 /* update block col-maps */
975 for (int bi = 0; bi < nr; bi++)
976 {
977 for (int bj = 0; bj < nc; bj++)
978 {
979 hypre_ParCSRMatrix *block = blocks[bi*nc + bj];
980 hypre_CSRMatrix *block_offd = hypre_ParCSRMatrixOffd(block);
981 HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols(block_offd);
982
983 HYPRE_BigInt *block_col_map = mfem_hypre_TAlloc_host(HYPRE_BigInt,
984 block_offd_cols);
985 for (i = j = 0; i < offd_cols; i++)
986 {
987 HYPRE_Int bn = offd_col_block_num_nc[i];
988 if (bn == bj) { block_col_map[j++] = offd_global_col[i]; }
989 }
990 hypre_assert(j == block_offd_cols);
991
992 hypre_ParCSRMatrixColMapOffd(block) = block_col_map;
993 }
994 }
995
996 mfem_hypre_TFree_host(offd_global_col);
997 mfem_hypre_TFree_host(offd_col_block_num_nc);
998
999 /* finish the new matrices, make them own all the stuff */
1000 for (i = 0; i < num_blocks; i++)
1001 {
1002 hypre_ParCSRMatrixSetNumNonzeros(blocks[i]);
1003 hypre_MatvecCommPkgCreate(blocks[i]);
1004
1005 hypre_ParCSRMatrixOwnsData(blocks[i]) = 1;
1006
1007#if MFEM_HYPRE_VERSION <= 22200
1008 /* only the first block will own the row/col_starts */
1009 hypre_ParCSRMatrixOwnsRowStarts(blocks[i]) = !i;
1010 hypre_ParCSRMatrixOwnsColStarts(blocks[i]) = !i;
1011#endif
1012 }
1013
1014#if MFEM_HYPRE_VERSION > 22200
1015 mfem_hypre_TFree_host(row_starts);
1016 mfem_hypre_TFree_host(col_starts);
1017#endif
1018}
1019
1020/* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
1021void hypre_CSRMatrixAbsMatvec(hypre_CSRMatrix *A,
1022 HYPRE_Real alpha,
1023 HYPRE_Real *x,
1024 HYPRE_Real beta,
1025 HYPRE_Real *y)
1026{
1027 HYPRE_Real *A_data = hypre_CSRMatrixData(A);
1028 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1029 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1030 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1031
1032 HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
1033 HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
1034
1035 HYPRE_Real *x_data = x;
1036 HYPRE_Real *y_data = y;
1037
1038 HYPRE_Real temp, tempx;
1039
1040 HYPRE_Int i, jj;
1041
1042 HYPRE_Int m;
1043
1044 HYPRE_Real xpar=0.7;
1045
1046 /*-----------------------------------------------------------------------
1047 * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
1048 *-----------------------------------------------------------------------*/
1049
1050 if (alpha == 0.0)
1051 {
1052 for (i = 0; i < num_rows; i++)
1053 {
1054 y_data[i] *= beta;
1055 }
1056 return;
1057 }
1058
1059 /*-----------------------------------------------------------------------
1060 * y = (beta/alpha)*y
1061 *-----------------------------------------------------------------------*/
1062
1063 temp = beta / alpha;
1064
1065 if (temp != 1.0)
1066 {
1067 if (temp == 0.0)
1068 {
1069 for (i = 0; i < num_rows; i++)
1070 {
1071 y_data[i] = 0.0;
1072 }
1073 }
1074 else
1075 {
1076 for (i = 0; i < num_rows; i++)
1077 {
1078 y_data[i] *= temp;
1079 }
1080 }
1081 }
1082
1083 /*-----------------------------------------------------------------
1084 * y += abs(A)*x
1085 *-----------------------------------------------------------------*/
1086
1087 /* use rownnz pointer to do the abs(A)*x multiplication
1088 when num_rownnz is smaller than num_rows */
1089
1090 if (num_rownnz < xpar*(num_rows))
1091 {
1092 for (i = 0; i < num_rownnz; i++)
1093 {
1094 m = A_rownnz[i];
1095
1096 tempx = 0;
1097 for (jj = A_i[m]; jj < A_i[m+1]; jj++)
1098 {
1099 tempx += std::abs(A_data[jj])*x_data[A_j[jj]];
1100 }
1101 y_data[m] += tempx;
1102 }
1103 }
1104 else
1105 {
1106 for (i = 0; i < num_rows; i++)
1107 {
1108 tempx = 0;
1109 for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1110 {
1111 tempx += std::abs(A_data[jj])*x_data[A_j[jj]];
1112 }
1113 y_data[i] += tempx;
1114 }
1115 }
1116
1117 /*-----------------------------------------------------------------
1118 * y = alpha*y
1119 *-----------------------------------------------------------------*/
1120
1121 if (alpha != 1.0)
1122 {
1123 for (i = 0; i < num_rows; i++)
1124 {
1125 y_data[i] *= alpha;
1126 }
1127 }
1128
1129}
1130
1131/* Based on hypre_CSRMatrixMatvecT in hypre's csr_matvec.c */
1132void hypre_CSRMatrixAbsMatvecT(hypre_CSRMatrix *A,
1133 HYPRE_Real alpha,
1134 HYPRE_Real *x,
1135 HYPRE_Real beta,
1136 HYPRE_Real *y)
1137{
1138 HYPRE_Real *A_data = hypre_CSRMatrixData(A);
1139 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1140 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1141 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1142 HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
1143
1144 HYPRE_Real *x_data = x;
1145 HYPRE_Real *y_data = y;
1146
1147 HYPRE_Int i, j, jj;
1148
1149 HYPRE_Real temp;
1150
1151 if (alpha == 0.0)
1152 {
1153 for (i = 0; i < num_cols; i++)
1154 {
1155 y_data[i] *= beta;
1156 }
1157 return;
1158 }
1159
1160 /*-----------------------------------------------------------------------
1161 * y = (beta/alpha)*y
1162 *-----------------------------------------------------------------------*/
1163
1164 temp = beta / alpha;
1165
1166 if (temp != 1.0)
1167 {
1168 if (temp == 0.0)
1169 {
1170 for (i = 0; i < num_cols; i++)
1171 {
1172 y_data[i] = 0.0;
1173 }
1174 }
1175 else
1176 {
1177 for (i = 0; i < num_cols; i++)
1178 {
1179 y_data[i] *= temp;
1180 }
1181 }
1182 }
1183
1184 /*-----------------------------------------------------------------
1185 * y += abs(A)^T*x
1186 *-----------------------------------------------------------------*/
1187
1188 for (i = 0; i < num_rows; i++)
1189 {
1190 for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1191 {
1192 j = A_j[jj];
1193 y_data[j] += std::abs(A_data[jj]) * x_data[i];
1194 }
1195 }
1196
1197 /*-----------------------------------------------------------------
1198 * y = alpha*y
1199 *-----------------------------------------------------------------*/
1200
1201 if (alpha != 1.0)
1202 {
1203 for (i = 0; i < num_cols; i++)
1204 {
1205 y_data[i] *= alpha;
1206 }
1207 }
1208}
1209
1210/* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
1211void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,
1212 HYPRE_Bool alpha,
1213 HYPRE_Bool *x,
1214 HYPRE_Bool beta,
1215 HYPRE_Bool *y)
1216{
1217 /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1218 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1219 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1220 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1221
1222 HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
1223 HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
1224
1225 HYPRE_Bool *x_data = x;
1226 HYPRE_Bool *y_data = y;
1227
1228 HYPRE_Bool temp, tempx;
1229
1230 HYPRE_Int i, jj;
1231
1232 HYPRE_Int m;
1233
1234 HYPRE_Real xpar=0.7;
1235
1236 /*-----------------------------------------------------------------------
1237 * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
1238 *-----------------------------------------------------------------------*/
1239
1240 if (alpha == 0)
1241 {
1242#ifdef HYPRE_USING_OPENMP
1243 #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1244#endif
1245 for (i = 0; i < num_rows; i++)
1246 {
1247 y_data[i] = y_data[i] && beta;
1248 }
1249 return;
1250 }
1251
1252 /*-----------------------------------------------------------------------
1253 * y = (beta/alpha)*y
1254 *-----------------------------------------------------------------------*/
1255
1256 if (beta == 0)
1257 {
1258#ifdef HYPRE_USING_OPENMP
1259 #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1260#endif
1261 for (i = 0; i < num_rows; i++)
1262 {
1263 y_data[i] = 0;
1264 }
1265 }
1266 else
1267 {
1268 /* beta is true -> no change to y_data */
1269 }
1270
1271 /*-----------------------------------------------------------------
1272 * y += A*x
1273 *-----------------------------------------------------------------*/
1274
1275 /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows */
1276
1277 if (num_rownnz < xpar*(num_rows))
1278 {
1279#ifdef HYPRE_USING_OPENMP
1280 #pragma omp parallel for private(i,jj,m,tempx) HYPRE_SMP_SCHEDULE
1281#endif
1282 for (i = 0; i < num_rownnz; i++)
1283 {
1284 m = A_rownnz[i];
1285
1286 tempx = 0;
1287 for (jj = A_i[m]; jj < A_i[m+1]; jj++)
1288 {
1289 /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1290 tempx = tempx || x_data[A_j[jj]];
1291 }
1292 y_data[m] = y_data[m] || tempx;
1293 }
1294 }
1295 else
1296 {
1297#ifdef HYPRE_USING_OPENMP
1298 #pragma omp parallel for private(i,jj,temp) HYPRE_SMP_SCHEDULE
1299#endif
1300 for (i = 0; i < num_rows; i++)
1301 {
1302 temp = 0;
1303 for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1304 {
1305 /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1306 temp = temp || x_data[A_j[jj]];
1307 }
1308 y_data[i] = y_data[i] || temp;
1309 }
1310 }
1311
1312 /*-----------------------------------------------------------------
1313 * y = alpha*y
1314 *-----------------------------------------------------------------*/
1315 /* alpha is true */
1316}
1317
1318/* Based on hypre_CSRMatrixMatvecT in hypre's csr_matvec.c */
1319void hypre_CSRMatrixBooleanMatvecT(hypre_CSRMatrix *A,
1320 HYPRE_Bool alpha,
1321 HYPRE_Bool *x,
1322 HYPRE_Bool beta,
1323 HYPRE_Bool *y)
1324{
1325 /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1326 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1327 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1328 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1329 HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
1330
1331 HYPRE_Bool *x_data = x;
1332 HYPRE_Bool *y_data = y;
1333
1334 HYPRE_Int i, j, jj;
1335
1336 /*-----------------------------------------------------------------------
1337 * y = beta*y
1338 *-----------------------------------------------------------------------*/
1339
1340 if (beta == 0)
1341 {
1342 for (i = 0; i < num_cols; i++)
1343 {
1344 y_data[i] = 0;
1345 }
1346 }
1347 else
1348 {
1349 /* beta is true -> no change to y_data */
1350 }
1351
1352 /*-----------------------------------------------------------------------
1353 * Check if (alpha == 0)
1354 *-----------------------------------------------------------------------*/
1355
1356 if (alpha == 0)
1357 {
1358 return;
1359 }
1360
1361 /* alpha is true */
1362
1363 /*-----------------------------------------------------------------
1364 * y += A^T*x
1365 *-----------------------------------------------------------------*/
1366 for (i = 0; i < num_rows; i++)
1367 {
1368 if (x_data[i] != 0)
1369 {
1370 for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1371 {
1372 j = A_j[jj];
1373 /* y_data[j] += A_data[jj] * x_data[i]; */
1374 y_data[j] = 1;
1375 }
1376 }
1377 }
1378}
1379
1380/* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
1381 input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
1382hypre_ParCSRCommHandle *
1383hypre_ParCSRCommHandleCreate_bool(HYPRE_Int job,
1384 hypre_ParCSRCommPkg *comm_pkg,
1385 HYPRE_Bool *send_data,
1386 HYPRE_Bool *recv_data)
1387{
1388 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1389 HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
1390 MPI_Comm comm = hypre_ParCSRCommPkgComm(comm_pkg);
1391
1392 hypre_ParCSRCommHandle *comm_handle;
1393 HYPRE_Int num_requests;
1394 hypre_MPI_Request *requests;
1395
1396 HYPRE_Int i, j;
1397 HYPRE_Int my_id, num_procs;
1398 HYPRE_Int ip, vec_start, vec_len;
1399
1400 num_requests = num_sends + num_recvs;
1401 requests = mfem_hypre_CTAlloc_host(hypre_MPI_Request, num_requests);
1402
1403 hypre_MPI_Comm_size(comm, &num_procs);
1404 hypre_MPI_Comm_rank(comm, &my_id);
1405
1406 j = 0;
1407 switch (job)
1408 {
1409 case 1:
1410 {
1411 HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1412 HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1413 for (i = 0; i < num_recvs; i++)
1414 {
1415 ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1416 vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1417 vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1418 hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1419 ip, 0, comm, &requests[j++]);
1420 }
1421 for (i = 0; i < num_sends; i++)
1422 {
1423 vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1424 vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1425 ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1426 hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1427 ip, 0, comm, &requests[j++]);
1428 }
1429 break;
1430 }
1431 case 2:
1432 {
1433 HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1434 HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1435 for (i = 0; i < num_sends; i++)
1436 {
1437 vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1438 vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1439 ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1440 hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1441 ip, 0, comm, &requests[j++]);
1442 }
1443 for (i = 0; i < num_recvs; i++)
1444 {
1445 ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1446 vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1447 vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1448 hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1449 ip, 0, comm, &requests[j++]);
1450 }
1451 break;
1452 }
1453 }
1454 /*--------------------------------------------------------------------
1455 * set up comm_handle and return
1456 *--------------------------------------------------------------------*/
1457
1458 comm_handle = mfem_hypre_CTAlloc_host(hypre_ParCSRCommHandle, 1);
1459
1460 hypre_ParCSRCommHandleCommPkg(comm_handle) = comm_pkg;
1461 hypre_ParCSRCommHandleSendData(comm_handle) = send_data;
1462 hypre_ParCSRCommHandleRecvData(comm_handle) = recv_data;
1463 hypre_ParCSRCommHandleNumRequests(comm_handle) = num_requests;
1464 hypre_ParCSRCommHandleRequests(comm_handle) = requests;
1465
1466 return comm_handle;
1467}
1468
1469/* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1470void hypre_ParCSRMatrixAbsMatvec(hypre_ParCSRMatrix *A,
1471 HYPRE_Real alpha,
1472 HYPRE_Real *x,
1473 HYPRE_Real beta,
1474 HYPRE_Real *y)
1475{
1476 hypre_ParCSRCommHandle *comm_handle;
1477 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1478 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1479 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1480
1481 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1482 HYPRE_Int num_sends, i, j, index;
1483
1484 HYPRE_Real *x_tmp, *x_buf;
1485
1486 x_tmp = mfem_hypre_CTAlloc_host(HYPRE_Real, num_cols_offd);
1487
1488 /*---------------------------------------------------------------------
1489 * If there exists no CommPkg for A, a CommPkg is generated using
1490 * equally load balanced partitionings
1491 *--------------------------------------------------------------------*/
1492 if (!comm_pkg)
1493 {
1494 hypre_MatvecCommPkgCreate(A);
1495 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1496 }
1497
1498 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1499 x_buf = mfem_hypre_CTAlloc_host(
1500 HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1501
1502 index = 0;
1503 for (i = 0; i < num_sends; i++)
1504 {
1505 j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1506 for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1507 {
1508 x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1509 }
1510 }
1511
1512 comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg, x_buf, x_tmp);
1513
1514 hypre_CSRMatrixAbsMatvec(diag, alpha, x, beta, y);
1515
1516 hypre_ParCSRCommHandleDestroy(comm_handle);
1517
1518 if (num_cols_offd)
1519 {
1520 hypre_CSRMatrixAbsMatvec(offd, alpha, x_tmp, 1.0, y);
1521 }
1522
1523 mfem_hypre_TFree_host(x_buf);
1524 mfem_hypre_TFree_host(x_tmp);
1525}
1526
1527/* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
1528void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
1529 HYPRE_Real alpha,
1530 HYPRE_Real *x,
1531 HYPRE_Real beta,
1532 HYPRE_Real *y)
1533{
1534 hypre_ParCSRCommHandle *comm_handle;
1535 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1536 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1537 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1538 HYPRE_Real *y_tmp;
1539 HYPRE_Real *y_buf;
1540
1541 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1542
1543 HYPRE_Int i, j, jj, end, num_sends;
1544
1545 y_tmp = mfem_hypre_TAlloc_host(HYPRE_Real, num_cols_offd);
1546
1547 /*---------------------------------------------------------------------
1548 * If there exists no CommPkg for A, a CommPkg is generated using
1549 * equally load balanced partitionings
1550 *--------------------------------------------------------------------*/
1551 if (!comm_pkg)
1552 {
1553 hypre_MatvecCommPkgCreate(A);
1554 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1555 }
1556
1557 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1558 y_buf = mfem_hypre_CTAlloc_host(
1559 HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1560
1561 if (num_cols_offd)
1562 {
1563 // Disable the use of offdT for now, until we implement
1564 // hypre_CSRMatrixAbsMatvec on device.
1565#if MFEM_HYPRE_VERSION >= 21100 && 0
1566 if (A->offdT)
1567 {
1568 // offdT is optional. Used only if it's present.
1569 hypre_CSRMatrixAbsMatvec(A->offdT, alpha, x, 0., y_tmp);
1570 }
1571 else
1572#endif
1573 {
1574 hypre_CSRMatrixAbsMatvecT(offd, alpha, x, 0., y_tmp);
1575 }
1576 }
1577
1578 comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg, y_tmp, y_buf);
1579
1580 // Disable the use of diagT for now, until we implement
1581 // hypre_CSRMatrixAbsMatvec on device.
1582#if MFEM_HYPRE_VERSION >= 21100 && 0
1583 if (A->diagT)
1584 {
1585 // diagT is optional. Used only if it's present.
1586 hypre_CSRMatrixAbsMatvec(A->diagT, alpha, x, beta, y);
1587 }
1588 else
1589#endif
1590 {
1591 hypre_CSRMatrixAbsMatvecT(diag, alpha, x, beta, y);
1592 }
1593
1594 hypre_ParCSRCommHandleDestroy(comm_handle);
1595
1596 for (i = 0; i < num_sends; i++)
1597 {
1598 end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
1599 for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
1600 {
1601 jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
1602 y[jj] += y_buf[j];
1603 }
1604 }
1605
1606 mfem_hypre_TFree_host(y_buf);
1607 mfem_hypre_TFree_host(y_tmp);
1608}
1609
1610/* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1611void hypre_ParCSRMatrixBooleanMatvec(hypre_ParCSRMatrix *A,
1612 HYPRE_Bool alpha,
1613 HYPRE_Bool *x,
1614 HYPRE_Bool beta,
1615 HYPRE_Bool *y)
1616{
1617 hypre_ParCSRCommHandle *comm_handle;
1618 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1619 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1620 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1621
1622 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1623 HYPRE_Int num_sends, i, j, index;
1624
1625 HYPRE_Bool *x_tmp, *x_buf;
1626
1627 x_tmp = mfem_hypre_CTAlloc_host(HYPRE_Bool, num_cols_offd);
1628
1629 /*---------------------------------------------------------------------
1630 * If there exists no CommPkg for A, a CommPkg is generated using
1631 * equally load balanced partitionings
1632 *--------------------------------------------------------------------*/
1633 if (!comm_pkg)
1634 {
1635 hypre_MatvecCommPkgCreate(A);
1636 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1637 }
1638
1639 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1640 x_buf = mfem_hypre_CTAlloc_host(
1641 HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1642
1643 index = 0;
1644 for (i = 0; i < num_sends; i++)
1645 {
1646 j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1647 for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1648 {
1649 x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1650 }
1651 }
1652
1653 comm_handle = hypre_ParCSRCommHandleCreate_bool(1, comm_pkg, x_buf, x_tmp);
1654
1655 hypre_CSRMatrixBooleanMatvec(diag, alpha, x, beta, y);
1656
1657 hypre_ParCSRCommHandleDestroy(comm_handle);
1658
1659 if (num_cols_offd)
1660 {
1661 hypre_CSRMatrixBooleanMatvec(offd, alpha, x_tmp, 1, y);
1662 }
1663
1664 mfem_hypre_TFree_host(x_buf);
1665 mfem_hypre_TFree_host(x_tmp);
1666}
1667
1668/* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
1669void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
1670 HYPRE_Bool alpha,
1671 HYPRE_Bool *x,
1672 HYPRE_Bool beta,
1673 HYPRE_Bool *y)
1674{
1675 hypre_ParCSRCommHandle *comm_handle;
1676 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1677 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1678 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1679 HYPRE_Bool *y_tmp;
1680 HYPRE_Bool *y_buf;
1681
1682 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1683
1684 HYPRE_Int i, j, jj, end, num_sends;
1685
1686 y_tmp = mfem_hypre_TAlloc_host(HYPRE_Bool, num_cols_offd);
1687
1688 /*---------------------------------------------------------------------
1689 * If there exists no CommPkg for A, a CommPkg is generated using
1690 * equally load balanced partitionings
1691 *--------------------------------------------------------------------*/
1692 if (!comm_pkg)
1693 {
1694 hypre_MatvecCommPkgCreate(A);
1695 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1696 }
1697
1698 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1699 y_buf = mfem_hypre_CTAlloc_host(
1700 HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1701
1702 if (num_cols_offd)
1703 {
1704 // Disable the use of offdT for now, until we implement
1705 // hypre_CSRMatrixBooleanMatvec on device.
1706#if MFEM_HYPRE_VERSION >= 21100 && 0
1707 if (A->offdT)
1708 {
1709 // offdT is optional. Used only if it's present.
1710 hypre_CSRMatrixBooleanMatvec(A->offdT, alpha, x, 0, y_tmp);
1711 }
1712 else
1713#endif
1714 {
1715 hypre_CSRMatrixBooleanMatvecT(offd, alpha, x, 0, y_tmp);
1716 }
1717 }
1718
1719 comm_handle = hypre_ParCSRCommHandleCreate_bool(2, comm_pkg, y_tmp, y_buf);
1720
1721 // Disable the use of diagT for now, until we implement
1722 // hypre_CSRMatrixBooleanMatvec on device.
1723#if MFEM_HYPRE_VERSION >= 21100 && 0
1724 if (A->diagT)
1725 {
1726 // diagT is optional. Used only if it's present.
1727 hypre_CSRMatrixBooleanMatvec(A->diagT, alpha, x, beta, y);
1728 }
1729 else
1730#endif
1731 {
1732 hypre_CSRMatrixBooleanMatvecT(diag, alpha, x, beta, y);
1733 }
1734
1735 hypre_ParCSRCommHandleDestroy(comm_handle);
1736
1737 for (i = 0; i < num_sends; i++)
1738 {
1739 end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
1740 for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
1741 {
1742 jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
1743 y[jj] = y[jj] || y_buf[j];
1744 }
1745 }
1746
1747 mfem_hypre_TFree_host(y_buf);
1748 mfem_hypre_TFree_host(y_tmp);
1749}
1750
1751HYPRE_Int
1752hypre_CSRMatrixSum(hypre_CSRMatrix *A,
1753 HYPRE_Complex beta,
1754 hypre_CSRMatrix *B)
1755{
1756 HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1757 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1758 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1759 HYPRE_Int nrows_A = hypre_CSRMatrixNumRows(A);
1760 HYPRE_Int ncols_A = hypre_CSRMatrixNumCols(A);
1761 HYPRE_Complex *B_data = hypre_CSRMatrixData(B);
1762 HYPRE_Int *B_i = hypre_CSRMatrixI(B);
1763 HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
1764 HYPRE_Int nrows_B = hypre_CSRMatrixNumRows(B);
1765 HYPRE_Int ncols_B = hypre_CSRMatrixNumCols(B);
1766
1767 HYPRE_Int ia, j, pos;
1768 HYPRE_Int *marker;
1769
1770 if (nrows_A != nrows_B || ncols_A != ncols_B)
1771 {
1772 return -1; /* error: incompatible matrix dimensions */
1773 }
1774
1775 marker = mfem_hypre_CTAlloc_host(HYPRE_Int, ncols_A);
1776 for (ia = 0; ia < ncols_A; ia++)
1777 {
1778 marker[ia] = -1;
1779 }
1780
1781 for (ia = 0; ia < nrows_A; ia++)
1782 {
1783 for (j = A_i[ia]; j < A_i[ia+1]; j++)
1784 {
1785 marker[A_j[j]] = j;
1786 }
1787
1788 for (j = B_i[ia]; j < B_i[ia+1]; j++)
1789 {
1790 pos = marker[B_j[j]];
1791 if (pos < A_i[ia])
1792 {
1793 return -2; /* error: found an entry in B that is not present in A */
1794 }
1795 A_data[pos] += beta * B_data[j];
1796 }
1797 }
1798
1799 mfem_hypre_TFree_host(marker);
1800 return 0;
1801}
1802
1803hypre_ParCSRMatrix *
1804hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
1805 hypre_ParCSRMatrix *B)
1806{
1807 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
1808 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1809 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1810 HYPRE_BigInt *A_cmap = hypre_ParCSRMatrixColMapOffd(A);
1811 HYPRE_Int A_cmap_size = hypre_CSRMatrixNumCols(A_offd);
1812 hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1813 hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1814 HYPRE_BigInt *B_cmap = hypre_ParCSRMatrixColMapOffd(B);
1815 HYPRE_Int B_cmap_size = hypre_CSRMatrixNumCols(B_offd);
1816 hypre_ParCSRMatrix *C;
1817 hypre_CSRMatrix *C_diag;
1818 hypre_CSRMatrix *C_offd;
1819 HYPRE_BigInt *C_cmap;
1820 HYPRE_Int im;
1821 HYPRE_Int cmap_differ;
1822
1823 /* Check if A_cmap and B_cmap are the same. */
1824 cmap_differ = 0;
1825 if (A_cmap_size != B_cmap_size)
1826 {
1827 cmap_differ = 1; /* A and B have different cmap_size */
1828 }
1829 else
1830 {
1831 for (im = 0; im < A_cmap_size; im++)
1832 {
1833 if (A_cmap[im] != B_cmap[im])
1834 {
1835 cmap_differ = 1; /* A and B have different cmap arrays */
1836 break;
1837 }
1838 }
1839 }
1840
1841 if ( cmap_differ == 0 )
1842 {
1843 /* A and B have the same column mapping for their off-diagonal blocks so
1844 we can sum the diagonal and off-diagonal blocks separately and reduce
1845 temporary memory usage. */
1846
1847 /* Add diagonals, off-diagonals, copy cmap. */
1848 C_diag = hypre_CSRMatrixAdd(A_diag, B_diag);
1849 if (!C_diag)
1850 {
1851 return NULL; /* error: A_diag and B_diag have different dimensions */
1852 }
1853 C_offd = hypre_CSRMatrixAdd(A_offd, B_offd);
1854 if (!C_offd)
1855 {
1856 hypre_CSRMatrixDestroy(C_diag);
1857 return NULL; /* error: A_offd and B_offd have different dimensions */
1858 }
1859 /* copy A_cmap -> C_cmap */
1860 C_cmap = mfem_hypre_TAlloc_host(HYPRE_BigInt, A_cmap_size);
1861 for (im = 0; im < A_cmap_size; im++)
1862 {
1863 C_cmap[im] = A_cmap[im];
1864 }
1865
1866 C = hypre_ParCSRMatrixCreate(comm,
1867 hypre_ParCSRMatrixGlobalNumRows(A),
1868 hypre_ParCSRMatrixGlobalNumCols(A),
1869 hypre_ParCSRMatrixRowStarts(A),
1870 hypre_ParCSRMatrixColStarts(A),
1871 hypre_CSRMatrixNumCols(C_offd),
1872 hypre_CSRMatrixNumNonzeros(C_diag),
1873 hypre_CSRMatrixNumNonzeros(C_offd));
1874
1875 /* In C, destroy diag/offd (allocated by Create) and replace them with
1876 C_diag/C_offd. */
1877 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(C));
1878 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(C));
1879 hypre_ParCSRMatrixDiag(C) = C_diag;
1880 hypre_ParCSRMatrixOffd(C) = C_offd;
1881
1882 hypre_ParCSRMatrixColMapOffd(C) = C_cmap;
1883 }
1884 else
1885 {
1886 /* A and B have different column mappings for their off-diagonal blocks so
1887 we need to use the column maps to create full-width CSR matrices. */
1888
1889 int ierr = 0;
1890 hypre_CSRMatrix * csr_A;
1891 hypre_CSRMatrix * csr_B;
1892 hypre_CSRMatrix * csr_C_temp;
1893
1894 /* merge diag and off-diag portions of A */
1895 csr_A = hypre_MergeDiagAndOffd(A);
1896
1897 /* merge diag and off-diag portions of B */
1898 csr_B = hypre_MergeDiagAndOffd(B);
1899
1900 /* add A and B */
1901 csr_C_temp = hypre_CSRMatrixAdd(csr_A, csr_B);
1902
1903 /* delete CSR versions of A and B */
1904 ierr += hypre_CSRMatrixDestroy(csr_A);
1905 ierr += hypre_CSRMatrixDestroy(csr_B);
1906
1907 /* create a new empty ParCSR matrix to contain the sum */
1908 C = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
1909 hypre_ParCSRMatrixGlobalNumRows(A),
1910 hypre_ParCSRMatrixGlobalNumCols(A),
1911 hypre_ParCSRMatrixRowStarts(A),
1912 hypre_ParCSRMatrixColStarts(A),
1913 0, 0, 0);
1914
1915 /* split C into diag and off-diag portions */
1916 /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
1917 number of columns in csr_C_temp which is the global number of columns
1918 in A and B. This does not scale well. */
1919 ierr += GenerateDiagAndOffd(csr_C_temp, C,
1920 hypre_ParCSRMatrixFirstColDiag(A),
1921 hypre_ParCSRMatrixLastColDiag(A));
1922
1923 /* delete CSR version of C */
1924 ierr += hypre_CSRMatrixDestroy(csr_C_temp);
1925
1926 MFEM_VERIFY(ierr == 0, "");
1927 }
1928
1929 /* hypre_ParCSRMatrixSetNumNonzeros(A); */
1930
1931 /* Make sure that the first entry in each row is the diagonal one. */
1932 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(C));
1933
1934 /* C owns diag, offd, and cmap. */
1935 hypre_ParCSRMatrixSetDataOwner(C, 1);
1936
1937#if MFEM_HYPRE_VERSION <= 22200
1938 /* C does not own row and column starts. */
1939 hypre_ParCSRMatrixSetRowStartsOwner(C, 0);
1940 hypre_ParCSRMatrixSetColStartsOwner(C, 0);
1941#endif
1942
1943 return C;
1944}
1945
1946HYPRE_Int
1947hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
1948 HYPRE_Complex beta,
1949 hypre_ParCSRMatrix *B)
1950{
1951 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1952 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1953 hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1954 hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1955 HYPRE_Int ncols_B_offd = hypre_CSRMatrixNumCols(B_offd);
1956 HYPRE_Int error;
1957
1958 error = hypre_CSRMatrixSum(A_diag, beta, B_diag);
1959 if (ncols_B_offd > 0) /* treat B_offd as zero if it has no columns */
1960 {
1961 error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
1962 }
1963
1964 return error;
1965}
1966
1967HYPRE_Int
1968hypre_CSRMatrixSetConstantValues(hypre_CSRMatrix *A,
1969 HYPRE_Complex value)
1970{
1971 HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1972 HYPRE_Int A_nnz = hypre_CSRMatrixNumNonzeros(A);
1973 HYPRE_Int ia;
1974
1975 for (ia = 0; ia < A_nnz; ia++)
1976 {
1977 A_data[ia] = value;
1978 }
1979
1980 return 0;
1981}
1982
1983HYPRE_Int
1984hypre_ParCSRMatrixSetConstantValues(hypre_ParCSRMatrix *A,
1985 HYPRE_Complex value)
1986{
1987 internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixDiag(A), value);
1988 internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixOffd(A), value);
1989
1990 return 0;
1991}
1992
1993} // namespace mfem::internal
1994
1995} // namespace mfem
1996
1997#endif // MFEM_USE_MPI
beta
Vector beta
Definition convection-diffusion.cpp:82
alpha
const real_t alpha
Definition ex15.cpp:369
hypre_parcsr.hpp
HYPRE_BigInt
HYPRE_Int HYPRE_BigInt
Definition hypre_parcsr.hpp:27
index
int index(int i, int j, int nx, int ny)
Definition life.cpp:235
a
real_t a
Definition lissajous.cpp:41
doxygen 1.11.0