MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
hypre_parcsr.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 
23 namespace mfem
24 {
25 namespace 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 */
39 void 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 */
100 void 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 */
139 void 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 */
185 void 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 i = 0; i < num_rowscols_to_elim; i++)
306  {
307  irow = rowscols_to_elim[i];
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 */
325 void 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 */
394 void 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 */
452 void 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 
491 void 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  hypre_ParCSRMatrixSetRowStartsOwner(*Ae, 0);
512  hypre_ParCSRMatrixSetColStartsOwner(*Ae, 0);
513 
514  hypre_CSRMatrix *Ae_diag = hypre_ParCSRMatrixDiag(*Ae);
515  hypre_CSRMatrix *Ae_offd = hypre_ParCSRMatrixOffd(*Ae);
516  HYPRE_Int Ae_offd_ncols;
517 
518  HYPRE_Int num_offd_cols_to_elim;
519  HYPRE_Int *offd_cols_to_elim;
520 
521  HYPRE_BigInt *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
522  HYPRE_BigInt *Ae_col_map_offd;
523 
524  HYPRE_Int *col_mark;
525  HYPRE_Int *col_remap;
526 
527  /* figure out which offd cols should be eliminated */
528  {
529  hypre_ParCSRCommHandle *comm_handle;
530  hypre_ParCSRCommPkg *comm_pkg;
531  HYPRE_Int num_sends, *int_buf_data;
532  HYPRE_Int index, start;
533 
534  HYPRE_Int *eliminate_diag_col = mfem_hypre_CTAlloc_host(HYPRE_Int,
535  A_diag_ncols);
536  HYPRE_Int *eliminate_offd_col = mfem_hypre_CTAlloc_host(HYPRE_Int,
537  A_offd_ncols);
538 
539  /* make sure A has a communication package */
540  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
541  if (!comm_pkg)
542  {
543  hypre_MatvecCommPkgCreate(A);
544  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
545  }
546 
547  /* which of the local rows are to be eliminated */
548  for (i = 0; i < A_diag_ncols; i++)
549  {
550  eliminate_diag_col[i] = 0;
551  }
552  for (i = 0; i < num_rowscols_to_elim; i++)
553  {
554  eliminate_diag_col[rowscols_to_elim[i]] = 1;
555  }
556 
557  /* use a Matvec communication pattern to find (in eliminate_col)
558  which of the local offd columns are to be eliminated */
559  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
560  int_buf_data = mfem_hypre_CTAlloc_host(
561  HYPRE_Int,
562  hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
563  index = 0;
564  for (i = 0; i < num_sends; i++)
565  {
566  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
567  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
568  {
569  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
570  int_buf_data[index++] = eliminate_diag_col[k];
571  }
572  }
573  comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
574  int_buf_data,
575  eliminate_offd_col);
576 
577  /* eliminate diagonal part, overlapping it with communication */
578  if (ignore_rows)
579  {
580  hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
581  0, nullptr,
582  num_rowscols_to_elim, rowscols_to_elim,
583  NULL);
584 
585  hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
586  0, nullptr,
587  num_rowscols_to_elim, rowscols_to_elim,
588  1, NULL);
589  }
590  else
591  {
592  hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
593  num_rowscols_to_elim, rowscols_to_elim,
594  num_rowscols_to_elim, rowscols_to_elim,
595  NULL);
596 
597  hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
598  num_rowscols_to_elim, rowscols_to_elim,
599  num_rowscols_to_elim, rowscols_to_elim,
600  1, NULL);
601  }
602 
603  hypre_CSRMatrixReorder(Ae_diag);
604 
605  /* finish the communication */
606  hypre_ParCSRCommHandleDestroy(comm_handle);
607 
608  /* received eliminate_col[], count offd columns to eliminate */
609  num_offd_cols_to_elim = 0;
610  for (i = 0; i < A_offd_ncols; i++)
611  {
612  if (eliminate_offd_col[i]) { num_offd_cols_to_elim++; }
613  }
614 
615  offd_cols_to_elim = mfem_hypre_CTAlloc_host(HYPRE_Int,
616  num_offd_cols_to_elim);
617 
618  /* get a list of offd column indices and coefs */
619  num_offd_cols_to_elim = 0;
620  for (i = 0; i < A_offd_ncols; i++)
621  {
622  if (eliminate_offd_col[i])
623  {
624  offd_cols_to_elim[num_offd_cols_to_elim++] = i;
625  }
626  }
627 
628  mfem_hypre_TFree_host(int_buf_data);
629  mfem_hypre_TFree_host(eliminate_offd_col);
630  mfem_hypre_TFree_host(eliminate_diag_col);
631  }
632 
633  /* eliminate the off-diagonal part */
634  col_mark = mfem_hypre_CTAlloc_host(HYPRE_Int, A_offd_ncols);
635  col_remap = mfem_hypre_CTAlloc_host(HYPRE_Int, A_offd_ncols);
636 
637  if (ignore_rows)
638  {
639  hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
640  0, nullptr,
641  num_offd_cols_to_elim, offd_cols_to_elim,
642  col_mark);
643 
644  for (i = k = 0; i < A_offd_ncols; i++)
645  {
646  if (col_mark[i]) { col_remap[i] = k++; }
647  }
648 
649  hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
650  0, nullptr,
651  num_offd_cols_to_elim, offd_cols_to_elim,
652  0, col_remap);
653  }
654  else
655  {
656  hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
657  num_rowscols_to_elim, rowscols_to_elim,
658  num_offd_cols_to_elim, offd_cols_to_elim,
659  col_mark);
660 
661  for (i = k = 0; i < A_offd_ncols; i++)
662  {
663  if (col_mark[i]) { col_remap[i] = k++; }
664  }
665 
666  hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
667  num_rowscols_to_elim, rowscols_to_elim,
668  num_offd_cols_to_elim, offd_cols_to_elim,
669  0, col_remap);
670  }
671 
672  /* create col_map_offd for Ae */
673  Ae_offd_ncols = 0;
674  for (i = 0; i < A_offd_ncols; i++)
675  {
676  if (col_mark[i]) { Ae_offd_ncols++; }
677  }
678 
679  Ae_col_map_offd = mfem_hypre_CTAlloc_host(HYPRE_BigInt, Ae_offd_ncols);
680 
681  Ae_offd_ncols = 0;
682  for (i = 0; i < A_offd_ncols; i++)
683  {
684  if (col_mark[i])
685  {
686  Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
687  }
688  }
689 
690  hypre_ParCSRMatrixColMapOffd(*Ae) = Ae_col_map_offd;
691  hypre_CSRMatrixNumCols(Ae_offd) = Ae_offd_ncols;
692 
693  mfem_hypre_TFree_host(col_remap);
694  mfem_hypre_TFree_host(col_mark);
695  mfem_hypre_TFree_host(offd_cols_to_elim);
696 
697  hypre_ParCSRMatrixSetNumNonzeros(*Ae);
698  hypre_MatvecCommPkgCreate(*Ae);
699 }
700 
701 
702 // Eliminate rows from the diagonal and off-diagonal blocks of the matrix
703 void hypre_ParCSRMatrixEliminateRows(hypre_ParCSRMatrix *A,
704  HYPRE_Int num_rows_to_elim,
705  const HYPRE_Int *rows_to_elim)
706 {
707  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
708  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
709  hypre_CSRMatrixEliminateRows(A_diag, num_rows_to_elim, rows_to_elim);
710  hypre_CSRMatrixEliminateRows(A_offd, num_rows_to_elim, rows_to_elim);
711 }
712 
713 
714 /*--------------------------------------------------------------------------
715  * Split
716  *--------------------------------------------------------------------------*/
717 
718 void hypre_CSRMatrixSplit(hypre_CSRMatrix *A,
719  HYPRE_Int nr, HYPRE_Int nc,
720  HYPRE_Int *row_block_num, HYPRE_Int *col_block_num,
721  hypre_CSRMatrix **blocks)
722 {
723  HYPRE_Int i, j, k, bi, bj;
724 
725  HYPRE_Int* A_i = hypre_CSRMatrixI(A);
726  HYPRE_Int* A_j = hypre_CSRMatrixJ(A);
727  HYPRE_Complex* A_data = hypre_CSRMatrixData(A);
728 
729  HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
730  HYPRE_Int A_cols = hypre_CSRMatrixNumCols(A);
731 
732  HYPRE_Int *num_rows = mfem_hypre_CTAlloc_host(HYPRE_Int, nr);
733  HYPRE_Int *num_cols = mfem_hypre_CTAlloc_host(HYPRE_Int, nc);
734 
735  HYPRE_Int *block_row = mfem_hypre_TAlloc_host(HYPRE_Int, A_rows);
736  HYPRE_Int *block_col = mfem_hypre_TAlloc_host(HYPRE_Int, A_cols);
737 
738  for (i = 0; i < A_rows; i++)
739  {
740  block_row[i] = num_rows[row_block_num[i]]++;
741  }
742  for (j = 0; j < A_cols; j++)
743  {
744  block_col[j] = num_cols[col_block_num[j]]++;
745  }
746 
747  /* allocate the blocks */
748  for (i = 0; i < nr; i++)
749  {
750  for (j = 0; j < nc; j++)
751  {
752  hypre_CSRMatrix *B = hypre_CSRMatrixCreate(num_rows[i],
753  num_cols[j], 0);
754  hypre_CSRMatrixI(B) = mfem_hypre_CTAlloc_host(HYPRE_Int,
755  num_rows[i] + 1);
756 #if MFEM_HYPRE_VERSION >= 21800
757  hypre_CSRMatrixMemoryLocation(B) = HYPRE_MEMORY_HOST;
758 #endif
759  blocks[i*nc + j] = B;
760  }
761  }
762 
763  /* count block row nnz */
764  for (i = 0; i < A_rows; i++)
765  {
766  bi = row_block_num[i];
767  for (j = A_i[i]; j < A_i[i+1]; j++)
768  {
769  bj = col_block_num[A_j[j]];
770  hypre_CSRMatrix *B = blocks[bi*nc + bj];
771  hypre_CSRMatrixI(B)[block_row[i] + 1]++;
772  }
773  }
774 
775  /* count block nnz */
776  for (k = 0; k < nr*nc; k++)
777  {
778  hypre_CSRMatrix *B = blocks[k];
779  HYPRE_Int* B_i = hypre_CSRMatrixI(B);
780 
781  HYPRE_Int nnz = 0, rs;
782  for (int k = 1; k <= hypre_CSRMatrixNumRows(B); k++)
783  {
784  rs = B_i[k], B_i[k] = nnz, nnz += rs;
785  }
786 
787  hypre_CSRMatrixJ(B) = mfem_hypre_TAlloc_host(HYPRE_Int, nnz);
788  hypre_CSRMatrixData(B) = mfem_hypre_TAlloc_host(HYPRE_Complex, nnz);
789  hypre_CSRMatrixNumNonzeros(B) = nnz;
790  }
791 
792  /* populate blocks */
793  for (i = 0; i < A_rows; i++)
794  {
795  bi = row_block_num[i];
796  for (j = A_i[i]; j < A_i[i+1]; j++)
797  {
798  k = A_j[j];
799  bj = col_block_num[k];
800  hypre_CSRMatrix *B = blocks[bi*nc + bj];
801  HYPRE_Int *bii = hypre_CSRMatrixI(B) + block_row[i] + 1;
802  hypre_CSRMatrixJ(B)[*bii] = block_col[k];
803  hypre_CSRMatrixData(B)[*bii] = A_data[j];
804  (*bii)++;
805  }
806  }
807 
808  mfem_hypre_TFree_host(block_col);
809  mfem_hypre_TFree_host(block_row);
810 
811  mfem_hypre_TFree_host(num_cols);
812  mfem_hypre_TFree_host(num_rows);
813 }
814 
815 
816 void hypre_ParCSRMatrixSplit(hypre_ParCSRMatrix *A,
817  HYPRE_Int nr, HYPRE_Int nc,
818  hypre_ParCSRMatrix **blocks,
819  int interleaved_rows, int interleaved_cols)
820 {
821  HYPRE_Int i, j, k;
822 
823  MPI_Comm comm = hypre_ParCSRMatrixComm(A);
824 
825  hypre_CSRMatrix *Adiag = hypre_ParCSRMatrixDiag(A);
826  hypre_CSRMatrix *Aoffd = hypre_ParCSRMatrixOffd(A);
827 
828  HYPRE_BigInt global_rows = hypre_ParCSRMatrixGlobalNumRows(A);
829  HYPRE_BigInt global_cols = hypre_ParCSRMatrixGlobalNumCols(A);
830 
831  HYPRE_Int local_rows = hypre_CSRMatrixNumRows(Adiag);
832  HYPRE_Int local_cols = hypre_CSRMatrixNumCols(Adiag);
833  HYPRE_Int offd_cols = hypre_CSRMatrixNumCols(Aoffd);
834 
835  hypre_assert(local_rows % nr == 0 && local_cols % nc == 0);
836  hypre_assert(global_rows % nr == 0 && global_cols % nc == 0);
837 
838  HYPRE_Int block_rows = local_rows / nr;
839  HYPRE_Int block_cols = local_cols / nc;
840  HYPRE_Int num_blocks = nr * nc;
841 
842  /* mark local rows and columns with block number */
843  HYPRE_Int *row_block_num = mfem_hypre_TAlloc_host(HYPRE_Int, local_rows);
844  HYPRE_Int *col_block_num = mfem_hypre_TAlloc_host(HYPRE_Int, local_cols);
845 
846  for (i = 0; i < local_rows; i++)
847  {
848  row_block_num[i] = interleaved_rows ? (i % nr) : (i / block_rows);
849  }
850  for (i = 0; i < local_cols; i++)
851  {
852  col_block_num[i] = interleaved_cols ? (i % nc) : (i / block_cols);
853  }
854 
855  /* determine the block numbers for offd columns */
856  HYPRE_BigInt *offd_col_block_num = mfem_hypre_TAlloc_host(HYPRE_BigInt,
857  offd_cols);
858  hypre_ParCSRCommHandle *comm_handle;
859  HYPRE_BigInt *int_buf_data;
860  {
861  /* make sure A has a communication package */
862  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
863  if (!comm_pkg)
864  {
865  hypre_MatvecCommPkgCreate(A);
866  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
867  }
868 
869  /* calculate the final global column numbers for each block */
870  HYPRE_Int *count = mfem_hypre_CTAlloc_host(HYPRE_Int, nc);
871  HYPRE_BigInt *block_global_col = mfem_hypre_TAlloc_host(HYPRE_BigInt,
872  local_cols);
873  HYPRE_BigInt first_col = hypre_ParCSRMatrixFirstColDiag(A) / nc;
874  for (i = 0; i < local_cols; i++)
875  {
876  block_global_col[i] = first_col + count[col_block_num[i]]++;
877  }
878  mfem_hypre_TFree_host(count);
879 
880  /* use a Matvec communication pattern to determine offd_col_block_num */
881  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
882  int_buf_data = mfem_hypre_CTAlloc_host(
883  HYPRE_BigInt,
884  hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
885  HYPRE_Int start, index = 0;
886  for (i = 0; i < num_sends; i++)
887  {
888  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
889  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
890  {
891  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
892  int_buf_data[index++] = col_block_num[k] + nc*block_global_col[k];
893  }
894  }
895 
896  mfem_hypre_TFree_host(block_global_col);
897 
898 #if MFEM_HYPRE_VERSION < 21600
899  const int job = 11;
900 #else
901  const int job = 21;
902 #endif
903  comm_handle = hypre_ParCSRCommHandleCreate(job, comm_pkg, int_buf_data,
904  offd_col_block_num);
905  }
906 
907  /* create the block matrices */
908  HYPRE_Int num_procs = 1;
909  if (!hypre_ParCSRMatrixAssumedPartition(A))
910  {
911  hypre_MPI_Comm_size(comm, &num_procs);
912  }
913 
914  HYPRE_BigInt *row_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, num_procs+1);
915  HYPRE_BigInt *col_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, num_procs+1);
916  for (i = 0; i <= num_procs; i++)
917  {
918  row_starts[i] = hypre_ParCSRMatrixRowStarts(A)[i] / nr;
919  col_starts[i] = hypre_ParCSRMatrixColStarts(A)[i] / nc;
920  }
921 
922  for (i = 0; i < num_blocks; i++)
923  {
924  blocks[i] = hypre_ParCSRMatrixCreate(comm,
925  global_rows/nr, global_cols/nc,
926  row_starts, col_starts, 0, 0, 0);
927  }
928 
929  /* split diag part */
930  hypre_CSRMatrix **csr_blocks = mfem_hypre_TAlloc_host(hypre_CSRMatrix*,
931  nr*nc);
932  hypre_CSRMatrixSplit(Adiag, nr, nc, row_block_num, col_block_num,
933  csr_blocks);
934 
935  for (i = 0; i < num_blocks; i++)
936  {
937  mfem_hypre_TFree_host(hypre_ParCSRMatrixDiag(blocks[i]));
938  hypre_ParCSRMatrixDiag(blocks[i]) = csr_blocks[i];
939  }
940 
941  /* finish communication, receive offd_col_block_num */
942  hypre_ParCSRCommHandleDestroy(comm_handle);
943  mfem_hypre_TFree_host(int_buf_data);
944 
945  /* decode global offd column numbers */
946  HYPRE_Int *offd_col_block_num_nc = mfem_hypre_TAlloc_host(HYPRE_Int,
947  offd_cols);
948  HYPRE_BigInt* offd_global_col = mfem_hypre_TAlloc_host(HYPRE_BigInt,
949  offd_cols);
950  for (i = 0; i < offd_cols; i++)
951  {
952  offd_global_col[i] = offd_col_block_num[i] / nc;
953  offd_col_block_num_nc[i] = offd_col_block_num[i] % nc;
954  }
955 
956  mfem_hypre_TFree_host(offd_col_block_num);
957 
958  /* split offd part */
959  hypre_CSRMatrixSplit(Aoffd, nr, nc, row_block_num, offd_col_block_num_nc,
960  csr_blocks);
961 
962  for (i = 0; i < num_blocks; i++)
963  {
964  mfem_hypre_TFree_host(hypre_ParCSRMatrixOffd(blocks[i]));
965  hypre_ParCSRMatrixOffd(blocks[i]) = csr_blocks[i];
966  }
967 
968  mfem_hypre_TFree_host(csr_blocks);
969  mfem_hypre_TFree_host(col_block_num);
970  mfem_hypre_TFree_host(row_block_num);
971 
972  /* update block col-maps */
973  for (int bi = 0; bi < nr; bi++)
974  {
975  for (int bj = 0; bj < nc; bj++)
976  {
977  hypre_ParCSRMatrix *block = blocks[bi*nc + bj];
978  hypre_CSRMatrix *block_offd = hypre_ParCSRMatrixOffd(block);
979  HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols(block_offd);
980 
981  HYPRE_BigInt *block_col_map = mfem_hypre_TAlloc_host(HYPRE_BigInt,
982  block_offd_cols);
983  for (i = j = 0; i < offd_cols; i++)
984  {
985  HYPRE_Int bn = offd_col_block_num_nc[i];
986  if (bn == bj) { block_col_map[j++] = offd_global_col[i]; }
987  }
988  hypre_assert(j == block_offd_cols);
989 
990  hypre_ParCSRMatrixColMapOffd(block) = block_col_map;
991  }
992  }
993 
994  mfem_hypre_TFree_host(offd_global_col);
995  mfem_hypre_TFree_host(offd_col_block_num_nc);
996 
997  /* finish the new matrices, make them own all the stuff */
998  for (i = 0; i < num_blocks; i++)
999  {
1000  hypre_ParCSRMatrixSetNumNonzeros(blocks[i]);
1001  hypre_MatvecCommPkgCreate(blocks[i]);
1002 
1003  hypre_ParCSRMatrixOwnsData(blocks[i]) = 1;
1004 
1005  /* only the first block will own the row/col_starts */
1006  hypre_ParCSRMatrixOwnsRowStarts(blocks[i]) = !i;
1007  hypre_ParCSRMatrixOwnsColStarts(blocks[i]) = !i;
1008  }
1009 }
1010 
1011 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
1012 void hypre_CSRMatrixAbsMatvec(hypre_CSRMatrix *A,
1013  HYPRE_Real alpha,
1014  HYPRE_Real *x,
1015  HYPRE_Real beta,
1016  HYPRE_Real *y)
1017 {
1018  HYPRE_Real *A_data = hypre_CSRMatrixData(A);
1019  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1020  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1021  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1022 
1023  HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
1024  HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
1025 
1026  HYPRE_Real *x_data = x;
1027  HYPRE_Real *y_data = y;
1028 
1029  HYPRE_Real temp, tempx;
1030 
1031  HYPRE_Int i, jj;
1032 
1033  HYPRE_Int m;
1034 
1035  HYPRE_Real xpar=0.7;
1036 
1037  /*-----------------------------------------------------------------------
1038  * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
1039  *-----------------------------------------------------------------------*/
1040 
1041  if (alpha == 0.0)
1042  {
1043  for (i = 0; i < num_rows; i++)
1044  {
1045  y_data[i] *= beta;
1046  }
1047  return;
1048  }
1049 
1050  /*-----------------------------------------------------------------------
1051  * y = (beta/alpha)*y
1052  *-----------------------------------------------------------------------*/
1053 
1054  temp = beta / alpha;
1055 
1056  if (temp != 1.0)
1057  {
1058  if (temp == 0.0)
1059  {
1060  for (i = 0; i < num_rows; i++)
1061  {
1062  y_data[i] = 0.0;
1063  }
1064  }
1065  else
1066  {
1067  for (i = 0; i < num_rows; i++)
1068  {
1069  y_data[i] *= temp;
1070  }
1071  }
1072  }
1073 
1074  /*-----------------------------------------------------------------
1075  * y += abs(A)*x
1076  *-----------------------------------------------------------------*/
1077 
1078  /* use rownnz pointer to do the abs(A)*x multiplication
1079  when num_rownnz is smaller than num_rows */
1080 
1081  if (num_rownnz < xpar*(num_rows))
1082  {
1083  for (i = 0; i < num_rownnz; i++)
1084  {
1085  m = A_rownnz[i];
1086 
1087  tempx = 0;
1088  for (jj = A_i[m]; jj < A_i[m+1]; jj++)
1089  {
1090  tempx += std::abs(A_data[jj])*x_data[A_j[jj]];
1091  }
1092  y_data[m] += tempx;
1093  }
1094  }
1095  else
1096  {
1097  for (i = 0; i < num_rows; i++)
1098  {
1099  tempx = 0;
1100  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1101  {
1102  tempx += std::abs(A_data[jj])*x_data[A_j[jj]];
1103  }
1104  y_data[i] += tempx;
1105  }
1106  }
1107 
1108  /*-----------------------------------------------------------------
1109  * y = alpha*y
1110  *-----------------------------------------------------------------*/
1111 
1112  if (alpha != 1.0)
1113  {
1114  for (i = 0; i < num_rows; i++)
1115  {
1116  y_data[i] *= alpha;
1117  }
1118  }
1119 
1120 }
1121 
1122 /* Based on hypre_CSRMatrixMatvecT in hypre's csr_matvec.c */
1123 void hypre_CSRMatrixAbsMatvecT(hypre_CSRMatrix *A,
1124  HYPRE_Real alpha,
1125  HYPRE_Real *x,
1126  HYPRE_Real beta,
1127  HYPRE_Real *y)
1128 {
1129  HYPRE_Real *A_data = hypre_CSRMatrixData(A);
1130  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1131  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1132  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1133  HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
1134 
1135  HYPRE_Real *x_data = x;
1136  HYPRE_Real *y_data = y;
1137 
1138  HYPRE_Int i, j, jj;
1139 
1140  HYPRE_Real temp;
1141 
1142  if (alpha == 0.0)
1143  {
1144  for (i = 0; i < num_cols; i++)
1145  {
1146  y_data[i] *= beta;
1147  }
1148  return;
1149  }
1150 
1151  /*-----------------------------------------------------------------------
1152  * y = (beta/alpha)*y
1153  *-----------------------------------------------------------------------*/
1154 
1155  temp = beta / alpha;
1156 
1157  if (temp != 1.0)
1158  {
1159  if (temp == 0.0)
1160  {
1161  for (i = 0; i < num_cols; i++)
1162  {
1163  y_data[i] = 0.0;
1164  }
1165  }
1166  else
1167  {
1168  for (i = 0; i < num_cols; i++)
1169  {
1170  y_data[i] *= temp;
1171  }
1172  }
1173  }
1174 
1175  /*-----------------------------------------------------------------
1176  * y += abs(A)^T*x
1177  *-----------------------------------------------------------------*/
1178 
1179  for (i = 0; i < num_rows; i++)
1180  {
1181  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1182  {
1183  j = A_j[jj];
1184  y_data[j] += std::abs(A_data[jj]) * x_data[i];
1185  }
1186  }
1187 
1188  /*-----------------------------------------------------------------
1189  * y = alpha*y
1190  *-----------------------------------------------------------------*/
1191 
1192  if (alpha != 1.0)
1193  {
1194  for (i = 0; i < num_cols; i++)
1195  {
1196  y_data[i] *= alpha;
1197  }
1198  }
1199 }
1200 
1201 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
1202 void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,
1203  HYPRE_Bool alpha,
1204  HYPRE_Bool *x,
1205  HYPRE_Bool beta,
1206  HYPRE_Bool *y)
1207 {
1208  /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1209  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1210  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1211  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1212 
1213  HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
1214  HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
1215 
1216  HYPRE_Bool *x_data = x;
1217  HYPRE_Bool *y_data = y;
1218 
1219  HYPRE_Bool temp, tempx;
1220 
1221  HYPRE_Int i, jj;
1222 
1223  HYPRE_Int m;
1224 
1225  HYPRE_Real xpar=0.7;
1226 
1227  /*-----------------------------------------------------------------------
1228  * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
1229  *-----------------------------------------------------------------------*/
1230 
1231  if (alpha == 0)
1232  {
1233 #ifdef HYPRE_USING_OPENMP
1234  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1235 #endif
1236  for (i = 0; i < num_rows; i++)
1237  {
1238  y_data[i] = y_data[i] && beta;
1239  }
1240  return;
1241  }
1242 
1243  /*-----------------------------------------------------------------------
1244  * y = (beta/alpha)*y
1245  *-----------------------------------------------------------------------*/
1246 
1247  if (beta == 0)
1248  {
1249 #ifdef HYPRE_USING_OPENMP
1250  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1251 #endif
1252  for (i = 0; i < num_rows; i++)
1253  {
1254  y_data[i] = 0;
1255  }
1256  }
1257  else
1258  {
1259  /* beta is true -> no change to y_data */
1260  }
1261 
1262  /*-----------------------------------------------------------------
1263  * y += A*x
1264  *-----------------------------------------------------------------*/
1265 
1266  /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows */
1267 
1268  if (num_rownnz < xpar*(num_rows))
1269  {
1270 #ifdef HYPRE_USING_OPENMP
1271  #pragma omp parallel for private(i,jj,m,tempx) HYPRE_SMP_SCHEDULE
1272 #endif
1273  for (i = 0; i < num_rownnz; i++)
1274  {
1275  m = A_rownnz[i];
1276 
1277  tempx = 0;
1278  for (jj = A_i[m]; jj < A_i[m+1]; jj++)
1279  {
1280  /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1281  tempx = tempx || x_data[A_j[jj]];
1282  }
1283  y_data[m] = y_data[m] || tempx;
1284  }
1285  }
1286  else
1287  {
1288 #ifdef HYPRE_USING_OPENMP
1289  #pragma omp parallel for private(i,jj,temp) HYPRE_SMP_SCHEDULE
1290 #endif
1291  for (i = 0; i < num_rows; i++)
1292  {
1293  temp = 0;
1294  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1295  {
1296  /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1297  temp = temp || x_data[A_j[jj]];
1298  }
1299  y_data[i] = y_data[i] || temp;
1300  }
1301  }
1302 
1303  /*-----------------------------------------------------------------
1304  * y = alpha*y
1305  *-----------------------------------------------------------------*/
1306  /* alpha is true */
1307 }
1308 
1309 /* Based on hypre_CSRMatrixMatvecT in hypre's csr_matvec.c */
1310 void hypre_CSRMatrixBooleanMatvecT(hypre_CSRMatrix *A,
1311  HYPRE_Bool alpha,
1312  HYPRE_Bool *x,
1313  HYPRE_Bool beta,
1314  HYPRE_Bool *y)
1315 {
1316  /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1317  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1318  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1319  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1320  HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
1321 
1322  HYPRE_Bool *x_data = x;
1323  HYPRE_Bool *y_data = y;
1324 
1325  HYPRE_Int i, j, jj;
1326 
1327  /*-----------------------------------------------------------------------
1328  * y = beta*y
1329  *-----------------------------------------------------------------------*/
1330 
1331  if (beta == 0)
1332  {
1333  for (i = 0; i < num_cols; i++)
1334  {
1335  y_data[i] = 0;
1336  }
1337  }
1338  else
1339  {
1340  /* beta is true -> no change to y_data */
1341  }
1342 
1343  /*-----------------------------------------------------------------------
1344  * Check if (alpha == 0)
1345  *-----------------------------------------------------------------------*/
1346 
1347  if (alpha == 0)
1348  {
1349  return;
1350  }
1351 
1352  /* alpha is true */
1353 
1354  /*-----------------------------------------------------------------
1355  * y += A^T*x
1356  *-----------------------------------------------------------------*/
1357  for (i = 0; i < num_rows; i++)
1358  {
1359  if (x_data[i] != 0)
1360  {
1361  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1362  {
1363  j = A_j[jj];
1364  /* y_data[j] += A_data[jj] * x_data[i]; */
1365  y_data[j] = 1;
1366  }
1367  }
1368  }
1369 }
1370 
1371 /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
1372  input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
1373 hypre_ParCSRCommHandle *
1374 hypre_ParCSRCommHandleCreate_bool(HYPRE_Int job,
1375  hypre_ParCSRCommPkg *comm_pkg,
1376  HYPRE_Bool *send_data,
1377  HYPRE_Bool *recv_data)
1378 {
1379  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1380  HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
1381  MPI_Comm comm = hypre_ParCSRCommPkgComm(comm_pkg);
1382 
1383  hypre_ParCSRCommHandle *comm_handle;
1384  HYPRE_Int num_requests;
1385  hypre_MPI_Request *requests;
1386 
1387  HYPRE_Int i, j;
1388  HYPRE_Int my_id, num_procs;
1389  HYPRE_Int ip, vec_start, vec_len;
1390 
1391  num_requests = num_sends + num_recvs;
1392  requests = mfem_hypre_CTAlloc_host(hypre_MPI_Request, num_requests);
1393 
1394  hypre_MPI_Comm_size(comm, &num_procs);
1395  hypre_MPI_Comm_rank(comm, &my_id);
1396 
1397  j = 0;
1398  switch (job)
1399  {
1400  case 1:
1401  {
1402  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1403  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1404  for (i = 0; i < num_recvs; i++)
1405  {
1406  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1407  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1408  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1409  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1410  ip, 0, comm, &requests[j++]);
1411  }
1412  for (i = 0; i < num_sends; i++)
1413  {
1414  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1415  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1416  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1417  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1418  ip, 0, comm, &requests[j++]);
1419  }
1420  break;
1421  }
1422  case 2:
1423  {
1424  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1425  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1426  for (i = 0; i < num_sends; i++)
1427  {
1428  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1429  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1430  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1431  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1432  ip, 0, comm, &requests[j++]);
1433  }
1434  for (i = 0; i < num_recvs; i++)
1435  {
1436  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1437  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1438  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1439  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1440  ip, 0, comm, &requests[j++]);
1441  }
1442  break;
1443  }
1444  }
1445  /*--------------------------------------------------------------------
1446  * set up comm_handle and return
1447  *--------------------------------------------------------------------*/
1448 
1449  comm_handle = mfem_hypre_CTAlloc_host(hypre_ParCSRCommHandle, 1);
1450 
1451  hypre_ParCSRCommHandleCommPkg(comm_handle) = comm_pkg;
1452  hypre_ParCSRCommHandleSendData(comm_handle) = send_data;
1453  hypre_ParCSRCommHandleRecvData(comm_handle) = recv_data;
1454  hypre_ParCSRCommHandleNumRequests(comm_handle) = num_requests;
1455  hypre_ParCSRCommHandleRequests(comm_handle) = requests;
1456 
1457  return comm_handle;
1458 }
1459 
1460 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1461 void hypre_ParCSRMatrixAbsMatvec(hypre_ParCSRMatrix *A,
1462  HYPRE_Real alpha,
1463  HYPRE_Real *x,
1464  HYPRE_Real beta,
1465  HYPRE_Real *y)
1466 {
1467  hypre_ParCSRCommHandle *comm_handle;
1468  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1469  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1470  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1471 
1472  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1473  HYPRE_Int num_sends, i, j, index;
1474 
1475  HYPRE_Real *x_tmp, *x_buf;
1476 
1477  x_tmp = mfem_hypre_CTAlloc_host(HYPRE_Real, num_cols_offd);
1478 
1479  /*---------------------------------------------------------------------
1480  * If there exists no CommPkg for A, a CommPkg is generated using
1481  * equally load balanced partitionings
1482  *--------------------------------------------------------------------*/
1483  if (!comm_pkg)
1484  {
1485  hypre_MatvecCommPkgCreate(A);
1486  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1487  }
1488 
1489  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1490  x_buf = mfem_hypre_CTAlloc_host(
1491  HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1492 
1493  index = 0;
1494  for (i = 0; i < num_sends; i++)
1495  {
1496  j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1497  for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1498  {
1499  x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1500  }
1501  }
1502 
1503  comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg, x_buf, x_tmp);
1504 
1505  hypre_CSRMatrixAbsMatvec(diag, alpha, x, beta, y);
1506 
1507  hypre_ParCSRCommHandleDestroy(comm_handle);
1508 
1509  if (num_cols_offd)
1510  {
1511  hypre_CSRMatrixAbsMatvec(offd, alpha, x_tmp, 1.0, y);
1512  }
1513 
1514  mfem_hypre_TFree_host(x_buf);
1515  mfem_hypre_TFree_host(x_tmp);
1516 }
1517 
1518 /* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
1519 void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
1520  HYPRE_Real alpha,
1521  HYPRE_Real *x,
1522  HYPRE_Real beta,
1523  HYPRE_Real *y)
1524 {
1525  hypre_ParCSRCommHandle *comm_handle;
1526  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1527  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1528  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1529  HYPRE_Real *y_tmp;
1530  HYPRE_Real *y_buf;
1531 
1532  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1533 
1534  HYPRE_Int i, j, jj, end, num_sends;
1535 
1536  y_tmp = mfem_hypre_TAlloc_host(HYPRE_Real, num_cols_offd);
1537 
1538  /*---------------------------------------------------------------------
1539  * If there exists no CommPkg for A, a CommPkg is generated using
1540  * equally load balanced partitionings
1541  *--------------------------------------------------------------------*/
1542  if (!comm_pkg)
1543  {
1544  hypre_MatvecCommPkgCreate(A);
1545  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1546  }
1547 
1548  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1549  y_buf = mfem_hypre_CTAlloc_host(
1550  HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1551 
1552  if (num_cols_offd)
1553  {
1554 #if MFEM_HYPRE_VERSION >= 21100
1555  if (A->offdT)
1556  {
1557  // offdT is optional. Used only if it's present.
1558  hypre_CSRMatrixAbsMatvec(A->offdT, alpha, x, 0., y_tmp);
1559  }
1560  else
1561 #endif
1562  {
1563  hypre_CSRMatrixAbsMatvecT(offd, alpha, x, 0., y_tmp);
1564  }
1565  }
1566 
1567  comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg, y_tmp, y_buf);
1568 
1569 #if MFEM_HYPRE_VERSION >= 21100
1570  if (A->diagT)
1571  {
1572  // diagT is optional. Used only if it's present.
1573  hypre_CSRMatrixAbsMatvec(A->diagT, alpha, x, beta, y);
1574  }
1575  else
1576 #endif
1577  {
1578  hypre_CSRMatrixAbsMatvecT(diag, alpha, x, beta, y);
1579  }
1580 
1581  hypre_ParCSRCommHandleDestroy(comm_handle);
1582 
1583  for (i = 0; i < num_sends; i++)
1584  {
1585  end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
1586  for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
1587  {
1588  jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
1589  y[jj] += y_buf[j];
1590  }
1591  }
1592 
1593  mfem_hypre_TFree_host(y_buf);
1594  mfem_hypre_TFree_host(y_tmp);
1595 }
1596 
1597 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1598 void hypre_ParCSRMatrixBooleanMatvec(hypre_ParCSRMatrix *A,
1599  HYPRE_Bool alpha,
1600  HYPRE_Bool *x,
1601  HYPRE_Bool beta,
1602  HYPRE_Bool *y)
1603 {
1604  hypre_ParCSRCommHandle *comm_handle;
1605  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1606  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1607  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1608 
1609  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1610  HYPRE_Int num_sends, i, j, index;
1611 
1612  HYPRE_Bool *x_tmp, *x_buf;
1613 
1614  x_tmp = mfem_hypre_CTAlloc_host(HYPRE_Bool, num_cols_offd);
1615 
1616  /*---------------------------------------------------------------------
1617  * If there exists no CommPkg for A, a CommPkg is generated using
1618  * equally load balanced partitionings
1619  *--------------------------------------------------------------------*/
1620  if (!comm_pkg)
1621  {
1622  hypre_MatvecCommPkgCreate(A);
1623  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1624  }
1625 
1626  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1627  x_buf = mfem_hypre_CTAlloc_host(
1628  HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1629 
1630  index = 0;
1631  for (i = 0; i < num_sends; i++)
1632  {
1633  j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1634  for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1635  {
1636  x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1637  }
1638  }
1639 
1640  comm_handle = hypre_ParCSRCommHandleCreate_bool(1, comm_pkg, x_buf, x_tmp);
1641 
1642  hypre_CSRMatrixBooleanMatvec(diag, alpha, x, beta, y);
1643 
1644  hypre_ParCSRCommHandleDestroy(comm_handle);
1645 
1646  if (num_cols_offd)
1647  {
1648  hypre_CSRMatrixBooleanMatvec(offd, alpha, x_tmp, 1, y);
1649  }
1650 
1651  mfem_hypre_TFree_host(x_buf);
1652  mfem_hypre_TFree_host(x_tmp);
1653 }
1654 
1655 /* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
1656 void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
1657  HYPRE_Bool alpha,
1658  HYPRE_Bool *x,
1659  HYPRE_Bool beta,
1660  HYPRE_Bool *y)
1661 {
1662  hypre_ParCSRCommHandle *comm_handle;
1663  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1664  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1665  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1666  HYPRE_Bool *y_tmp;
1667  HYPRE_Bool *y_buf;
1668 
1669  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1670 
1671  HYPRE_Int i, j, jj, end, num_sends;
1672 
1673  y_tmp = mfem_hypre_TAlloc_host(HYPRE_Bool, num_cols_offd);
1674 
1675  /*---------------------------------------------------------------------
1676  * If there exists no CommPkg for A, a CommPkg is generated using
1677  * equally load balanced partitionings
1678  *--------------------------------------------------------------------*/
1679  if (!comm_pkg)
1680  {
1681  hypre_MatvecCommPkgCreate(A);
1682  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1683  }
1684 
1685  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1686  y_buf = mfem_hypre_CTAlloc_host(
1687  HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1688 
1689  if (num_cols_offd)
1690  {
1691 #if MFEM_HYPRE_VERSION >= 21100
1692  if (A->offdT)
1693  {
1694  // offdT is optional. Used only if it's present.
1695  hypre_CSRMatrixBooleanMatvec(A->offdT, alpha, x, 0, y_tmp);
1696  }
1697  else
1698 #endif
1699  {
1700  hypre_CSRMatrixBooleanMatvecT(offd, alpha, x, 0, y_tmp);
1701  }
1702  }
1703 
1704  comm_handle = hypre_ParCSRCommHandleCreate_bool(2, comm_pkg, y_tmp, y_buf);
1705 
1706 #if MFEM_HYPRE_VERSION >= 21100
1707  if (A->diagT)
1708  {
1709  // diagT is optional. Used only if it's present.
1710  hypre_CSRMatrixBooleanMatvec(A->diagT, alpha, x, beta, y);
1711  }
1712  else
1713 #endif
1714  {
1715  hypre_CSRMatrixBooleanMatvecT(diag, alpha, x, beta, y);
1716  }
1717 
1718  hypre_ParCSRCommHandleDestroy(comm_handle);
1719 
1720  for (i = 0; i < num_sends; i++)
1721  {
1722  end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
1723  for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
1724  {
1725  jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
1726  y[jj] = y[jj] || y_buf[j];
1727  }
1728  }
1729 
1730  mfem_hypre_TFree_host(y_buf);
1731  mfem_hypre_TFree_host(y_tmp);
1732 }
1733 
1734 HYPRE_Int
1735 hypre_CSRMatrixSum(hypre_CSRMatrix *A,
1736  HYPRE_Complex beta,
1737  hypre_CSRMatrix *B)
1738 {
1739  HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1740  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1741  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1742  HYPRE_Int nrows_A = hypre_CSRMatrixNumRows(A);
1743  HYPRE_Int ncols_A = hypre_CSRMatrixNumCols(A);
1744  HYPRE_Complex *B_data = hypre_CSRMatrixData(B);
1745  HYPRE_Int *B_i = hypre_CSRMatrixI(B);
1746  HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
1747  HYPRE_Int nrows_B = hypre_CSRMatrixNumRows(B);
1748  HYPRE_Int ncols_B = hypre_CSRMatrixNumCols(B);
1749 
1750  HYPRE_Int ia, j, pos;
1751  HYPRE_Int *marker;
1752 
1753  if (nrows_A != nrows_B || ncols_A != ncols_B)
1754  {
1755  return -1; /* error: incompatible matrix dimensions */
1756  }
1757 
1758  marker = mfem_hypre_CTAlloc_host(HYPRE_Int, ncols_A);
1759  for (ia = 0; ia < ncols_A; ia++)
1760  {
1761  marker[ia] = -1;
1762  }
1763 
1764  for (ia = 0; ia < nrows_A; ia++)
1765  {
1766  for (j = A_i[ia]; j < A_i[ia+1]; j++)
1767  {
1768  marker[A_j[j]] = j;
1769  }
1770 
1771  for (j = B_i[ia]; j < B_i[ia+1]; j++)
1772  {
1773  pos = marker[B_j[j]];
1774  if (pos < A_i[ia])
1775  {
1776  return -2; /* error: found an entry in B that is not present in A */
1777  }
1778  A_data[pos] += beta * B_data[j];
1779  }
1780  }
1781 
1782  mfem_hypre_TFree_host(marker);
1783  return 0;
1784 }
1785 
1786 hypre_ParCSRMatrix *
1787 hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
1788  hypre_ParCSRMatrix *B)
1789 {
1790  MPI_Comm comm = hypre_ParCSRMatrixComm(A);
1791  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1792  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1793  HYPRE_BigInt *A_cmap = hypre_ParCSRMatrixColMapOffd(A);
1794  HYPRE_Int A_cmap_size = hypre_CSRMatrixNumCols(A_offd);
1795  hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1796  hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1797  HYPRE_BigInt *B_cmap = hypre_ParCSRMatrixColMapOffd(B);
1798  HYPRE_Int B_cmap_size = hypre_CSRMatrixNumCols(B_offd);
1799  hypre_ParCSRMatrix *C;
1800  hypre_CSRMatrix *C_diag;
1801  hypre_CSRMatrix *C_offd;
1802  HYPRE_BigInt *C_cmap;
1803  HYPRE_Int im;
1804  HYPRE_Int cmap_differ;
1805 
1806  /* Check if A_cmap and B_cmap are the same. */
1807  cmap_differ = 0;
1808  if (A_cmap_size != B_cmap_size)
1809  {
1810  cmap_differ = 1; /* A and B have different cmap_size */
1811  }
1812  else
1813  {
1814  for (im = 0; im < A_cmap_size; im++)
1815  {
1816  if (A_cmap[im] != B_cmap[im])
1817  {
1818  cmap_differ = 1; /* A and B have different cmap arrays */
1819  break;
1820  }
1821  }
1822  }
1823 
1824  if ( cmap_differ == 0 )
1825  {
1826  /* A and B have the same column mapping for their off-diagonal blocks so
1827  we can sum the diagonal and off-diagonal blocks separately and reduce
1828  temporary memory usage. */
1829 
1830  /* Add diagonals, off-diagonals, copy cmap. */
1831  C_diag = hypre_CSRMatrixAdd(A_diag, B_diag);
1832  if (!C_diag)
1833  {
1834  return NULL; /* error: A_diag and B_diag have different dimensions */
1835  }
1836  C_offd = hypre_CSRMatrixAdd(A_offd, B_offd);
1837  if (!C_offd)
1838  {
1839  hypre_CSRMatrixDestroy(C_diag);
1840  return NULL; /* error: A_offd and B_offd have different dimensions */
1841  }
1842  /* copy A_cmap -> C_cmap */
1843  C_cmap = mfem_hypre_TAlloc_host(HYPRE_BigInt, A_cmap_size);
1844  for (im = 0; im < A_cmap_size; im++)
1845  {
1846  C_cmap[im] = A_cmap[im];
1847  }
1848 
1849  C = hypre_ParCSRMatrixCreate(comm,
1850  hypre_ParCSRMatrixGlobalNumRows(A),
1851  hypre_ParCSRMatrixGlobalNumCols(A),
1852  hypre_ParCSRMatrixRowStarts(A),
1853  hypre_ParCSRMatrixColStarts(A),
1854  hypre_CSRMatrixNumCols(C_offd),
1855  hypre_CSRMatrixNumNonzeros(C_diag),
1856  hypre_CSRMatrixNumNonzeros(C_offd));
1857 
1858  /* In C, destroy diag/offd (allocated by Create) and replace them with
1859  C_diag/C_offd. */
1860  hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(C));
1861  hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(C));
1862  hypre_ParCSRMatrixDiag(C) = C_diag;
1863  hypre_ParCSRMatrixOffd(C) = C_offd;
1864 
1865  hypre_ParCSRMatrixColMapOffd(C) = C_cmap;
1866  }
1867  else
1868  {
1869  /* A and B have different column mappings for their off-diagonal blocks so
1870  we need to use the column maps to create full-width CSR matrices. */
1871 
1872  int ierr = 0;
1873  hypre_CSRMatrix * csr_A;
1874  hypre_CSRMatrix * csr_B;
1875  hypre_CSRMatrix * csr_C_temp;
1876 
1877  /* merge diag and off-diag portions of A */
1878  csr_A = hypre_MergeDiagAndOffd(A);
1879 
1880  /* merge diag and off-diag portions of B */
1881  csr_B = hypre_MergeDiagAndOffd(B);
1882 
1883  /* add A and B */
1884  csr_C_temp = hypre_CSRMatrixAdd(csr_A, csr_B);
1885 
1886  /* delete CSR versions of A and B */
1887  ierr += hypre_CSRMatrixDestroy(csr_A);
1888  ierr += hypre_CSRMatrixDestroy(csr_B);
1889 
1890  /* create a new empty ParCSR matrix to contain the sum */
1891  C = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
1892  hypre_ParCSRMatrixGlobalNumRows(A),
1893  hypre_ParCSRMatrixGlobalNumCols(A),
1894  hypre_ParCSRMatrixRowStarts(A),
1895  hypre_ParCSRMatrixColStarts(A),
1896  0, 0, 0);
1897 
1898  /* split C into diag and off-diag portions */
1899  /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
1900  number of columns in csr_C_temp which is the global number of columns
1901  in A and B. This does not scale well. */
1902  ierr += GenerateDiagAndOffd(csr_C_temp, C,
1903  hypre_ParCSRMatrixFirstColDiag(A),
1904  hypre_ParCSRMatrixLastColDiag(A));
1905 
1906  /* delete CSR version of C */
1907  ierr += hypre_CSRMatrixDestroy(csr_C_temp);
1908 
1909  MFEM_VERIFY(ierr == 0, "");
1910  }
1911 
1912  /* hypre_ParCSRMatrixSetNumNonzeros(A); */
1913 
1914  /* Make sure that the first entry in each row is the diagonal one. */
1915  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(C));
1916 
1917  /* C owns diag, offd, and cmap. */
1918  hypre_ParCSRMatrixSetDataOwner(C, 1);
1919  /* C does not own row and column starts. */
1920  hypre_ParCSRMatrixSetRowStartsOwner(C, 0);
1921  hypre_ParCSRMatrixSetColStartsOwner(C, 0);
1922 
1923  return C;
1924 }
1925 
1926 HYPRE_Int
1927 hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
1928  HYPRE_Complex beta,
1929  hypre_ParCSRMatrix *B)
1930 {
1931  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1932  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1933  hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1934  hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1935  HYPRE_Int ncols_B_offd = hypre_CSRMatrixNumCols(B_offd);
1936  HYPRE_Int error;
1937 
1938  error = hypre_CSRMatrixSum(A_diag, beta, B_diag);
1939  if (ncols_B_offd > 0) /* treat B_offd as zero if it has no columns */
1940  {
1941  error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
1942  }
1943 
1944  return error;
1945 }
1946 
1947 HYPRE_Int
1948 hypre_CSRMatrixSetConstantValues(hypre_CSRMatrix *A,
1949  HYPRE_Complex value)
1950 {
1951  HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1952  HYPRE_Int A_nnz = hypre_CSRMatrixNumNonzeros(A);
1953  HYPRE_Int ia;
1954 
1955  for (ia = 0; ia < A_nnz; ia++)
1956  {
1957  A_data[ia] = value;
1958  }
1959 
1960  return 0;
1961 }
1962 
1963 HYPRE_Int
1964 hypre_ParCSRMatrixSetConstantValues(hypre_ParCSRMatrix *A,
1965  HYPRE_Complex value)
1966 {
1967  internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixDiag(A), value);
1968  internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixOffd(A), value);
1969 
1970  return 0;
1971 }
1972 
1973 } // namespace mfem::internal
1974 
1975 } // namespace mfem
1976 
1977 #endif // MFEM_USE_MPI
hypre_parcsr.hpp
HYPRE_BigInt
HYPRE_Int HYPRE_BigInt
Definition: hypre_parcsr.hpp:28
a
double a
Definition: lissajous.cpp:41
index
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
alpha
const double alpha
Definition: ex15.cpp:369
Generated on Thu Jul 29 2021 15:58:39 for MFEM by   doxygen 1.8.5