MFEM  v4.5.1
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-2022, 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 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 */
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 #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
705 void 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 
720 void 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 
818 void 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 */
1021 void 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 */
1132 void 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 */
1211 void 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 */
1319 void 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. */
1382 hypre_ParCSRCommHandle *
1383 hypre_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 */
1470 void 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 */
1528 void 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 */
1611 void 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 */
1669 void 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 
1751 HYPRE_Int
1752 hypre_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 
1803 hypre_ParCSRMatrix *
1804 hypre_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 
1946 HYPRE_Int
1947 hypre_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 
1967 HYPRE_Int
1968 hypre_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 
1983 HYPRE_Int
1984 hypre_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
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:235
alpha
const double alpha
Definition: ex15.cpp:369
Generated on Wed Dec 7 2022 10:41:47 for MFEM by   doxygen 1.8.5