MFEM  v3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "hypre_parcsr.hpp"
17 #include <limits>
18 
19 namespace mfem
20 {
21 namespace internal
22 {
23 
24 /*--------------------------------------------------------------------------
25  * A*X = B style elimination
26  *--------------------------------------------------------------------------*/
27 
28 /*
29  Function: hypre_CSRMatrixEliminateAXB
30 
31  Eliminate the rows and columns of A corresponding to the
32  given sorted (!) list of rows. Put I on the eliminated diagonal,
33  subtract columns times X from B.
34 */
35 void hypre_CSRMatrixEliminateAXB(hypre_CSRMatrix *A,
36  HYPRE_Int nrows_to_eliminate,
37  HYPRE_Int *rows_to_eliminate,
38  hypre_Vector *X,
39  hypre_Vector *B)
40 {
41  HYPRE_Int i, j;
42  HYPRE_Int irow, jcol, ibeg, iend, pos;
43  HYPRE_Real a;
44 
45  HYPRE_Int *Ai = hypre_CSRMatrixI(A);
46  HYPRE_Int *Aj = hypre_CSRMatrixJ(A);
47  HYPRE_Real *Adata = hypre_CSRMatrixData(A);
48  HYPRE_Int nrows = hypre_CSRMatrixNumRows(A);
49 
50  HYPRE_Real *Xdata = hypre_VectorData(X);
51  HYPRE_Real *Bdata = hypre_VectorData(B);
52 
53  /* eliminate the columns */
54  for (i = 0; i < nrows; i++)
55  {
56  ibeg = Ai[i];
57  iend = Ai[i+1];
58  for (j = ibeg; j < iend; j++)
59  {
60  jcol = Aj[j];
61  pos = hypre_BinarySearch(rows_to_eliminate, jcol, nrows_to_eliminate);
62  if (pos != -1)
63  {
64  a = Adata[j];
65  Adata[j] = 0.0;
66  Bdata[i] -= a * Xdata[jcol];
67  }
68  }
69  }
70 
71  /* remove the rows and set the diagonal equal to 1 */
72  for (i = 0; i < nrows_to_eliminate; i++)
73  {
74  irow = rows_to_eliminate[i];
75  ibeg = Ai[irow];
76  iend = Ai[irow+1];
77  for (j = ibeg; j < iend; j++)
78  {
79  if (Aj[j] == irow)
80  {
81  Adata[j] = 1.0;
82  }
83  else
84  {
85  Adata[j] = 0.0;
86  }
87  }
88  }
89 }
90 
91 /*
92  Function: hypre_CSRMatrixEliminateOffdColsAXB
93 
94  Eliminate the given sorted (!) list of columns of A, subtract them from B.
95 */
96 void hypre_CSRMatrixEliminateOffdColsAXB(hypre_CSRMatrix *A,
97  HYPRE_Int ncols_to_eliminate,
98  HYPRE_Int *eliminate_cols,
99  HYPRE_Real *eliminate_coefs,
100  hypre_Vector *B)
101 {
102  HYPRE_Int i, j;
103  HYPRE_Int ibeg, iend, pos;
104  HYPRE_Real a;
105 
106  HYPRE_Int *Ai = hypre_CSRMatrixI(A);
107  HYPRE_Int *Aj = hypre_CSRMatrixJ(A);
108  HYPRE_Real *Adata = hypre_CSRMatrixData(A);
109  HYPRE_Int nrows = hypre_CSRMatrixNumRows(A);
110 
111  HYPRE_Real *Bdata = hypre_VectorData(B);
112 
113  for (i = 0; i < nrows; i++)
114  {
115  ibeg = Ai[i];
116  iend = Ai[i+1];
117  for (j = ibeg; j < iend; j++)
118  {
119  pos = hypre_BinarySearch(eliminate_cols, Aj[j], ncols_to_eliminate);
120  if (pos != -1)
121  {
122  a = Adata[j];
123  Adata[j] = 0.0;
124  Bdata[i] -= a * eliminate_coefs[pos];
125  }
126  }
127  }
128 }
129 
130 /*
131  Function: hypre_CSRMatrixEliminateOffdRowsAXB
132 
133  Eliminate (zero) the given list of rows of A.
134 */
135 void hypre_CSRMatrixEliminateOffdRowsAXB(hypre_CSRMatrix *A,
136  HYPRE_Int nrows_to_eliminate,
137  HYPRE_Int *rows_to_eliminate)
138 {
139  HYPRE_Int *Ai = hypre_CSRMatrixI(A);
140  HYPRE_Real *Adata = hypre_CSRMatrixData(A);
141 
142  HYPRE_Int i, j;
143  HYPRE_Int irow, ibeg, iend;
144 
145  for (i = 0; i < nrows_to_eliminate; i++)
146  {
147  irow = rows_to_eliminate[i];
148  ibeg = Ai[irow];
149  iend = Ai[irow+1];
150  for (j = ibeg; j < iend; j++)
151  {
152  Adata[j] = 0.0;
153  }
154  }
155 }
156 
157 /*
158  Function: hypre_ParCSRMatrixEliminateAXB
159 
160  This function eliminates the global rows and columns of a matrix
161  A corresponding to given lists of sorted (!) local row numbers,
162  so that the solution to the system A*X = B is X_b for the given rows.
163 
164  The elimination is done as follows:
165 
166  (input) (output)
167 
168  / A_ii | A_ib \ / A_ii | 0 \
169  A = | -----+----- | ---> | -----+----- |
170  \ A_bi | A_bb / \ 0 | I /
171 
172  / X_i \ / X_i \
173  X = | --- | ---> | --- | (no change)
174  \ X_b / \ X_b /
175 
176  / B_i \ / B_i - A_ib * X_b \
177  B = | --- | ---> | ---------------- |
178  \ B_b / \ X_b /
179 
180 */
181 void hypre_ParCSRMatrixEliminateAXB(hypre_ParCSRMatrix *A,
182  HYPRE_Int num_rowscols_to_elim,
183  HYPRE_Int *rowscols_to_elim,
184  hypre_ParVector *X,
185  hypre_ParVector *B)
186 {
187  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
188  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
189  HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows(diag);
190  HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
191 
192  hypre_Vector *Xlocal = hypre_ParVectorLocalVector(X);
193  hypre_Vector *Blocal = hypre_ParVectorLocalVector(B);
194 
195  HYPRE_Real *Bdata = hypre_VectorData(Blocal);
196  HYPRE_Real *Xdata = hypre_VectorData(Xlocal);
197 
198  HYPRE_Int num_offd_cols_to_elim;
199  HYPRE_Int *offd_cols_to_elim;
200  HYPRE_Real *eliminate_coefs;
201 
202  /* figure out which offd cols should be eliminated and with what coef */
203  hypre_ParCSRCommHandle *comm_handle;
204  hypre_ParCSRCommPkg *comm_pkg;
205  HYPRE_Int num_sends;
206  HYPRE_Int index, start;
207  HYPRE_Int i, j, k, irow;
208 
209  HYPRE_Real *eliminate_row = hypre_CTAlloc(HYPRE_Real, diag_nrows);
210  HYPRE_Real *eliminate_col = hypre_CTAlloc(HYPRE_Real, offd_ncols);
211  HYPRE_Real *buf_data, coef;
212 
213  /* make sure A has a communication package */
214  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
215  if (!comm_pkg)
216  {
217  hypre_MatvecCommPkgCreate(A);
218  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
219  }
220 
221  /* HACK: rows that shouldn't be eliminated are marked with quiet NaN;
222  those that should are set to the boundary value from X; this is to
223  avoid sending complex type (int+double) or communicating twice. */
224  for (i = 0; i < diag_nrows; i++)
225  {
226  eliminate_row[i] = std::numeric_limits<HYPRE_Real>::quiet_NaN();
227  }
228  for (i = 0; i < num_rowscols_to_elim; i++)
229  {
230  irow = rowscols_to_elim[i];
231  eliminate_row[irow] = Xdata[irow];
232  }
233 
234  /* use a Matvec communication pattern to find (in eliminate_col)
235  which of the local offd columns are to be eliminated */
236  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
237  buf_data = hypre_CTAlloc(HYPRE_Real,
238  hypre_ParCSRCommPkgSendMapStart(comm_pkg,
239  num_sends));
240  index = 0;
241  for (i = 0; i < num_sends; i++)
242  {
243  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
244  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
245  {
246  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
247  buf_data[index++] = eliminate_row[k];
248  }
249  }
250  comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg,
251  buf_data, eliminate_col);
252 
253  /* do sequential part of the elimination while stuff is getting sent */
254  hypre_CSRMatrixEliminateAXB(diag, num_rowscols_to_elim, rowscols_to_elim,
255  Xlocal, Blocal);
256 
257  /* finish the communication */
258  hypre_ParCSRCommHandleDestroy(comm_handle);
259 
260  /* received eliminate_col[], count offd columns to eliminate */
261  num_offd_cols_to_elim = 0;
262  for (i = 0; i < offd_ncols; i++)
263  {
264  coef = eliminate_col[i];
265  if (coef == coef) // test for NaN
266  {
267  num_offd_cols_to_elim++;
268  }
269  }
270 
271  offd_cols_to_elim = hypre_CTAlloc(HYPRE_Int, num_offd_cols_to_elim);
272  eliminate_coefs = hypre_CTAlloc(HYPRE_Real, num_offd_cols_to_elim);
273 
274  /* get a list of offd column indices and coefs */
275  num_offd_cols_to_elim = 0;
276  for (i = 0; i < offd_ncols; i++)
277  {
278  coef = eliminate_col[i];
279  if (coef == coef) // test for NaN
280  {
281  offd_cols_to_elim[num_offd_cols_to_elim] = i;
282  eliminate_coefs[num_offd_cols_to_elim] = coef;
283  num_offd_cols_to_elim++;
284  }
285  }
286 
287  hypre_TFree(buf_data);
288  hypre_TFree(eliminate_row);
289  hypre_TFree(eliminate_col);
290 
291  /* eliminate the off-diagonal part */
292  hypre_CSRMatrixEliminateOffdColsAXB(offd, num_offd_cols_to_elim,
293  offd_cols_to_elim,
294  eliminate_coefs, Blocal);
295 
296  hypre_CSRMatrixEliminateOffdRowsAXB(offd, num_rowscols_to_elim,
297  rowscols_to_elim);
298 
299  /* set boundary values in the rhs */
300  for (int i = 0; i < num_rowscols_to_elim; i++)
301  {
302  irow = rowscols_to_elim[i];
303  Bdata[irow] = Xdata[irow];
304  }
305 
306  hypre_TFree(offd_cols_to_elim);
307  hypre_TFree(eliminate_coefs);
308 }
309 
310 
311 /*--------------------------------------------------------------------------
312  * (A + Ae) style elimination
313  *--------------------------------------------------------------------------*/
314 
315 /*
316  Function: hypre_CSRMatrixElimCreate
317 
318  Prepare the Ae matrix: count nnz, initialize I, allocate J and data.
319 */
320 void hypre_CSRMatrixElimCreate(hypre_CSRMatrix *A,
321  hypre_CSRMatrix *Ae,
322  HYPRE_Int nrows, HYPRE_Int *rows,
323  HYPRE_Int ncols, HYPRE_Int *cols,
324  HYPRE_Int *col_mark)
325 {
326  HYPRE_Int i, j, col;
327  HYPRE_Int A_beg, A_end;
328 
329  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
330  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
331  HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
332 
333  hypre_CSRMatrixI(Ae) = hypre_TAlloc(HYPRE_Int, A_rows+1);
334 
335  HYPRE_Int *Ae_i = hypre_CSRMatrixI(Ae);
336  HYPRE_Int nnz = 0;
337 
338  for (i = 0; i < A_rows; i++)
339  {
340  Ae_i[i] = nnz;
341 
342  A_beg = A_i[i];
343  A_end = A_i[i+1];
344 
345  if (hypre_BinarySearch(rows, i, nrows) >= 0)
346  {
347  /* full row */
348  nnz += A_end - A_beg;
349 
350  if (col_mark)
351  {
352  for (j = A_beg; j < A_end; j++)
353  {
354  col_mark[A_j[j]] = 1;
355  }
356  }
357  }
358  else
359  {
360  /* count columns */
361  for (j = A_beg; j < A_end; j++)
362  {
363  col = A_j[j];
364  if (hypre_BinarySearch(cols, col, ncols) >= 0)
365  {
366  nnz++;
367  if (col_mark) { col_mark[col] = 1; }
368  }
369  }
370  }
371  }
372  Ae_i[A_rows] = nnz;
373 
374  hypre_CSRMatrixJ(Ae) = hypre_TAlloc(HYPRE_Int, nnz);
375  hypre_CSRMatrixData(Ae) = hypre_TAlloc(HYPRE_Real, nnz);
376  hypre_CSRMatrixNumNonzeros(Ae) = nnz;
377 }
378 
379 /*
380  Function: hypre_CSRMatrixEliminateRowsCols
381 
382  Eliminate rows and columns of A, store eliminated values in Ae.
383  If 'diag' is nonzero, the eliminated diagonal of A is set to identity.
384  If 'col_remap' is not NULL it specifies renumbering of columns of Ae.
385 */
386 void hypre_CSRMatrixEliminateRowsCols(hypre_CSRMatrix *A,
387  hypre_CSRMatrix *Ae,
388  HYPRE_Int nrows, HYPRE_Int *rows,
389  HYPRE_Int ncols, HYPRE_Int *cols,
390  int diag, HYPRE_Int* col_remap)
391 {
392  HYPRE_Int i, j, k, col;
393  HYPRE_Int A_beg, Ae_beg, A_end;
394  HYPRE_Real a;
395 
396  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
397  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
398  HYPRE_Real *A_data = hypre_CSRMatrixData(A);
399  HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
400 
401  HYPRE_Int *Ae_i = hypre_CSRMatrixI(Ae);
402  HYPRE_Int *Ae_j = hypre_CSRMatrixJ(Ae);
403  HYPRE_Real *Ae_data = hypre_CSRMatrixData(Ae);
404 
405  for (i = 0; i < A_rows; i++)
406  {
407  A_beg = A_i[i];
408  A_end = A_i[i+1];
409  Ae_beg = Ae_i[i];
410 
411  if (hypre_BinarySearch(rows, i, nrows) >= 0)
412  {
413  /* eliminate row */
414  for (j = A_beg, k = Ae_beg; j < A_end; j++, k++)
415  {
416  col = A_j[j];
417  Ae_j[k] = col_remap ? col_remap[col] : col;
418  a = (diag && col == i) ? 1.0 : 0.0;
419  Ae_data[k] = A_data[j] - a;
420  A_data[j] = a;
421  }
422  }
423  else
424  {
425  /* eliminate columns */
426  for (j = A_beg, k = Ae_beg; j < A_end; j++)
427  {
428  col = A_j[j];
429  if (hypre_BinarySearch(cols, col, ncols) >= 0)
430  {
431  Ae_j[k] = col_remap ? col_remap[col] : col;
432  Ae_data[k] = A_data[j];
433  A_data[j] = 0.0;
434  k++;
435  }
436  }
437  }
438  }
439 }
440 
441 
442 /*
443  Function: hypre_ParCSRMatrixEliminateAAe
444 
445  (input) (output)
446 
447  / A_ii | A_ib \ / A_ii | 0 \
448  A = | -----+----- | ---> | -----+----- |
449  \ A_bi | A_bb / \ 0 | I /
450 
451 
452  / 0 | A_ib \
453  Ae = | -----+--------- |
454  \ A_bi | A_bb - I /
455 
456 */
457 void hypre_ParCSRMatrixEliminateAAe(hypre_ParCSRMatrix *A,
458  hypre_ParCSRMatrix **Ae,
459  HYPRE_Int num_rowscols_to_elim,
460  HYPRE_Int *rowscols_to_elim)
461 {
462  HYPRE_Int i, j, k;
463 
464  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
465  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
466  HYPRE_Int A_diag_nrows = hypre_CSRMatrixNumRows(A_diag);
467  HYPRE_Int A_offd_ncols = hypre_CSRMatrixNumCols(A_offd);
468 
469  *Ae = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
470  hypre_ParCSRMatrixGlobalNumRows(A),
471  hypre_ParCSRMatrixGlobalNumCols(A),
472  hypre_ParCSRMatrixRowStarts(A),
473  hypre_ParCSRMatrixColStarts(A),
474  0, 0, 0);
475 
476  hypre_ParCSRMatrixSetRowStartsOwner(*Ae, 0);
477  hypre_ParCSRMatrixSetColStartsOwner(*Ae, 0);
478 
479  hypre_CSRMatrix *Ae_diag = hypre_ParCSRMatrixDiag(*Ae);
480  hypre_CSRMatrix *Ae_offd = hypre_ParCSRMatrixOffd(*Ae);
481  HYPRE_Int Ae_offd_ncols;
482 
483  HYPRE_Int num_offd_cols_to_elim;
484  HYPRE_Int *offd_cols_to_elim;
485 
486  HYPRE_Int *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
487  HYPRE_Int *Ae_col_map_offd;
488 
489  HYPRE_Int *col_mark;
490  HYPRE_Int *col_remap;
491 
492  /* figure out which offd cols should be eliminated */
493  {
494  hypre_ParCSRCommHandle *comm_handle;
495  hypre_ParCSRCommPkg *comm_pkg;
496  HYPRE_Int num_sends, *int_buf_data;
497  HYPRE_Int index, start;
498 
499  HYPRE_Int *eliminate_row = hypre_CTAlloc(HYPRE_Int, A_diag_nrows);
500  HYPRE_Int *eliminate_col = hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
501 
502  /* make sure A has a communication package */
503  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
504  if (!comm_pkg)
505  {
506  hypre_MatvecCommPkgCreate(A);
507  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
508  }
509 
510  /* which of the local rows are to be eliminated */
511  for (i = 0; i < A_diag_nrows; i++)
512  {
513  eliminate_row[i] = 0;
514  }
515  for (i = 0; i < num_rowscols_to_elim; i++)
516  {
517  eliminate_row[rowscols_to_elim[i]] = 1;
518  }
519 
520  /* use a Matvec communication pattern to find (in eliminate_col)
521  which of the local offd columns are to be eliminated */
522  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
523  int_buf_data = hypre_CTAlloc(HYPRE_Int,
524  hypre_ParCSRCommPkgSendMapStart(comm_pkg,
525  num_sends));
526  index = 0;
527  for (i = 0; i < num_sends; i++)
528  {
529  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
530  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
531  {
532  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
533  int_buf_data[index++] = eliminate_row[k];
534  }
535  }
536  comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
537  int_buf_data, eliminate_col);
538 
539  /* eliminate diagonal part, overlapping it with communication */
540  hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
541  num_rowscols_to_elim, rowscols_to_elim,
542  num_rowscols_to_elim, rowscols_to_elim,
543  NULL);
544 
545  hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
546  num_rowscols_to_elim, rowscols_to_elim,
547  num_rowscols_to_elim, rowscols_to_elim,
548  1, NULL);
549  hypre_CSRMatrixReorder(Ae_diag);
550 
551  /* finish the communication */
552  hypre_ParCSRCommHandleDestroy(comm_handle);
553 
554  /* received eliminate_col[], count offd columns to eliminate */
555  num_offd_cols_to_elim = 0;
556  for (i = 0; i < A_offd_ncols; i++)
557  {
558  if (eliminate_col[i]) { num_offd_cols_to_elim++; }
559  }
560 
561  offd_cols_to_elim = hypre_CTAlloc(HYPRE_Int, num_offd_cols_to_elim);
562 
563  /* get a list of offd column indices and coefs */
564  num_offd_cols_to_elim = 0;
565  for (i = 0; i < A_offd_ncols; i++)
566  {
567  if (eliminate_col[i])
568  {
569  offd_cols_to_elim[num_offd_cols_to_elim++] = i;
570  }
571  }
572 
573  hypre_TFree(int_buf_data);
574  hypre_TFree(eliminate_row);
575  hypre_TFree(eliminate_col);
576  }
577 
578  /* eliminate the off-diagonal part */
579  col_mark = hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
580  col_remap = hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
581 
582  hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
583  num_rowscols_to_elim, rowscols_to_elim,
584  num_offd_cols_to_elim, offd_cols_to_elim,
585  col_mark);
586 
587  for (i = k = 0; i < A_offd_ncols; i++)
588  {
589  if (col_mark[i]) { col_remap[i] = k++; }
590  }
591 
592  hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
593  num_rowscols_to_elim, rowscols_to_elim,
594  num_offd_cols_to_elim, offd_cols_to_elim,
595  0, col_remap);
596 
597  /* create col_map_offd for Ae */
598  Ae_offd_ncols = 0;
599  for (i = 0; i < A_offd_ncols; i++)
600  {
601  if (col_mark[i]) { Ae_offd_ncols++; }
602  }
603 
604  Ae_col_map_offd = hypre_CTAlloc(HYPRE_Int, Ae_offd_ncols);
605 
606  Ae_offd_ncols = 0;
607  for (i = 0; i < A_offd_ncols; i++)
608  {
609  if (col_mark[i])
610  {
611  Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
612  }
613  }
614 
615  hypre_ParCSRMatrixColMapOffd(*Ae) = Ae_col_map_offd;
616  hypre_CSRMatrixNumCols(Ae_offd) = Ae_offd_ncols;
617 
618  hypre_TFree(col_remap);
619  hypre_TFree(col_mark);
620  hypre_TFree(offd_cols_to_elim);
621 
622  hypre_ParCSRMatrixSetNumNonzeros(*Ae);
623  hypre_MatvecCommPkgCreate(*Ae);
624 }
625 
626 
627 /*--------------------------------------------------------------------------
628  * Split
629  *--------------------------------------------------------------------------*/
630 
631 void hypre_CSRMatrixSplit(hypre_CSRMatrix *A,
632  HYPRE_Int nr, HYPRE_Int nc,
633  HYPRE_Int *row_block_num, HYPRE_Int *col_block_num,
634  hypre_CSRMatrix **blocks)
635 {
636  HYPRE_Int i, j, k, bi, bj;
637 
638  HYPRE_Int* A_i = hypre_CSRMatrixI(A);
639  HYPRE_Int* A_j = hypre_CSRMatrixJ(A);
640  HYPRE_Complex* A_data = hypre_CSRMatrixData(A);
641 
642  HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
643  HYPRE_Int A_cols = hypre_CSRMatrixNumCols(A);
644 
645  HYPRE_Int *num_rows = hypre_CTAlloc(HYPRE_Int, nr);
646  HYPRE_Int *num_cols = hypre_CTAlloc(HYPRE_Int, nc);
647 
648  HYPRE_Int *block_row = hypre_TAlloc(HYPRE_Int, A_rows);
649  HYPRE_Int *block_col = hypre_TAlloc(HYPRE_Int, A_cols);
650 
651  for (i = 0; i < A_rows; i++)
652  {
653  block_row[i] = num_rows[row_block_num[i]]++;
654  }
655  for (j = 0; j < A_cols; j++)
656  {
657  block_col[j] = num_cols[col_block_num[j]]++;
658  }
659 
660  /* allocate the blocks */
661  for (i = 0; i < nr; i++)
662  {
663  for (j = 0; j < nc; j++)
664  {
665  hypre_CSRMatrix *B = hypre_CSRMatrixCreate(num_rows[i], num_cols[j], 0);
666  hypre_CSRMatrixI(B) = hypre_CTAlloc(HYPRE_Int, num_rows[i] + 1);
667  blocks[i*nc + j] = B;
668  }
669  }
670 
671  /* count block row nnz */
672  for (i = 0; i < A_rows; i++)
673  {
674  bi = row_block_num[i];
675  for (j = A_i[i]; j < A_i[i+1]; j++)
676  {
677  bj = col_block_num[A_j[j]];
678  hypre_CSRMatrix *B = blocks[bi*nc + bj];
679  hypre_CSRMatrixI(B)[block_row[i] + 1]++;
680  }
681  }
682 
683  /* count block nnz */
684  for (k = 0; k < nr*nc; k++)
685  {
686  hypre_CSRMatrix *B = blocks[k];
687  HYPRE_Int* B_i = hypre_CSRMatrixI(B);
688 
689  HYPRE_Int nnz = 0, rs;
690  for (int k = 1; k <= hypre_CSRMatrixNumRows(B); k++)
691  {
692  rs = B_i[k], B_i[k] = nnz, nnz += rs;
693  }
694 
695  hypre_CSRMatrixJ(B) = hypre_TAlloc(HYPRE_Int, nnz);
696  hypre_CSRMatrixData(B) = hypre_TAlloc(HYPRE_Complex, nnz);
697  hypre_CSRMatrixNumNonzeros(B) = nnz;
698  }
699 
700  /* populate blocks */
701  for (i = 0; i < A_rows; i++)
702  {
703  bi = row_block_num[i];
704  for (j = A_i[i]; j < A_i[i+1]; j++)
705  {
706  k = A_j[j];
707  bj = col_block_num[k];
708  hypre_CSRMatrix *B = blocks[bi*nc + bj];
709  HYPRE_Int *bii = hypre_CSRMatrixI(B) + block_row[i] + 1;
710  hypre_CSRMatrixJ(B)[*bii] = block_col[k];
711  hypre_CSRMatrixData(B)[*bii] = A_data[j];
712  (*bii)++;
713  }
714  }
715 
716  hypre_TFree(block_col);
717  hypre_TFree(block_row);
718 
719  hypre_TFree(num_cols);
720  hypre_TFree(num_rows);
721 }
722 
723 
724 void hypre_ParCSRMatrixSplit(hypre_ParCSRMatrix *A,
725  HYPRE_Int nr, HYPRE_Int nc,
726  hypre_ParCSRMatrix **blocks,
727  int interleaved_rows, int interleaved_cols)
728 {
729  HYPRE_Int i, j, k;
730 
731  MPI_Comm comm = hypre_ParCSRMatrixComm(A);
732 
733  hypre_CSRMatrix *Adiag = hypre_ParCSRMatrixDiag(A);
734  hypre_CSRMatrix *Aoffd = hypre_ParCSRMatrixOffd(A);
735 
736  HYPRE_Int global_rows = hypre_ParCSRMatrixGlobalNumRows(A);
737  HYPRE_Int global_cols = hypre_ParCSRMatrixGlobalNumCols(A);
738 
739  HYPRE_Int local_rows = hypre_CSRMatrixNumRows(Adiag);
740  HYPRE_Int local_cols = hypre_CSRMatrixNumCols(Adiag);
741  HYPRE_Int offd_cols = hypre_CSRMatrixNumCols(Aoffd);
742 
743  hypre_assert(local_rows % nr == 0 && local_cols % nc == 0);
744  hypre_assert(global_rows % nr == 0 && global_cols % nc == 0);
745 
746  HYPRE_Int block_rows = local_rows / nr;
747  HYPRE_Int block_cols = local_cols / nc;
748  HYPRE_Int num_blocks = nr * nc;
749 
750  /* mark local rows and columns with block number */
751  HYPRE_Int *row_block_num = hypre_TAlloc(HYPRE_Int, local_rows);
752  HYPRE_Int *col_block_num = hypre_TAlloc(HYPRE_Int, local_cols);
753 
754  for (i = 0; i < local_rows; i++)
755  {
756  row_block_num[i] = interleaved_rows ? (i % nr) : (i / block_rows);
757  }
758  for (i = 0; i < local_cols; i++)
759  {
760  col_block_num[i] = interleaved_cols ? (i % nc) : (i / block_cols);
761  }
762 
763  /* determine the block numbers for offd columns */
764  HYPRE_Int* offd_col_block_num = hypre_TAlloc(HYPRE_Int, offd_cols);
765  hypre_ParCSRCommHandle *comm_handle;
766  HYPRE_Int *int_buf_data;
767  {
768  /* make sure A has a communication package */
769  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
770  if (!comm_pkg)
771  {
772  hypre_MatvecCommPkgCreate(A);
773  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
774  }
775 
776  /* calculate the final global column numbers for each block */
777  HYPRE_Int *count = hypre_CTAlloc(HYPRE_Int, nc);
778  HYPRE_Int *block_global_col = hypre_TAlloc(HYPRE_Int, local_cols);
779  HYPRE_Int first_col = hypre_ParCSRMatrixFirstColDiag(A) / nc;
780  for (i = 0; i < local_cols; i++)
781  {
782  block_global_col[i] = first_col + count[col_block_num[i]]++;
783  }
784  hypre_TFree(count);
785 
786  /* use a Matvec communication pattern to determine offd_col_block_num */
787  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
788  int_buf_data = hypre_CTAlloc(HYPRE_Int,
789  hypre_ParCSRCommPkgSendMapStart(comm_pkg,
790  num_sends));
791  HYPRE_Int start, index = 0;
792  for (i = 0; i < num_sends; i++)
793  {
794  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
795  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
796  {
797  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
798  int_buf_data[index++] = col_block_num[k] + nc*block_global_col[k];
799  }
800  }
801  hypre_TFree(block_global_col);
802 
803  comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
804  offd_col_block_num);
805  }
806 
807  /* create the block matrices */
808  HYPRE_Int num_procs = 1;
809  if (!hypre_ParCSRMatrixAssumedPartition(A))
810  {
811  hypre_MPI_Comm_size(comm, &num_procs);
812  }
813 
814  HYPRE_Int *row_starts = hypre_TAlloc(HYPRE_Int, num_procs+1);
815  HYPRE_Int *col_starts = hypre_TAlloc(HYPRE_Int, num_procs+1);
816  for (i = 0; i <= num_procs; i++)
817  {
818  row_starts[i] = hypre_ParCSRMatrixRowStarts(A)[i] / nr;
819  col_starts[i] = hypre_ParCSRMatrixColStarts(A)[i] / nc;
820  }
821 
822  for (i = 0; i < num_blocks; i++)
823  {
824  blocks[i] = hypre_ParCSRMatrixCreate(comm,
825  global_rows/nr, global_cols/nc,
826  row_starts, col_starts, 0, 0, 0);
827  }
828 
829  /* split diag part */
830  hypre_CSRMatrix **csr_blocks = hypre_TAlloc(hypre_CSRMatrix*, nr*nc);
831  hypre_CSRMatrixSplit(Adiag, nr, nc, row_block_num, col_block_num,
832  csr_blocks);
833 
834  for (i = 0; i < num_blocks; i++)
835  {
836  hypre_TFree(hypre_ParCSRMatrixDiag(blocks[i]));
837  hypre_ParCSRMatrixDiag(blocks[i]) = csr_blocks[i];
838  }
839 
840  /* finish communication, receive offd_col_block_num */
841  hypre_ParCSRCommHandleDestroy(comm_handle);
842  hypre_TFree(int_buf_data);
843 
844  /* decode global offd column numbers */
845  HYPRE_Int* offd_global_col = hypre_TAlloc(HYPRE_Int, offd_cols);
846  for (i = 0; i < offd_cols; i++)
847  {
848  offd_global_col[i] = offd_col_block_num[i] / nc;
849  offd_col_block_num[i] %= nc;
850  }
851 
852  /* split offd part */
853  hypre_CSRMatrixSplit(Aoffd, nr, nc, row_block_num, offd_col_block_num,
854  csr_blocks);
855 
856  for (i = 0; i < num_blocks; i++)
857  {
858  hypre_TFree(hypre_ParCSRMatrixOffd(blocks[i]));
859  hypre_ParCSRMatrixOffd(blocks[i]) = csr_blocks[i];
860  }
861 
862  hypre_TFree(csr_blocks);
863  hypre_TFree(col_block_num);
864  hypre_TFree(row_block_num);
865 
866  /* update block col-maps */
867  for (int bi = 0; bi < nr; bi++)
868  {
869  for (int bj = 0; bj < nc; bj++)
870  {
871  hypre_ParCSRMatrix *block = blocks[bi*nc + bj];
872  hypre_CSRMatrix *block_offd = hypre_ParCSRMatrixOffd(block);
873  HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols(block_offd);
874 
875  HYPRE_Int *block_col_map = hypre_TAlloc(HYPRE_Int, block_offd_cols);
876  for (i = j = 0; i < offd_cols; i++)
877  {
878  HYPRE_Int bn = offd_col_block_num[i];
879  if (bn == bj) { block_col_map[j++] = offd_global_col[i]; }
880  }
881  hypre_assert(j == block_offd_cols);
882 
883  hypre_ParCSRMatrixColMapOffd(block) = block_col_map;
884  }
885  }
886 
887  hypre_TFree(offd_global_col);
888  hypre_TFree(offd_col_block_num);
889 
890  /* finish the new matrices, make them own all the stuff */
891  for (i = 0; i < num_blocks; i++)
892  {
893  hypre_ParCSRMatrixSetNumNonzeros(blocks[i]);
894  hypre_MatvecCommPkgCreate(blocks[i]);
895 
896  hypre_ParCSRMatrixOwnsData(blocks[i]) = 1;
897 
898  /* only the first block will own the row/col_starts */
899  hypre_ParCSRMatrixOwnsRowStarts(blocks[i]) = !i;
900  hypre_ParCSRMatrixOwnsColStarts(blocks[i]) = !i;
901  }
902 }
903 
904 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
905 void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,
906  HYPRE_Bool alpha,
907  HYPRE_Bool *x,
908  HYPRE_Bool beta,
909  HYPRE_Bool *y)
910 {
911  /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
912  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
913  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
914  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
915 
916  HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
917  HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
918 
919  HYPRE_Bool *x_data = x;
920  HYPRE_Bool *y_data = y;
921 
922  HYPRE_Bool temp, tempx;
923 
924  HYPRE_Int i, jj;
925 
926  HYPRE_Int m;
927 
928  HYPRE_Real xpar=0.7;
929 
930  /*-----------------------------------------------------------------------
931  * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
932  *-----------------------------------------------------------------------*/
933 
934  if (alpha == 0)
935  {
936 #ifdef HYPRE_USING_OPENMP
937  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
938 #endif
939  for (i = 0; i < num_rows; i++)
940  {
941  y_data[i] = y_data[i] && beta;
942  }
943  return;
944  }
945 
946  /*-----------------------------------------------------------------------
947  * y = (beta/alpha)*y
948  *-----------------------------------------------------------------------*/
949 
950  if (beta == 0)
951  {
952 #ifdef HYPRE_USING_OPENMP
953  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
954 #endif
955  for (i = 0; i < num_rows; i++)
956  {
957  y_data[i] = 0;
958  }
959  }
960  else
961  {
962  /* beta is true -> no change to y_data */
963  }
964 
965  /*-----------------------------------------------------------------
966  * y += A*x
967  *-----------------------------------------------------------------*/
968 
969  /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows */
970 
971  if (num_rownnz < xpar*(num_rows))
972  {
973 #ifdef HYPRE_USING_OPENMP
974  #pragma omp parallel for private(i,jj,m,tempx) HYPRE_SMP_SCHEDULE
975 #endif
976  for (i = 0; i < num_rownnz; i++)
977  {
978  m = A_rownnz[i];
979 
980  tempx = 0;
981  for (jj = A_i[m]; jj < A_i[m+1]; jj++)
982  {
983  /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
984  tempx = tempx || x_data[A_j[jj]];
985  }
986  y_data[m] = y_data[m] || tempx;
987  }
988  }
989  else
990  {
991 #ifdef HYPRE_USING_OPENMP
992  #pragma omp parallel for private(i,jj,temp) HYPRE_SMP_SCHEDULE
993 #endif
994  for (i = 0; i < num_rows; i++)
995  {
996  temp = 0;
997  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
998  {
999  /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1000  temp = temp || x_data[A_j[jj]];
1001  }
1002  y_data[i] = y_data[i] || temp;
1003  }
1004  }
1005 
1006  /*-----------------------------------------------------------------
1007  * y = alpha*y
1008  *-----------------------------------------------------------------*/
1009  /* alpha is true */
1010 }
1011 
1012 /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
1013  input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
1014 hypre_ParCSRCommHandle *
1015 hypre_ParCSRCommHandleCreate_bool(HYPRE_Int job,
1016  hypre_ParCSRCommPkg *comm_pkg,
1017  HYPRE_Bool *send_data,
1018  HYPRE_Bool *recv_data)
1019 {
1020  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1021  HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
1022  MPI_Comm comm = hypre_ParCSRCommPkgComm(comm_pkg);
1023 
1024  hypre_ParCSRCommHandle *comm_handle;
1025  HYPRE_Int num_requests;
1026  hypre_MPI_Request *requests;
1027 
1028  HYPRE_Int i, j;
1029  HYPRE_Int my_id, num_procs;
1030  HYPRE_Int ip, vec_start, vec_len;
1031 
1032  num_requests = num_sends + num_recvs;
1033  requests = hypre_CTAlloc(hypre_MPI_Request, num_requests);
1034 
1035  hypre_MPI_Comm_size(comm, &num_procs);
1036  hypre_MPI_Comm_rank(comm, &my_id);
1037 
1038  j = 0;
1039  switch (job)
1040  {
1041  case 1:
1042  {
1043  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1044  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1045  for (i = 0; i < num_recvs; i++)
1046  {
1047  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1048  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1049  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1050  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1051  ip, 0, comm, &requests[j++]);
1052  }
1053  for (i = 0; i < num_sends; i++)
1054  {
1055  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1056  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1057  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1058  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1059  ip, 0, comm, &requests[j++]);
1060  }
1061  break;
1062  }
1063  case 2:
1064  {
1065  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1066  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1067  for (i = 0; i < num_sends; i++)
1068  {
1069  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1070  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1071  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1072  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1073  ip, 0, comm, &requests[j++]);
1074  }
1075  for (i = 0; i < num_recvs; i++)
1076  {
1077  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1078  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1079  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1080  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1081  ip, 0, comm, &requests[j++]);
1082  }
1083  break;
1084  }
1085  }
1086  /*--------------------------------------------------------------------
1087  * set up comm_handle and return
1088  *--------------------------------------------------------------------*/
1089 
1090  comm_handle = hypre_CTAlloc(hypre_ParCSRCommHandle, 1);
1091 
1092  hypre_ParCSRCommHandleCommPkg(comm_handle) = comm_pkg;
1093  hypre_ParCSRCommHandleSendData(comm_handle) = send_data;
1094  hypre_ParCSRCommHandleRecvData(comm_handle) = recv_data;
1095  hypre_ParCSRCommHandleNumRequests(comm_handle) = num_requests;
1096  hypre_ParCSRCommHandleRequests(comm_handle) = requests;
1097 
1098  return comm_handle;
1099 }
1100 
1101 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1102 void hypre_ParCSRMatrixBooleanMatvec(hypre_ParCSRMatrix *A,
1103  HYPRE_Bool alpha,
1104  HYPRE_Bool *x,
1105  HYPRE_Bool beta,
1106  HYPRE_Bool *y)
1107 {
1108  hypre_ParCSRCommHandle *comm_handle;
1109  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1110  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1111  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1112 
1113  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1114  HYPRE_Int num_sends, i, j, index;
1115 
1116  HYPRE_Bool *x_tmp, *x_buf;
1117 
1118  x_tmp = hypre_CTAlloc(HYPRE_Bool, num_cols_offd);
1119 
1120  /*---------------------------------------------------------------------
1121  * If there exists no CommPkg for A, a CommPkg is generated using
1122  * equally load balanced partitionings
1123  *--------------------------------------------------------------------*/
1124  if (!comm_pkg)
1125  {
1126  hypre_MatvecCommPkgCreate(A);
1127  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1128  }
1129 
1130  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1131  x_buf = hypre_CTAlloc(HYPRE_Bool,
1132  hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1133 
1134  index = 0;
1135  for (i = 0; i < num_sends; i++)
1136  {
1137  j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1138  for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1139  {
1140  x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1141  }
1142  }
1143 
1144  comm_handle = hypre_ParCSRCommHandleCreate_bool(1, comm_pkg, x_buf, x_tmp);
1145 
1146  hypre_CSRMatrixBooleanMatvec(diag, alpha, x, beta, y);
1147 
1148  hypre_ParCSRCommHandleDestroy(comm_handle);
1149 
1150  if (num_cols_offd)
1151  {
1152  hypre_CSRMatrixBooleanMatvec(offd, alpha, x_tmp, 1, y);
1153  }
1154 
1155  hypre_TFree(x_buf);
1156  hypre_TFree(x_tmp);
1157 }
1158 
1159 } // namespace mfem::internal
1160 
1161 } // namespace mfem
1162 
1163 #endif // MFEM_USE_MPI
hypre_parcsr.hpp
alpha
const double alpha
Definition: ex15.cpp:337
Generated on Fri Feb 18 2022 19:00:29 for MFEM by   doxygen 1.8.5