MFEM  v4.5.2
Finite element discretization library
sparsemat.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // Implementation of sparse matrix
13 
14 #include "linalg.hpp"
15 #include "../general/forall.hpp"
16 #include "../general/table.hpp"
17 #include "../general/sort_pairs.hpp"
18 #include "../general/backends.hpp"
19 
20 #include <iostream>
21 #include <iomanip>
22 #include <cmath>
23 #include <algorithm>
24 #include <limits>
25 #include <cstring>
26 
27 #if defined(MFEM_USE_CUDA)
28 #define MFEM_cu_or_hip(stub) cu##stub
29 #define MFEM_Cu_or_Hip(stub) Cu##stub
30 #define MFEM_CU_or_HIP(stub) CU##stub
31 #define MFEM_CUDA_or_HIP(stub) CUDA##stub
32 
33 #if CUSPARSE_VERSION >= 11400
34 #define MFEM_GPUSPARSE_ALG CUSPARSE_SPMV_CSR_ALG1
35 #else // CUSPARSE_VERSION >= 11400
36 #define MFEM_GPUSPARSE_ALG CUSPARSE_CSRMV_ALG1
37 #endif // CUSPARSE_VERSION >= 11400
38 
39 #elif defined(MFEM_USE_HIP)
40 #define MFEM_cu_or_hip(stub) hip##stub
41 #define MFEM_Cu_or_Hip(stub) Hip##stub
42 #define MFEM_CU_or_HIP(stub) HIP##stub
43 #define MFEM_CUDA_or_HIP(stub) HIP##stub
44 
45 // https://hipsparse.readthedocs.io/en/latest/usermanual.html#hipsparsespmvalg-t
46 #define MFEM_GPUSPARSE_ALG HIPSPARSE_CSRMV_ALG1
47 #endif // defined(MFEM_USE_CUDA)
48 
49 namespace mfem
50 {
51 
52 using namespace std;
53 
54 #ifdef MFEM_USE_CUDA_OR_HIP
56 // doxygen doesn't like the macro-assisted typename so let's skip parsing it:
57 // \cond false
58 MFEM_cu_or_hip(sparseHandle_t) SparseMatrix::handle = nullptr;
59 // \endcond
60 size_t SparseMatrix::bufferSize = 0;
61 void * SparseMatrix::dBuffer = nullptr;
62 #endif // MFEM_USE_CUDA_OR_HIP
63 
65 {
66  // Initialize cuSPARSE/hipSPARSE library
67 #ifdef MFEM_USE_CUDA_OR_HIP
69  {
70  if (!handle) { MFEM_cu_or_hip(sparseCreate)(&handle); }
71  useGPUSparse=true;
72  SparseMatrixCount++;
73  }
74  else
75  {
76  useGPUSparse=false;
77  }
78 #endif // MFEM_USE_CUDA_OR_HIP
79 }
80 
82 {
83 #ifdef MFEM_USE_CUDA_OR_HIP
84  if (initBuffers)
85  {
86 #if CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
87  MFEM_cu_or_hip(sparseDestroySpMat)(matA_descr);
88  MFEM_cu_or_hip(sparseDestroyDnVec)(vecX_descr);
89  MFEM_cu_or_hip(sparseDestroyDnVec)(vecY_descr);
90 #else
91  cusparseDestroyMatDescr(matA_descr);
92 #endif // CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
93  initBuffers = false;
94  }
95 #endif // MFEM_USE_CUDA_OR_HIP
96 }
97 
98 SparseMatrix::SparseMatrix(int nrows, int ncols)
99  : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
100  Rows(new RowNode *[nrows]),
101  current_row(-1),
102  ColPtrJ(NULL),
103  ColPtrNode(NULL),
104  At(NULL),
105  isSorted(false)
106 {
107  // We probably do not need to set the ownership flags here.
108  I.SetHostPtrOwner(true);
109  J.SetHostPtrOwner(true);
110  A.SetHostPtrOwner(true);
111 
112  for (int i = 0; i < nrows; i++)
113  {
114  Rows[i] = NULL;
115  }
116 
117 #ifdef MFEM_USE_MEMALLOC
118  NodesMem = new RowNodeAlloc;
119 #endif
120 
121  InitGPUSparse();
122 }
123 
124 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
125  : AbstractSparseMatrix(m, n),
126  Rows(NULL),
127  ColPtrJ(NULL),
128  ColPtrNode(NULL),
129  At(NULL),
130  isSorted(false)
131 {
132  I.Wrap(i, height+1, true);
133  J.Wrap(j, I[height], true);
134  A.Wrap(data, I[height], true);
135 
136 #ifdef MFEM_USE_MEMALLOC
137  NodesMem = NULL;
138 #endif
139 
140  InitGPUSparse();
141 }
142 
143 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
144  bool ownij, bool owna, bool issorted)
145  : AbstractSparseMatrix(m, n),
146  Rows(NULL),
147  ColPtrJ(NULL),
148  ColPtrNode(NULL),
149  At(NULL),
150  isSorted(issorted)
151 {
152  I.Wrap(i, height+1, ownij);
153  J.Wrap(j, I[height], ownij);
154 
155 #ifdef MFEM_USE_MEMALLOC
156  NodesMem = NULL;
157 #endif
158 
159  if (data)
160  {
161  A.Wrap(data, I[height], owna);
162  }
163  else
164  {
165  const int nnz = I[height];
166  A.New(nnz);
167  for (int ii=0; ii<nnz; ++ii)
168  {
169  A[ii] = 0.0;
170  }
171  }
172 
173  InitGPUSparse();
174 }
175 
176 SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
177  : AbstractSparseMatrix(nrows, ncols)
178  , Rows(NULL)
179  , ColPtrJ(NULL)
180  , ColPtrNode(NULL)
181  , At(NULL)
182  , isSorted(false)
183 {
184 #ifdef MFEM_USE_MEMALLOC
185  NodesMem = NULL;
186 #endif
187  I.New(nrows + 1);
188  J.New(nrows * rowsize);
189  A.New(nrows * rowsize);
190 
191  for (int i = 0; i <= nrows; i++)
192  {
193  I[i] = i * rowsize;
194  }
195 
196  InitGPUSparse();
197 }
198 
199 SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph,
200  MemoryType mt)
201  : AbstractSparseMatrix(mat.Height(), mat.Width())
202 {
203  if (mat.Finalized())
204  {
205  mat.HostReadI();
206  const int nnz = mat.I[height];
207  if (copy_graph)
208  {
209  I.New(height+1, mt == MemoryType::PRESERVE ? mat.I.GetMemoryType() : mt);
210  J.New(nnz, mt == MemoryType::PRESERVE ? mat.J.GetMemoryType() : mt);
211  I.CopyFrom(mat.I, height+1);
212  J.CopyFrom(mat.J, nnz);
213  }
214  else
215  {
216  I = mat.I;
217  J = mat.J;
218  I.ClearOwnerFlags();
219  J.ClearOwnerFlags();
220  }
221  A.New(nnz, mt == MemoryType::PRESERVE ? mat.A.GetMemoryType() : mt);
222  A.CopyFrom(mat.A, nnz);
223 
224  Rows = NULL;
225 #ifdef MFEM_USE_MEMALLOC
226  NodesMem = NULL;
227 #endif
228  }
229  else
230  {
231 #ifdef MFEM_USE_MEMALLOC
232  NodesMem = new RowNodeAlloc;
233 #endif
234  Rows = new RowNode *[height];
235  for (int i = 0; i < height; i++)
236  {
237  RowNode **node_pp = &Rows[i];
238  for (RowNode *node_p = mat.Rows[i]; node_p; node_p = node_p->Prev)
239  {
240 #ifdef MFEM_USE_MEMALLOC
241  RowNode *new_node_p = NodesMem->Alloc();
242 #else
243  RowNode *new_node_p = new RowNode;
244 #endif
245  new_node_p->Value = node_p->Value;
246  new_node_p->Column = node_p->Column;
247  *node_pp = new_node_p;
248  node_pp = &new_node_p->Prev;
249  }
250  *node_pp = NULL;
251  }
252 
253  // We probably do not need to set the ownership flags here.
254  I.SetHostPtrOwner(true);
255  J.SetHostPtrOwner(true);
256  A.SetHostPtrOwner(true);
257  }
258 
259  current_row = -1;
260  ColPtrJ = NULL;
261  ColPtrNode = NULL;
262  At = NULL;
263  isSorted = mat.isSorted;
264 
265  InitGPUSparse();
266 }
267 
269  : AbstractSparseMatrix(v.Size(), v.Size())
270  , Rows(NULL)
271  , ColPtrJ(NULL)
272  , ColPtrNode(NULL)
273  , At(NULL)
274  , isSorted(true)
275 {
276 #ifdef MFEM_USE_MEMALLOC
277  NodesMem = NULL;
278 #endif
279  I.New(height + 1);
280  J.New(height);
281  A.New(height);
282 
283  for (int i = 0; i <= height; i++)
284  {
285  I[i] = i;
286  }
287 
288  for (int r=0; r<height; r++)
289  {
290  J[r] = r;
291  A[r] = v[r];
292  }
293 
294  InitGPUSparse();
295 }
296 
297 void SparseMatrix::OverrideSize(int height_, int width_)
298 {
299  height = height_;
300  width = width_;
301 }
302 
304 {
305  Clear();
306 
307  SparseMatrix copy(rhs);
308  Swap(copy);
309 
310  return *this;
311 }
312 
314 {
315  MFEM_ASSERT(master.Finalized(), "'master' must be finalized");
316  Clear();
317  height = master.Height();
318  width = master.Width();
319  I = master.I; I.ClearOwnerFlags();
320  J = master.J; J.ClearOwnerFlags();
321  A = master.A; A.ClearOwnerFlags();
322  isSorted = master.isSorted;
323 }
324 
326 {
327  height = width = 0;
328  I.Reset();
329  J.Reset();
330  A.Reset();
331  Rows = NULL;
332  current_row = -1;
333  ColPtrJ = NULL;
334  ColPtrNode = NULL;
335  At = NULL;
336 #ifdef MFEM_USE_MEMALLOC
337  NodesMem = NULL;
338 #endif
339  isSorted = false;
340 
341  ClearGPUSparse();
342 }
343 
344 int SparseMatrix::RowSize(const int i) const
345 {
346  int gi = i;
347  if (gi < 0)
348  {
349  gi = -1-gi;
350  }
351 
352  if (I)
353  {
354  return I[gi+1]-I[gi];
355  }
356 
357  int s = 0;
358  RowNode *row = Rows[gi];
359  for ( ; row != NULL; row = row->Prev)
360  if (row->Value != 0.0)
361  {
362  s++;
363  }
364  return s;
365 }
366 
368 {
369  int max_row_size=0;
370  int rowSize=0;
371  if (I)
372  {
373  for (int i=0; i < height; ++i)
374  {
375  rowSize = I[i+1]-I[i];
376  max_row_size = (max_row_size > rowSize) ? max_row_size : rowSize;
377  }
378  }
379  else
380  {
381  for (int i=0; i < height; ++i)
382  {
383  rowSize = RowSize(i);
384  max_row_size = (max_row_size > rowSize) ? max_row_size : rowSize;
385  }
386  }
387 
388  return max_row_size;
389 }
390 
391 int *SparseMatrix::GetRowColumns(const int row)
392 {
393  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
394 
395  return J + I[row];
396 }
397 
398 const int *SparseMatrix::GetRowColumns(const int row) const
399 {
400  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
401 
402  return J + I[row];
403 }
404 
405 double *SparseMatrix::GetRowEntries(const int row)
406 {
407  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
408 
409  return A + I[row];
410 }
411 
412 const double *SparseMatrix::GetRowEntries(const int row) const
413 {
414  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
415 
416  return A + I[row];
417 }
418 
419 void SparseMatrix::SetWidth(int newWidth)
420 {
421  if (newWidth == width)
422  {
423  // Nothing to be done here
424  return;
425  }
426  else if (newWidth == -1)
427  {
428  // Compute the actual width
429  width = ActualWidth();
430  // No need to reset the ColPtr, since the new ColPtr will be shorter.
431  }
432  else if (newWidth > width)
433  {
434  // We need to reset ColPtr, since now we may have additional columns.
435  if (Rows != NULL)
436  {
437  delete [] ColPtrNode;
438  ColPtrNode = static_cast<RowNode **>(NULL);
439  }
440  else
441  {
442  delete [] ColPtrJ;
443  ColPtrJ = static_cast<int *>(NULL);
444  }
445  width = newWidth;
446  }
447  else
448  {
449  // Check that the new width is bigger or equal to the actual width.
450  MFEM_ASSERT(newWidth >= ActualWidth(),
451  "The new width needs to be bigger or equal to the actual width");
452  width = newWidth;
453  }
454 }
455 
456 
458 {
459  MFEM_VERIFY(Finalized(), "Matrix is not Finalized!");
460 
461  if (isSorted)
462  {
463  return;
464  }
465 
466 #ifdef MFEM_USE_CUDA_OR_HIP
468  {
469 #if defined(MFEM_USE_CUDA)
470  size_t pBufferSizeInBytes = 0;
471  void *pBuffer = NULL;
472 
473  const int n = Height();
474  const int m = Width();
475  const int nnzA = J.Capacity();
476  double * d_a_sorted = ReadWriteData();
477  const int * d_ia = ReadI();
478  int * d_ja_sorted = ReadWriteJ();
479  csru2csrInfo_t sortInfoA;
480 
481  cusparseMatDescr_t matA_descr;
482  cusparseCreateMatDescr( &matA_descr );
483  cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
484  cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
485 
486  cusparseCreateCsru2csrInfo( &sortInfoA );
487 
488  cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
489  d_ja_sorted, sortInfoA,
490  &pBufferSizeInBytes);
491 
492  CuMemAlloc( &pBuffer, pBufferSizeInBytes );
493 
494  cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
495  d_ja_sorted, sortInfoA, pBuffer);
496 
497  // The above call is (at least in some cases) asynchronous, so we need to
498  // wait for it to finish before we can free device temporaries.
499  MFEM_STREAM_SYNC;
500 
501  cusparseDestroyCsru2csrInfo( sortInfoA );
502  cusparseDestroyMatDescr( matA_descr );
503 
504  CuMemFree( pBuffer );
505 #endif
506  }
507  else if ( Device::Allows( Backend::HIP_MASK ))
508  {
509 #if defined(MFEM_USE_HIP)
510  size_t pBufferSizeInBytes = 0;
511  void *pBuffer = NULL;
512  int *P = NULL;
513 
514  const int n = Height();
515  const int m = Width();
516  const int nnzA = J.Capacity();
517  double * d_a_sorted = ReadWriteData();
518  const int * d_ia = ReadI();
519  int * d_ja_sorted = ReadWriteJ();
520 
521  hipsparseMatDescr_t descrA;
522  hipsparseCreateMatDescr( &descrA );
523  // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
524  // a temporary copy of the data for gthr, sort that, and then copy the sorted values
525  // back to the array being returned. Where there is an in-place version available,
526  // we should use it.
527  Array< double > a_tmp( nnzA );
528  double *d_a_tmp = a_tmp.Write();
529 
530  hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
531  &pBufferSizeInBytes);
532 
533  HipMemAlloc( &pBuffer, pBufferSizeInBytes );
534  HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
535 
536  hipsparseCreateIdentityPermutation(handle, nnzA, P);
537  hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
538 
539  hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
540  HIPSPARSE_INDEX_BASE_ZERO);
541 
542  A.CopyFrom( a_tmp.GetMemory(), nnzA );
543  hipsparseDestroyMatDescr( descrA );
544 
545  HipMemFree( pBuffer );
546  HipMemFree( P );
547 #endif
548  }
549  else
550 #endif // MFEM_USE_CUDA_OR_HIP
551  {
552  const int * Ip=HostReadI();
553  HostReadWriteJ();
555 
557  for (int j = 0, i = 0; i < height; i++)
558  {
559  int end = Ip[i+1];
560  row.SetSize(end - j);
561  for (int k = 0; k < row.Size(); k++)
562  {
563  row[k].one = J[j+k];
564  row[k].two = A[j+k];
565  }
566  row.Sort();
567  for (int k = 0; k < row.Size(); k++, j++)
568  {
569  J[j] = row[k].one;
570  A[j] = row[k].two;
571  }
572  }
573  }
574  isSorted = true;
575 }
576 
578 {
579  MFEM_VERIFY(Finalized(), "Matrix is not Finalized!");
580 
581  for (int row = 0, end = 0; row < height; row++)
582  {
583  int start = end, j;
584  end = I[row+1];
585  for (j = start; true; j++)
586  {
587  MFEM_VERIFY(j < end, "diagonal entry not found in row = " << row);
588  if (J[j] == row) { break; }
589  }
590  const double diag = A[j];
591  for ( ; j > start; j--)
592  {
593  J[j] = J[j-1];
594  A[j] = A[j-1];
595  }
596  J[start] = row;
597  A[start] = diag;
598  }
599 }
600 
601 double &SparseMatrix::Elem(int i, int j)
602 {
603  return operator()(i,j);
604 }
605 
606 const double &SparseMatrix::Elem(int i, int j) const
607 {
608  return operator()(i,j);
609 }
610 
611 double &SparseMatrix::operator()(int i, int j)
612 {
613  MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
614  "Trying to access element outside of the matrix. "
615  << "height = " << height << ", "
616  << "width = " << width << ", "
617  << "i = " << i << ", "
618  << "j = " << j);
619 
620  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
621 
622  for (int k = I[i], end = I[i+1]; k < end; k++)
623  {
624  if (J[k] == j)
625  {
626  return A[k];
627  }
628  }
629 
630  MFEM_ABORT("Did not find i = " << i << ", j = " << j << " in matrix.");
631  return A[0];
632 }
633 
634 const double &SparseMatrix::operator()(int i, int j) const
635 {
636  static const double zero = 0.0;
637 
638  MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
639  "Trying to access element outside of the matrix. "
640  << "height = " << height << ", "
641  << "width = " << width << ", "
642  << "i = " << i << ", "
643  << "j = " << j);
644 
645  if (Finalized())
646  {
647  for (int k = I[i], end = I[i+1]; k < end; k++)
648  {
649  if (J[k] == j)
650  {
651  return A[k];
652  }
653  }
654  }
655  else
656  {
657  for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
658  {
659  if (node_p->Column == j)
660  {
661  return node_p->Value;
662  }
663  }
664  }
665 
666  return zero;
667 }
668 
670 {
671  MFEM_VERIFY(height == width, "Matrix must be square, not height = "
672  << height << ", width = " << width);
673  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
674 
675  d.SetSize(height);
676 
677  const auto II = this->ReadI();
678  const auto JJ = this->ReadJ();
679  const auto AA = this->ReadData();
680  auto dd = d.Write();
681 
682  MFEM_FORALL(i, height,
683  {
684  const int begin = II[i];
685  const int end = II[i+1];
686  int j;
687  for (j = begin; j < end; j++)
688  {
689  if (JJ[j] == i)
690  {
691  dd[i] = AA[j];
692  break;
693  }
694  }
695  if (j == end)
696  {
697  dd[i] = 0.;
698  }
699  });
700 }
701 
702 /// Produces a DenseMatrix from a SparseMatrix
704 {
705  int num_rows = this->Height();
706  int num_cols = this->Width();
707 
708  DenseMatrix * B = new DenseMatrix(num_rows, num_cols);
709 
710  this->ToDenseMatrix(*B);
711 
712  return B;
713 }
714 
715 /// Produces a DenseMatrix from a SparseMatrix
717 {
718  B.SetSize(height, width);
719  B = 0.0;
720 
721  for (int r=0; r<height; r++)
722  {
723  const int * col = this->GetRowColumns(r);
724  const double * val = this->GetRowEntries(r);
725 
726  for (int cj=0; cj<this->RowSize(r); cj++)
727  {
728  B(r, col[cj]) = val[cj];
729  }
730  }
731 }
732 
733 void SparseMatrix::Mult(const Vector &x, Vector &y) const
734 {
735  if (Finalized()) { y.UseDevice(true); }
736  y = 0.0;
737  AddMult(x, y);
738 }
739 
740 void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
741 {
742  MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
743  << ") must match matrix width (" << width << ")");
744  MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
745  << ") must match matrix height (" << height << ")");
746 
747  if (!Finalized())
748  {
749  const double *xp = x.HostRead();
750  double *yp = y.HostReadWrite();
751 
752  // The matrix is not finalized, but multiplication is still possible
753  for (int i = 0; i < height; i++)
754  {
755  RowNode *row = Rows[i];
756  double b = 0.0;
757  for ( ; row != NULL; row = row->Prev)
758  {
759  b += row->Value * xp[row->Column];
760  }
761  *yp += a * b;
762  yp++;
763  }
764  return;
765  }
766 
767 #ifndef MFEM_USE_LEGACY_OPENMP
768  const int height = this->height;
769  const int nnz = J.Capacity();
770  auto d_I = Read(I, height+1);
771  auto d_J = Read(J, nnz);
772  auto d_A = Read(A, nnz);
773  auto d_x = x.Read();
774  auto d_y = y.ReadWrite();
775 
776  // Skip if matrix has no non-zeros
777  if (nnz == 0) {return;}
779  {
780 #ifdef MFEM_USE_CUDA_OR_HIP
781  const double alpha = a;
782  const double beta = 1.0;
783 
784  // Setup descriptors
785  if (!initBuffers)
786  {
787 #if CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
788  // Setup matrix descriptor
789  MFEM_cu_or_hip(sparseCreateCsr)(
790  &matA_descr,Height(),
791  Width(),
792  J.Capacity(),
793  const_cast<int *>(d_I),
794  const_cast<int *>(d_J),
795  const_cast<double *>(d_A),
796  MFEM_CU_or_HIP(SPARSE_INDEX_32I),
797  MFEM_CU_or_HIP(SPARSE_INDEX_32I),
798  MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
799  MFEM_CUDA_or_HIP(_R_64F));
800 
801  // Create handles for input/output vectors
802  MFEM_cu_or_hip(sparseCreateDnVec)(&vecX_descr,
803  x.Size(),
804  const_cast<double *>(d_x),
805  MFEM_CUDA_or_HIP(_R_64F));
806  MFEM_cu_or_hip(sparseCreateDnVec)(&vecY_descr, y.Size(), d_y,
807  MFEM_CUDA_or_HIP(_R_64F));
808 #else
809  cusparseCreateMatDescr(&matA_descr);
810  cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
811  cusparseSetMatType(matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL);
812 #endif // CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
813  initBuffers = true;
814  }
815  // Allocate kernel space. Buffer is shared between different sparsemats
816  size_t newBufferSize = 0;
817 
818  MFEM_cu_or_hip(sparseSpMV_bufferSize)(
819  handle,
820  MFEM_CU_or_HIP(SPARSE_OPERATION_NON_TRANSPOSE),
821  &alpha,
822  matA_descr,
823  vecX_descr,
824  &beta,
825  vecY_descr,
826  MFEM_CUDA_or_HIP(_R_64F),
827  MFEM_GPUSPARSE_ALG,
828  &newBufferSize);
829 
830  // Check if we need to resize
831  if (newBufferSize > bufferSize)
832  {
833  bufferSize = newBufferSize;
834  if (dBuffer != nullptr) { MFEM_Cu_or_Hip(MemFree)(dBuffer); }
835  MFEM_Cu_or_Hip(MemAlloc)(&dBuffer, bufferSize);
836  }
837 
838 #if CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
839  // Update input/output vectors
840  MFEM_cu_or_hip(sparseDnVecSetValues)(vecX_descr,
841  const_cast<double *>(d_x));
842  MFEM_cu_or_hip(sparseDnVecSetValues)(vecY_descr, d_y);
843 
844  // Y = alpha A * X + beta * Y
845  MFEM_cu_or_hip(sparseSpMV)(
846  handle,
847  MFEM_CU_or_HIP(SPARSE_OPERATION_NON_TRANSPOSE),
848  &alpha,
849  matA_descr,
850  vecX_descr,
851  &beta,
852  vecY_descr,
853  MFEM_CUDA_or_HIP(_R_64F),
854  MFEM_GPUSPARSE_ALG,
855  dBuffer);
856 #else
857  cusparseDcsrmv(handle,
858  CUSPARSE_OPERATION_NON_TRANSPOSE,
859  Height(),
860  Width(),
861  J.Capacity(),
862  &alpha,
863  matA_descr,
864  const_cast<double *>(d_A),
865  const_cast<int *>(d_I),
866  const_cast<int *>(d_J),
867  const_cast<double *>(d_x),
868  &beta,
869  d_y);
870 #endif // CUDA_VERSION >= 10010 || defined(MFEM_USE_HIP)
871 #endif // MFEM_USE_CUDA_OR_HIP
872  }
873  else
874  {
875  // Native version
876  MFEM_FORALL(i, height,
877  {
878  double d = 0.0;
879  const int end = d_I[i+1];
880  for (int j = d_I[i]; j < end; j++)
881  {
882  d += d_A[j] * d_x[d_J[j]];
883  }
884  d_y[i] += a * d;
885  });
886 
887  }
888 
889 #else // MFEM_USE_LEGACY_OPENMP
890  const double *Ap = A, *xp = x.GetData();
891  double *yp = y.GetData();
892  const int *Jp = J, *Ip = I;
893 
894  #pragma omp parallel for
895  for (int i = 0; i < height; i++)
896  {
897  double d = 0.0;
898  const int end = Ip[i+1];
899  for (int j = Ip[i]; j < end; j++)
900  {
901  d += Ap[j] * xp[Jp[j]];
902  }
903  yp[i] += a * d;
904  }
905 #endif // MFEM_USE_LEGACY_OPENMP
906 }
907 
908 void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
909 {
910  if (Finalized()) { y.UseDevice(true); }
911  y = 0.0;
912  AddMultTranspose(x, y);
913 }
914 
916  const double a) const
917 {
918  MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
919  << ") must match matrix height (" << height << ")");
920  MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
921  << ") must match matrix width (" << width << ")");
922 
923  if (!Finalized())
924  {
925  double *yp = y.HostReadWrite();
926  const double *xp = x.HostRead();
927  // The matrix is not finalized, but multiplication is still possible
928  for (int i = 0; i < height; i++)
929  {
930  RowNode *row = Rows[i];
931  double b = a * xp[i];
932  for ( ; row != NULL; row = row->Prev)
933  {
934  yp[row->Column] += row->Value * b;
935  }
936  }
937  return;
938  }
939 
941  if (At)
942  {
943  At->AddMult(x, y, a);
944  }
945  else
946  {
947  for (int i = 0; i < height; i++)
948  {
949  const double xi = a * x[i];
950  const int end = I[i+1];
951  for (int j = I[i]; j < end; j++)
952  {
953  const int Jj = J[j];
954  y[Jj] += A[j] * xi;
955  }
956  }
957  }
958 }
959 
961 {
962  if (At == NULL)
963  {
964  At = Transpose(*this);
965  }
966 }
967 
969 {
970  delete At;
971  At = NULL;
972 }
973 
975 {
977  {
978  BuildTranspose();
979  }
980 }
981 
983  const Array<int> &rows, const Vector &x, Vector &y) const
984 {
985  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
986 
987  const int n = rows.Size();
988  const int nnz = J.Capacity();
989  auto d_rows = rows.Read();
990  auto d_I = Read(I, height+1);
991  auto d_J = Read(J, nnz);
992  auto d_A = Read(A, nnz);
993  auto d_x = x.Read();
994  auto d_y = y.Write();
995  MFEM_FORALL(i, n,
996  {
997  const int r = d_rows[i];
998  const int end = d_I[r + 1];
999  double a = 0.0;
1000  for (int j = d_I[r]; j < end; j++)
1001  {
1002  a += d_A[j] * d_x[d_J[j]];
1003  }
1004  d_y[r] = a;
1005  });
1006 }
1007 
1009  const Array<int> &rows, const Vector &x, Vector &y, const double a) const
1010 {
1011  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
1012 
1013  for (int i = 0; i < rows.Size(); i++)
1014  {
1015  int r = rows[i];
1016  int end = I[r + 1];
1017  double val = 0.0;
1018  for (int j = I[r]; j < end; j++)
1019  {
1020  val += A[j] * x(J[j]);
1021  }
1022  y(r) += a * val;
1023  }
1024 }
1025 
1027 {
1028  MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
1029  MFEM_ASSERT(x.Size() == Width(), "Input vector size (" << x.Size()
1030  << ") must match matrix width (" << Width() << ")");
1031 
1033 
1034  const int height = Height();
1035  const int nnz = J.Capacity();
1036  auto d_I = Read(I, height+1);
1037  auto d_J = Read(J, nnz);
1038  auto d_x = Read(x.GetMemory(), x.Size());
1039  auto d_y = Write(y.GetMemory(), y.Size());
1040  MFEM_FORALL(i, height,
1041  {
1042  bool d_yi = false;
1043  const int end = d_I[i+1];
1044  for (int j = d_I[i]; j < end; j++)
1045  {
1046  if (d_x[d_J[j]])
1047  {
1048  d_yi = true;
1049  break;
1050  }
1051  }
1052  d_y[i] = d_yi;
1053  });
1054 }
1055 
1057  Array<int> &y) const
1058 {
1059  MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
1060  MFEM_ASSERT(x.Size() == Height(), "Input vector size (" << x.Size()
1061  << ") must match matrix height (" << Height() << ")");
1062 
1063  y.SetSize(Width());
1064  y = 0;
1065 
1066  for (int i = 0; i < Height(); i++)
1067  {
1068  if (x[i])
1069  {
1070  int end = I[i+1];
1071  for (int j = I[i]; j < end; j++)
1072  {
1073  y[J[j]] = x[i];
1074  }
1075  }
1076  }
1077 }
1078 
1079 void SparseMatrix::AbsMult(const Vector &x, Vector &y) const
1080 {
1081  MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
1082  << ") must match matrix width (" << width << ")");
1083  MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
1084  << ") must match matrix height (" << height << ")");
1085 
1086  if (Finalized()) { y.UseDevice(true); }
1087  y = 0.0;
1088 
1089  if (!Finalized())
1090  {
1091  const double *xp = x.HostRead();
1092  double *yp = y.HostReadWrite();
1093 
1094  // The matrix is not finalized, but multiplication is still possible
1095  for (int i = 0; i < height; i++)
1096  {
1097  RowNode *row = Rows[i];
1098  double b = 0.0;
1099  for ( ; row != NULL; row = row->Prev)
1100  {
1101  b += std::abs(row->Value) * xp[row->Column];
1102  }
1103  *yp += b;
1104  yp++;
1105  }
1106  return;
1107  }
1108 
1109  const int height = this->height;
1110  const int nnz = J.Capacity();
1111  auto d_I = Read(I, height+1);
1112  auto d_J = Read(J, nnz);
1113  auto d_A = Read(A, nnz);
1114  auto d_x = x.Read();
1115  auto d_y = y.ReadWrite();
1116  MFEM_FORALL(i, height,
1117  {
1118  double d = 0.0;
1119  const int end = d_I[i+1];
1120  for (int j = d_I[i]; j < end; j++)
1121  {
1122  d += std::abs(d_A[j]) * d_x[d_J[j]];
1123  }
1124  d_y[i] += d;
1125  });
1126 }
1127 
1129 {
1130  MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
1131  << ") must match matrix height (" << height << ")");
1132  MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
1133  << ") must match matrix width (" << width << ")");
1134 
1135  y = 0.0;
1136 
1137  if (!Finalized())
1138  {
1139  double *yp = y.GetData();
1140  // The matrix is not finalized, but multiplication is still possible
1141  for (int i = 0; i < height; i++)
1142  {
1143  RowNode *row = Rows[i];
1144  double b = x(i);
1145  for ( ; row != NULL; row = row->Prev)
1146  {
1147  yp[row->Column] += fabs(row->Value) * b;
1148  }
1149  }
1150  return;
1151  }
1152 
1154  if (At)
1155  {
1156  At->AbsMult(x, y);
1157  }
1158  else
1159  {
1160  for (int i = 0; i < height; i++)
1161  {
1162  const double xi = x[i];
1163  const int end = I[i+1];
1164  for (int j = I[i]; j < end; j++)
1165  {
1166  const int Jj = J[j];
1167  y[Jj] += std::abs(A[j]) * xi;
1168  }
1169  }
1170  }
1171 }
1172 
1173 double SparseMatrix::InnerProduct(const Vector &x, const Vector &y) const
1174 {
1175  MFEM_ASSERT(x.Size() == Width(), "x.Size() = " << x.Size()
1176  << " must be equal to Width() = " << Width());
1177  MFEM_ASSERT(y.Size() == Height(), "y.Size() = " << y.Size()
1178  << " must be equal to Height() = " << Height());
1179 
1180  x.HostRead();
1181  y.HostRead();
1182  if (Finalized())
1183  {
1184  const int nnz = J.Capacity();
1185  HostRead(I, height+1);
1186  HostRead(J, nnz);
1187  HostRead(A, nnz);
1188  }
1189 
1190  double prod = 0.0;
1191  for (int i = 0; i < height; i++)
1192  {
1193  double a = 0.0;
1194  if (A)
1195  {
1196  for (int j = I[i], end = I[i+1]; j < end; j++)
1197  {
1198  a += A[j] * x(J[j]);
1199  }
1200  }
1201  else
1202  {
1203  for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
1204  {
1205  a += np->Value * x(np->Column);
1206  }
1207  }
1208  prod += a * y(i);
1209  }
1210 
1211  return prod;
1212 }
1213 
1215 {
1216  for (int i = 0; i < height; i++)
1217  {
1218  double a = 0.0;
1219  if (A)
1220  {
1221  for (int j = I[i], end = I[i+1]; j < end; j++)
1222  {
1223  a += A[j];
1224  }
1225  }
1226  else
1227  {
1228  for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
1229  {
1230  a += np->Value;
1231  }
1232  }
1233  x(i) = a;
1234  }
1235 }
1236 
1237 double SparseMatrix::GetRowNorml1(int irow) const
1238 {
1239  MFEM_VERIFY(irow < height,
1240  "row " << irow << " not in matrix with height " << height);
1241 
1242  double a = 0.0;
1243  if (A)
1244  {
1245  for (int j = I[irow], end = I[irow+1]; j < end; j++)
1246  {
1247  a += fabs(A[j]);
1248  }
1249  }
1250  else
1251  {
1252  for (RowNode *np = Rows[irow]; np != NULL; np = np->Prev)
1253  {
1254  a += fabs(np->Value);
1255  }
1256  }
1257 
1258  return a;
1259 }
1260 
1261 void SparseMatrix::Threshold(double tol, bool fix_empty_rows)
1262 {
1263  MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
1264  double atol;
1265  atol = std::abs(tol);
1266 
1267  fix_empty_rows = height == width ? fix_empty_rows : false;
1268 
1269  double *newA;
1270  int *newI, *newJ;
1271  int i, j, nz;
1272 
1273  newI = Memory<int>(height+1);
1274  newI[0] = 0;
1275  for (i = 0, nz = 0; i < height; i++)
1276  {
1277  bool found = false;
1278  for (j = I[i]; j < I[i+1]; j++)
1279  if (std::abs(A[j]) > atol)
1280  {
1281  found = true;
1282  nz++;
1283  }
1284  if (fix_empty_rows && !found) { nz++; }
1285  newI[i+1] = nz;
1286  }
1287 
1288  newJ = Memory<int>(nz);
1289  newA = Memory<double>(nz);
1290  // Assume we're sorted until we find out otherwise
1291  isSorted = true;
1292  for (i = 0, nz = 0; i < height; i++)
1293  {
1294  bool found = false;
1295  int lastCol = -1;
1296  for (j = I[i]; j < I[i+1]; j++)
1297  if (std::abs(A[j]) > atol)
1298  {
1299  found = true;
1300  newJ[nz] = J[j];
1301  newA[nz] = A[j];
1302  if ( lastCol > newJ[nz] )
1303  {
1304  isSorted = false;
1305  }
1306  lastCol = newJ[nz];
1307  nz++;
1308  }
1309  if (fix_empty_rows && !found)
1310  {
1311  newJ[nz] = i;
1312  newA[nz] = 0.0;
1313  nz++;
1314  }
1315  }
1316  Destroy();
1317  I.Wrap(newI, height+1, true);
1318  J.Wrap(newJ, I[height], true);
1319  A.Wrap(newA, I[height], true);
1320 }
1321 
1322 void SparseMatrix::Finalize(int skip_zeros, bool fix_empty_rows)
1323 {
1324  int i, j, nr, nz;
1325  RowNode *aux;
1326 
1327  if (Finalized())
1328  {
1329  return;
1330  }
1331 
1332  delete [] ColPtrNode;
1333  ColPtrNode = NULL;
1334 
1335  I.New(height+1);
1336  I[0] = 0;
1337  for (i = 1; i <= height; i++)
1338  {
1339  nr = 0;
1340  for (aux = Rows[i-1]; aux != NULL; aux = aux->Prev)
1341  {
1342  if (skip_zeros && aux->Value == 0.0)
1343  {
1344  if (skip_zeros == 2) { continue; }
1345  if ((i-1) != aux->Column) { continue; }
1346 
1347  bool found = false;
1348  double found_val;
1349  for (RowNode *other = Rows[aux->Column]; other != NULL; other = other->Prev)
1350  {
1351  if (other->Column == (i-1))
1352  {
1353  found = true;
1354  found_val = other->Value;
1355  break;
1356  }
1357  }
1358  if (found && found_val == 0.0) { continue; }
1359 
1360  }
1361  nr++;
1362  }
1363  if (fix_empty_rows && !nr) { nr = 1; }
1364  I[i] = I[i-1] + nr;
1365  }
1366 
1367  nz = I[height];
1368  J.New(nz);
1369  A.New(nz);
1370  // Assume we're sorted until we find out otherwise
1371  isSorted = true;
1372  for (j = i = 0; i < height; i++)
1373  {
1374  int lastCol = -1;
1375  nr = 0;
1376  for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
1377  {
1378  if (skip_zeros && aux->Value == 0.0)
1379  {
1380  if (skip_zeros == 2) { continue; }
1381  if (i != aux->Column) { continue; }
1382 
1383  bool found = false;
1384  double found_val;
1385  for (RowNode *other = Rows[aux->Column]; other != NULL; other = other->Prev)
1386  {
1387  if (other->Column == i)
1388  {
1389  found = true;
1390  found_val = other->Value;
1391  break;
1392  }
1393  }
1394  if (found && found_val == 0.0) { continue; }
1395  }
1396 
1397  J[j] = aux->Column;
1398  A[j] = aux->Value;
1399 
1400  if ( lastCol > J[j] )
1401  {
1402  isSorted = false;
1403  }
1404  lastCol = J[j];
1405 
1406  j++;
1407  nr++;
1408  }
1409  if (fix_empty_rows && !nr)
1410  {
1411  J[j] = i;
1412  A[j] = 1.0;
1413  j++;
1414  }
1415  }
1416 
1417 #ifdef MFEM_USE_MEMALLOC
1418  delete NodesMem;
1419  NodesMem = NULL;
1420 #else
1421  for (i = 0; i < height; i++)
1422  {
1423  RowNode *node_p = Rows[i];
1424  while (node_p != NULL)
1425  {
1426  aux = node_p;
1427  node_p = node_p->Prev;
1428  delete aux;
1429  }
1430  }
1431 #endif
1432 
1433  delete [] Rows;
1434  Rows = NULL;
1435 }
1436 
1438 {
1439  int br = blocks.NumRows(), bc = blocks.NumCols();
1440  int nr = (height + br - 1)/br, nc = (width + bc - 1)/bc;
1441 
1442  for (int j = 0; j < bc; j++)
1443  {
1444  for (int i = 0; i < br; i++)
1445  {
1446  int *bI = Memory<int>(nr + 1);
1447  for (int k = 0; k <= nr; k++)
1448  {
1449  bI[k] = 0;
1450  }
1451  blocks(i,j) = new SparseMatrix(bI, NULL, NULL, nr, nc);
1452  }
1453  }
1454 
1455  for (int gr = 0; gr < height; gr++)
1456  {
1457  int bi = gr/nr, i = gr%nr + 1;
1458  if (Finalized())
1459  {
1460  for (int j = I[gr]; j < I[gr+1]; j++)
1461  {
1462  if (A[j] != 0.0)
1463  {
1464  blocks(bi, J[j]/nc)->I[i]++;
1465  }
1466  }
1467  }
1468  else
1469  {
1470  for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
1471  {
1472  if (n_p->Value != 0.0)
1473  {
1474  blocks(bi, n_p->Column/nc)->I[i]++;
1475  }
1476  }
1477  }
1478  }
1479 
1480  for (int j = 0; j < bc; j++)
1481  {
1482  for (int i = 0; i < br; i++)
1483  {
1484  SparseMatrix &b = *blocks(i,j);
1485  int nnz = 0, rs;
1486  for (int k = 1; k <= nr; k++)
1487  {
1488  rs = b.I[k], b.I[k] = nnz, nnz += rs;
1489  }
1490  b.J.New(nnz);
1491  b.A.New(nnz);
1492  }
1493  }
1494 
1495  for (int gr = 0; gr < height; gr++)
1496  {
1497  int bi = gr/nr, i = gr%nr + 1;
1498  if (Finalized())
1499  {
1500  for (int j = I[gr]; j < I[gr+1]; j++)
1501  {
1502  if (A[j] != 0.0)
1503  {
1504  SparseMatrix &b = *blocks(bi, J[j]/nc);
1505  b.J[b.I[i]] = J[j] % nc;
1506  b.A[b.I[i]] = A[j];
1507  b.I[i]++;
1508  }
1509  }
1510  }
1511  else
1512  {
1513  for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
1514  {
1515  if (n_p->Value != 0.0)
1516  {
1517  SparseMatrix &b = *blocks(bi, n_p->Column/nc);
1518  b.J[b.I[i]] = n_p->Column % nc;
1519  b.A[b.I[i]] = n_p->Value;
1520  b.I[i]++;
1521  }
1522  }
1523  }
1524  }
1525 }
1526 
1528 {
1529  if (height != width)
1530  {
1531  return infinity();
1532  }
1533 
1534  double symm = 0.0;
1535  if (Empty())
1536  {
1537  // return 0.0;
1538  }
1539  else if (Finalized())
1540  {
1541  for (int i = 1; i < height; i++)
1542  {
1543  for (int j = I[i]; j < I[i+1]; j++)
1544  {
1545  if (J[j] < i)
1546  {
1547  symm = std::max(symm, std::abs(A[j]-(*this)(J[j],i)));
1548  }
1549  }
1550  }
1551  }
1552  else
1553  {
1554  for (int i = 0; i < height; i++)
1555  {
1556  for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
1557  {
1558  int col = node_p->Column;
1559  if (col < i)
1560  {
1561  symm = std::max(symm, std::abs(node_p->Value-(*this)(col,i)));
1562  }
1563  }
1564  }
1565  }
1566  return symm;
1567 }
1568 
1570 {
1571  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
1572 
1573  int i, j;
1574  for (i = 1; i < height; i++)
1575  {
1576  for (j = I[i]; j < I[i+1]; j++)
1577  {
1578  if (J[j] < i)
1579  {
1580  A[j] += (*this)(J[j],i);
1581  A[j] *= 0.5;
1582  (*this)(J[j],i) = A[j];
1583  }
1584  }
1585  }
1586 }
1587 
1589 {
1590  if (Finalized())
1591  {
1592  HostReadI();
1593  return I[height];
1594  }
1595  else
1596  {
1597  int nnz = 0;
1598 
1599  for (int i = 0; i < height; i++)
1600  {
1601  for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
1602  {
1603  nnz++;
1604  }
1605  }
1606 
1607  return nnz;
1608  }
1609 }
1610 
1612 {
1613  double m = 0.0;
1614 
1615  if (A)
1616  {
1617  int nnz = I[height];
1618  for (int j = 0; j < nnz; j++)
1619  {
1620  m = std::max(m, std::abs(A[j]));
1621  }
1622  }
1623  else
1624  {
1625  for (int i = 0; i < height; i++)
1626  {
1627  for (RowNode *n_p = Rows[i]; n_p != NULL; n_p = n_p->Prev)
1628  {
1629  m = std::max(m, std::abs(n_p->Value));
1630  }
1631  }
1632  }
1633  return m;
1634 }
1635 
1636 int SparseMatrix::CountSmallElems(double tol) const
1637 {
1638  int counter = 0;
1639 
1640  if (A)
1641  {
1642  const int nz = I[height];
1643  const double *Ap = A;
1644 
1645  for (int i = 0; i < nz; i++)
1646  {
1647  counter += (std::abs(Ap[i]) <= tol);
1648  }
1649  }
1650  else
1651  {
1652  for (int i = 0; i < height; i++)
1653  {
1654  for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1655  {
1656  counter += (std::abs(aux->Value) <= tol);
1657  }
1658  }
1659  }
1660 
1661  return counter;
1662 }
1663 
1665 {
1666  if (Empty())
1667  {
1668  return 0;
1669  }
1670  else if (Finalized())
1671  {
1672  return mfem::CheckFinite(A, I[height]);
1673  }
1674  else
1675  {
1676  int counter = 0;
1677  for (int i = 0; i < height; i++)
1678  {
1679  for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1680  {
1681  counter += !IsFinite(aux->Value);
1682  }
1683  }
1684  return counter;
1685  }
1686 }
1687 
1689 {
1690  return NULL;
1691 }
1692 
1693 void SparseMatrix::EliminateRow(int row, const double sol, Vector &rhs)
1694 {
1695  RowNode *aux;
1696 
1697  MFEM_ASSERT(row < height && row >= 0,
1698  "Row " << row << " not in matrix of height " << height);
1699 
1700  MFEM_VERIFY(!Finalized(), "Matrix must NOT be finalized.");
1701 
1702  for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
1703  {
1704  rhs(aux->Column) -= sol * aux->Value;
1705  aux->Value = 0.0;
1706  }
1707 }
1708 
1710 {
1711  RowNode *aux;
1712 
1713  MFEM_ASSERT(row < height && row >= 0,
1714  "Row " << row << " not in matrix of height " << height);
1715  MFEM_ASSERT(dpolicy != DIAG_KEEP, "Diagonal policy must not be DIAG_KEEP");
1716  MFEM_ASSERT(dpolicy != DIAG_ONE || height == width,
1717  "if dpolicy == DIAG_ONE, matrix must be square, not height = "
1718  << height << ", width = " << width);
1719 
1720  if (Rows == NULL)
1721  {
1722  for (int i=I[row]; i < I[row+1]; ++i)
1723  {
1724  A[i]=0.0;
1725  }
1726  }
1727  else
1728  {
1729  for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
1730  {
1731  aux->Value = 0.0;
1732  }
1733  }
1734 
1735  if (dpolicy == DIAG_ONE)
1736  {
1737  SearchRow(row, row) = 1.;
1738  }
1739 }
1740 
1742 {
1743  MFEM_ASSERT(col < width && col >= 0,
1744  "Col " << col << " not in matrix of width " << width);
1745  MFEM_ASSERT(dpolicy != DIAG_KEEP, "Diagonal policy must not be DIAG_KEEP");
1746  MFEM_ASSERT(dpolicy != DIAG_ONE || height == width,
1747  "if dpolicy == DIAG_ONE, matrix must be square, not height = "
1748  << height << ", width = " << width);
1749 
1750  if (Rows == NULL)
1751  {
1752  const int nnz = I[height];
1753  for (int jpos = 0; jpos != nnz; ++jpos)
1754  {
1755  if (J[jpos] == col)
1756  {
1757  A[jpos] = 0.0;
1758  }
1759  }
1760  }
1761  else
1762  {
1763  for (int i = 0; i < height; i++)
1764  {
1765  for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1766  {
1767  if (aux->Column == col)
1768  {
1769  aux->Value = 0.0;
1770  }
1771  }
1772  }
1773  }
1774 
1775  if (dpolicy == DIAG_ONE)
1776  {
1777  SearchRow(col, col) = 1.0;
1778  }
1779 }
1780 
1781 void SparseMatrix::EliminateCols(const Array<int> &cols, const Vector *x,
1782  Vector *b)
1783 {
1784  if (Rows == NULL)
1785  {
1786  for (int i = 0; i < height; i++)
1787  {
1788  for (int jpos = I[i]; jpos != I[i+1]; ++jpos)
1789  {
1790  if (cols[ J[jpos]])
1791  {
1792  if (x && b)
1793  {
1794  (*b)(i) -= A[jpos] * (*x)( J[jpos] );
1795  }
1796  A[jpos] = 0.0;
1797  }
1798  }
1799  }
1800  }
1801  else
1802  {
1803  for (int i = 0; i < height; i++)
1804  {
1805  for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1806  {
1807  if (cols[aux -> Column])
1808  {
1809  if (x && b)
1810  {
1811  (*b)(i) -= aux -> Value * (*x)(aux -> Column);
1812  }
1813  aux->Value = 0.0;
1814  }
1815  }
1816  }
1817  }
1818 }
1819 
1821 {
1822  if (Rows)
1823  {
1824  RowNode *nd;
1825  for (int row = 0; row < height; row++)
1826  {
1827  for (nd = Rows[row]; nd != NULL; nd = nd->Prev)
1828  {
1829  if (col_marker[nd->Column])
1830  {
1831  Ae.Add(row, nd->Column, nd->Value);
1832  nd->Value = 0.0;
1833  }
1834  }
1835  }
1836  }
1837  else
1838  {
1839  for (int row = 0; row < height; row++)
1840  {
1841  for (int j = I[row]; j < I[row+1]; j++)
1842  {
1843  if (col_marker[J[j]])
1844  {
1845  Ae.Add(row, J[j], A[j]);
1846  A[j] = 0.0;
1847  }
1848  }
1849  }
1850  }
1851 }
1852 
1853 
1854 void SparseMatrix::EliminateRowCol(int rc, const double sol, Vector &rhs,
1855  DiagonalPolicy dpolicy)
1856 {
1857  MFEM_ASSERT(rc < height && rc >= 0,
1858  "Row " << rc << " not in matrix of height " << height);
1859  HostReadWriteI();
1860  HostReadWriteJ();
1862 
1863  if (Rows == NULL)
1864  {
1865  for (int j = I[rc]; j < I[rc+1]; j++)
1866  {
1867  const int col = J[j];
1868  if (col == rc)
1869  {
1870  switch (dpolicy)
1871  {
1872  case DIAG_KEEP:
1873  rhs(rc) = A[j] * sol;
1874  break;
1875  case DIAG_ONE:
1876  A[j] = 1.0;
1877  rhs(rc) = sol;
1878  break;
1879  case DIAG_ZERO:
1880  A[j] = 0.;
1881  rhs(rc) = 0.;
1882  break;
1883  default:
1884  mfem_error("SparseMatrix::EliminateRowCol () #2");
1885  break;
1886  }
1887  }
1888  else
1889  {
1890  A[j] = 0.0;
1891  for (int k = I[col]; 1; k++)
1892  {
1893  if (k == I[col+1])
1894  {
1895  mfem_error("SparseMatrix::EliminateRowCol () #3");
1896  }
1897  else if (J[k] == rc)
1898  {
1899  rhs(col) -= sol * A[k];
1900  A[k] = 0.0;
1901  break;
1902  }
1903  }
1904  }
1905  }
1906  }
1907  else
1908  {
1909  for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
1910  {
1911  const int col = aux->Column;
1912  if (col == rc)
1913  {
1914  switch (dpolicy)
1915  {
1916  case DIAG_KEEP:
1917  rhs(rc) = aux->Value * sol;
1918  break;
1919  case DIAG_ONE:
1920  aux->Value = 1.0;
1921  rhs(rc) = sol;
1922  break;
1923  case DIAG_ZERO:
1924  aux->Value = 0.;
1925  rhs(rc) = 0.;
1926  break;
1927  default:
1928  mfem_error("SparseMatrix::EliminateRowCol () #4");
1929  break;
1930  }
1931  }
1932  else
1933  {
1934  aux->Value = 0.0;
1935  for (RowNode *node = Rows[col]; 1; node = node->Prev)
1936  {
1937  if (node == NULL)
1938  {
1939  mfem_error("SparseMatrix::EliminateRowCol () #5");
1940  }
1941  else if (node->Column == rc)
1942  {
1943  rhs(col) -= sol * node->Value;
1944  node->Value = 0.0;
1945  break;
1946  }
1947  }
1948  }
1949  }
1950  }
1951 }
1952 
1954  DenseMatrix &rhs,
1955  DiagonalPolicy dpolicy)
1956 {
1957  MFEM_ASSERT(rc < height && rc >= 0,
1958  "Row " << rc << " not in matrix of height " << height);
1959  MFEM_ASSERT(sol.Size() == rhs.Width(), "solution size (" << sol.Size()
1960  << ") must match rhs width (" << rhs.Width() << ")");
1961 
1962  const int num_rhs = rhs.Width();
1963  if (Rows == NULL)
1964  {
1965  for (int j = I[rc]; j < I[rc+1]; j++)
1966  {
1967  const int col = J[j];
1968  if (col == rc)
1969  {
1970  switch (dpolicy)
1971  {
1972  case DIAG_KEEP:
1973  for (int r = 0; r < num_rhs; r++)
1974  {
1975  rhs(rc,r) = A[j] * sol(r);
1976  }
1977  break;
1978  case DIAG_ONE:
1979  A[j] = 1.0;
1980  for (int r = 0; r < num_rhs; r++)
1981  {
1982  rhs(rc,r) = sol(r);
1983  }
1984  break;
1985  case DIAG_ZERO:
1986  A[j] = 0.;
1987  for (int r = 0; r < num_rhs; r++)
1988  {
1989  rhs(rc,r) = 0.;
1990  }
1991  break;
1992  default:
1993  mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #3");
1994  break;
1995  }
1996  }
1997  else
1998  {
1999  A[j] = 0.0;
2000  for (int k = I[col]; 1; k++)
2001  {
2002  if (k == I[col+1])
2003  {
2004  mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #4");
2005  }
2006  else if (J[k] == rc)
2007  {
2008  for (int r = 0; r < num_rhs; r++)
2009  {
2010  rhs(col,r) -= sol(r) * A[k];
2011  }
2012  A[k] = 0.0;
2013  break;
2014  }
2015  }
2016  }
2017  }
2018  }
2019  else
2020  {
2021  for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
2022  {
2023  const int col = aux->Column;
2024  if (col == rc)
2025  {
2026  switch (dpolicy)
2027  {
2028  case DIAG_KEEP:
2029  for (int r = 0; r < num_rhs; r++)
2030  {
2031  rhs(rc,r) = aux->Value * sol(r);
2032  }
2033  break;
2034  case DIAG_ONE:
2035  aux->Value = 1.0;
2036  for (int r = 0; r < num_rhs; r++)
2037  {
2038  rhs(rc,r) = sol(r);
2039  }
2040  break;
2041  case DIAG_ZERO:
2042  aux->Value = 0.;
2043  for (int r = 0; r < num_rhs; r++)
2044  {
2045  rhs(rc,r) = 0.;
2046  }
2047  break;
2048  default:
2049  mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #5");
2050  break;
2051  }
2052  }
2053  else
2054  {
2055  aux->Value = 0.0;
2056  for (RowNode *node = Rows[col]; 1; node = node->Prev)
2057  {
2058  if (node == NULL)
2059  {
2060  mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #6");
2061  }
2062  else if (node->Column == rc)
2063  {
2064  for (int r = 0; r < num_rhs; r++)
2065  {
2066  rhs(col,r) -= sol(r) * node->Value;
2067  }
2068  node->Value = 0.0;
2069  break;
2070  }
2071  }
2072  }
2073  }
2074  }
2075 }
2076 
2078 {
2079  MFEM_ASSERT(rc < height && rc >= 0,
2080  "Row " << rc << " not in matrix of height " << height);
2081 
2082  if (Rows == NULL)
2083  {
2084  const auto &II = this->I; // only use const access for I
2085  const auto &JJ = this->J; // only use const access for J
2086  for (int j = II[rc]; j < II[rc+1]; j++)
2087  {
2088  const int col = JJ[j];
2089  if (col == rc)
2090  {
2091  if (dpolicy == DIAG_ONE)
2092  {
2093  A[j] = 1.0;
2094  }
2095  else if (dpolicy == DIAG_ZERO)
2096  {
2097  A[j] = 0.0;
2098  }
2099  }
2100  else
2101  {
2102  A[j] = 0.0;
2103  for (int k = II[col]; 1; k++)
2104  {
2105  if (k == II[col+1])
2106  {
2107  mfem_error("SparseMatrix::EliminateRowCol() #2");
2108  }
2109  else if (JJ[k] == rc)
2110  {
2111  A[k] = 0.0;
2112  break;
2113  }
2114  }
2115  }
2116  }
2117  }
2118  else
2119  {
2120  RowNode *aux, *node;
2121 
2122  for (aux = Rows[rc]; aux != NULL; aux = aux->Prev)
2123  {
2124  const int col = aux->Column;
2125  if (col == rc)
2126  {
2127  if (dpolicy == DIAG_ONE)
2128  {
2129  aux->Value = 1.0;
2130  }
2131  else if (dpolicy == DIAG_ZERO)
2132  {
2133  aux->Value = 0.;
2134  }
2135  }
2136  else
2137  {
2138  aux->Value = 0.0;
2139  for (node = Rows[col]; 1; node = node->Prev)
2140  {
2141  if (node == NULL)
2142  {
2143  mfem_error("SparseMatrix::EliminateRowCol() #3");
2144  }
2145  else if (node->Column == rc)
2146  {
2147  node->Value = 0.0;
2148  break;
2149  }
2150  }
2151  }
2152  }
2153  }
2154 }
2155 
2156 // This is almost identical to EliminateRowCol(int, int), except for
2157 // the A[j] = value; and aux->Value = value; lines.
2158 void SparseMatrix::EliminateRowColDiag(int rc, double value)
2159 {
2160  MFEM_ASSERT(rc < height && rc >= 0,
2161  "Row " << rc << " not in matrix of height " << height);
2162 
2163  if (Rows == NULL)
2164  {
2165  for (int j = I[rc]; j < I[rc+1]; j++)
2166  {
2167  const int col = J[j];
2168  if (col == rc)
2169  {
2170  A[j] = value;
2171  }
2172  else
2173  {
2174  A[j] = 0.0;
2175  for (int k = I[col]; 1; k++)
2176  {
2177  if (k == I[col+1])
2178  {
2179  mfem_error("SparseMatrix::EliminateRowCol() #2");
2180  }
2181  else if (J[k] == rc)
2182  {
2183  A[k] = 0.0;
2184  break;
2185  }
2186  }
2187  }
2188  }
2189  }
2190  else
2191  {
2192  RowNode *aux, *node;
2193 
2194  for (aux = Rows[rc]; aux != NULL; aux = aux->Prev)
2195  {
2196  const int col = aux->Column;
2197  if (col == rc)
2198  {
2199  aux->Value = value;
2200  }
2201  else
2202  {
2203  aux->Value = 0.0;
2204  for (node = Rows[col]; 1; node = node->Prev)
2205  {
2206  if (node == NULL)
2207  {
2208  mfem_error("SparseMatrix::EliminateRowCol() #3");
2209  }
2210  else if (node->Column == rc)
2211  {
2212  node->Value = 0.0;
2213  break;
2214  }
2215  }
2216  }
2217  }
2218  }
2219 }
2220 
2222  DiagonalPolicy dpolicy)
2223 {
2224  if (Rows)
2225  {
2226  RowNode *nd, *nd2;
2227  for (nd = Rows[rc]; nd != NULL; nd = nd->Prev)
2228  {
2229  const int col = nd->Column;
2230  if (col == rc)
2231  {
2232  switch (dpolicy)
2233  {
2234  case DIAG_ONE:
2235  Ae.Add(rc, rc, nd->Value - 1.0);
2236  nd->Value = 1.0;
2237  break;
2238  case DIAG_ZERO:
2239  Ae.Add(rc, rc, nd->Value);
2240  nd->Value = 0.;
2241  break;
2242  case DIAG_KEEP:
2243  break;
2244  default:
2245  mfem_error("SparseMatrix::EliminateRowCol #1");
2246  break;
2247  }
2248  }
2249  else
2250  {
2251  Ae.Add(rc, col, nd->Value);
2252  nd->Value = 0.0;
2253  for (nd2 = Rows[col]; 1; nd2 = nd2->Prev)
2254  {
2255  if (nd2 == NULL)
2256  {
2257  mfem_error("SparseMatrix::EliminateRowCol #2");
2258  }
2259  else if (nd2->Column == rc)
2260  {
2261  Ae.Add(col, rc, nd2->Value);
2262  nd2->Value = 0.0;
2263  break;
2264  }
2265  }
2266  }
2267  }
2268  }
2269  else
2270  {
2271  for (int j = I[rc]; j < I[rc+1]; j++)
2272  {
2273  const int col = J[j];
2274  if (col == rc)
2275  {
2276  switch (dpolicy)
2277  {
2278  case DIAG_ONE:
2279  Ae.Add(rc, rc, A[j] - 1.0);
2280  A[j] = 1.0;
2281  break;
2282  case DIAG_ZERO:
2283  Ae.Add(rc, rc, A[j]);
2284  A[j] = 0.;
2285  break;
2286  case DIAG_KEEP:
2287  break;
2288  default:
2289  mfem_error("SparseMatrix::EliminateRowCol #3");
2290  break;
2291  }
2292  }
2293  else
2294  {
2295  Ae.Add(rc, col, A[j]);
2296  A[j] = 0.0;
2297  for (int k = I[col]; true; k++)
2298  {
2299  if (k == I[col+1])
2300  {
2301  mfem_error("SparseMatrix::EliminateRowCol #4");
2302  }
2303  else if (J[k] == rc)
2304  {
2305  Ae.Add(col, rc, A[k]);
2306  A[k] = 0.0;
2307  break;
2308  }
2309  }
2310  }
2311  }
2312  }
2313 }
2314 
2316  DiagonalPolicy diag_policy)
2317 {
2318  const int n_ess_dofs = ess_dofs.Size();
2319  const auto ess_dofs_d = ess_dofs.Read();
2320  const auto dI = ReadI();
2321  const auto dJ = ReadJ();
2322  auto dA = ReadWriteData();
2323 
2324  MFEM_FORALL(i, n_ess_dofs,
2325  {
2326  const int idof = ess_dofs_d[i];
2327  for (int j=dI[idof]; j<dI[idof+1]; ++j)
2328  {
2329  const int jdof = dJ[j];
2330  if (jdof != idof)
2331  {
2332  dA[j] = 0.0;
2333  for (int k=dI[jdof]; k<dI[jdof+1]; ++k)
2334  {
2335  if (dJ[k] == idof)
2336  {
2337  dA[k] = 0.0;
2338  break;
2339  }
2340  }
2341  }
2342  else
2343  {
2344  if (diag_policy == DiagonalPolicy::DIAG_ONE)
2345  {
2346  dA[j] = 1.0;
2347  }
2348  else if (diag_policy == DiagonalPolicy::DIAG_ZERO)
2349  {
2350  dA[j] = 0.0;
2351  }
2352  // else (diag_policy == DiagonalPolicy::DIAG_KEEP)
2353  }
2354  }
2355  });
2356 }
2357 
2359 {
2360  for (int i = 0; i < height; i++)
2361  {
2362  if (I[i+1] == I[i]+1 && fabs(A[I[i]]) < 1e-16)
2363  {
2364  A[I[i]] = 1.0;
2365  }
2366  }
2367 }
2368 
2369 void SparseMatrix::EliminateZeroRows(const double threshold)
2370 {
2371  for (int i = 0; i < height; i++)
2372  {
2373  double zero = 0.0;
2374  for (int j = I[i]; j < I[i+1]; j++)
2375  {
2376  zero += fabs(A[j]);
2377  }
2378  if (zero <= threshold)
2379  {
2380  for (int j = I[i]; j < I[i+1]; j++)
2381  {
2382  A[j] = (J[j] == i) ? 1.0 : 0.0;
2383  }
2384  }
2385  }
2386 }
2387 
2389 {
2390  if (!Finalized())
2391  {
2392  double *yp = y.GetData();
2393  const double *xp = x.GetData();
2394  RowNode *diag_p, *n_p, **R = Rows;
2395 
2396  const int s = height;
2397  for (int i = 0; i < s; i++)
2398  {
2399  double sum = 0.0;
2400  diag_p = NULL;
2401  for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
2402  {
2403  const int c = n_p->Column;
2404  if (c == i)
2405  {
2406  diag_p = n_p;
2407  }
2408  else
2409  {
2410  sum += n_p->Value * yp[c];
2411  }
2412  }
2413 
2414  if (diag_p != NULL && diag_p->Value != 0.0)
2415  {
2416  yp[i] = (xp[i] - sum) / diag_p->Value;
2417  }
2418  else if (xp[i] == sum)
2419  {
2420  yp[i] = sum;
2421  }
2422  else
2423  {
2424  mfem_error("SparseMatrix::Gauss_Seidel_forw()");
2425  }
2426  }
2427  }
2428  else
2429  {
2430  const int s = height;
2431  const int nnz = J.Capacity();
2432  const int *Ip = HostRead(I, s+1);
2433  const int *Jp = HostRead(J, nnz);
2434  const double *Ap = HostRead(A, nnz);
2435  double *yp = y.HostReadWrite();
2436  const double *xp = x.HostRead();
2437 
2438  for (int i = 0, j = Ip[0]; i < s; i++)
2439  {
2440  const int end = Ip[i+1];
2441  double sum = 0.0;
2442  int d = -1;
2443  for ( ; j < end; j++)
2444  {
2445  const int c = Jp[j];
2446  if (c == i)
2447  {
2448  d = j;
2449  }
2450  else
2451  {
2452  sum += Ap[j] * yp[c];
2453  }
2454  }
2455 
2456  if (d >= 0 && Ap[d] != 0.0)
2457  {
2458  yp[i] = (xp[i] - sum) / Ap[d];
2459  }
2460  else if (xp[i] == sum)
2461  {
2462  yp[i] = sum;
2463  }
2464  else
2465  {
2466  mfem_error("SparseMatrix::Gauss_Seidel_forw(...) #2");
2467  }
2468  }
2469  }
2470 }
2471 
2473 {
2474  if (!Finalized())
2475  {
2476  double *yp = y.GetData();
2477  const double *xp = x.GetData();
2478  RowNode *diag_p, *n_p, **R = Rows;
2479 
2480  for (int i = height-1; i >= 0; i--)
2481  {
2482  double sum = 0.;
2483  diag_p = NULL;
2484  for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
2485  {
2486  const int c = n_p->Column;
2487  if (c == i)
2488  {
2489  diag_p = n_p;
2490  }
2491  else
2492  {
2493  sum += n_p->Value * yp[c];
2494  }
2495  }
2496 
2497  if (diag_p != NULL && diag_p->Value != 0.0)
2498  {
2499  yp[i] = (xp[i] - sum) / diag_p->Value;
2500  }
2501  else if (xp[i] == sum)
2502  {
2503  yp[i] = sum;
2504  }
2505  else
2506  {
2507  mfem_error("SparseMatrix::Gauss_Seidel_back()");
2508  }
2509  }
2510  }
2511  else
2512  {
2513  const int s = height;
2514  const int nnz = J.Capacity();
2515  const int *Ip = HostRead(I, s+1);
2516  const int *Jp = HostRead(J, nnz);
2517  const double *Ap = HostRead(A, nnz);
2518  double *yp = y.HostReadWrite();
2519  const double *xp = x.HostRead();
2520 
2521  for (int i = s-1, j = Ip[s]-1; i >= 0; i--)
2522  {
2523  const int beg = Ip[i];
2524  double sum = 0.;
2525  int d = -1;
2526  for ( ; j >= beg; j--)
2527  {
2528  const int c = Jp[j];
2529  if (c == i)
2530  {
2531  d = j;
2532  }
2533  else
2534  {
2535  sum += Ap[j] * yp[c];
2536  }
2537  }
2538 
2539  if (d >= 0 && Ap[d] != 0.0)
2540  {
2541  yp[i] = (xp[i] - sum) / Ap[d];
2542  }
2543  else if (xp[i] == sum)
2544  {
2545  yp[i] = sum;
2546  }
2547  else
2548  {
2549  mfem_error("SparseMatrix::Gauss_Seidel_back(...) #2");
2550  }
2551  }
2552  }
2553 }
2554 
2556 {
2557  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2558 
2559  double sc = 1.0;
2560  for (int i = 0; i < height; i++)
2561  {
2562  int d = -1;
2563  double norm = 0.0;
2564  for (int j = I[i]; j < I[i+1]; j++)
2565  {
2566  if (J[j] == i)
2567  {
2568  d = j;
2569  }
2570  norm += fabs(A[j]);
2571  }
2572  if (d >= 0 && A[d] != 0.0)
2573  {
2574  double a = 1.8 * fabs(A[d]) / norm;
2575  if (a < sc)
2576  {
2577  sc = a;
2578  }
2579  }
2580  else
2581  {
2582  mfem_error("SparseMatrix::GetJacobiScaling() #2");
2583  }
2584  }
2585  return sc;
2586 }
2587 
2588 void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
2589  double sc, bool use_abs_diag) const
2590 {
2591  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2592 
2593  for (int i = 0; i < height; i++)
2594  {
2595  int d = -1;
2596  double sum = b(i);
2597  for (int j = I[i]; j < I[i+1]; j++)
2598  {
2599  if (J[j] == i)
2600  {
2601  d = j;
2602  }
2603  else
2604  {
2605  sum -= A[j] * x0(J[j]);
2606  }
2607  }
2608  if (d >= 0 && A[d] != 0.0)
2609  {
2610  const double diag = (use_abs_diag) ? fabs(A[d]) : A[d];
2611  x1(i) = sc * (sum / diag) + (1.0 - sc) * x0(i);
2612  }
2613  else
2614  {
2615  mfem_error("SparseMatrix::Jacobi(...) #2");
2616  }
2617  }
2618 }
2619 
2621  double sc, bool use_abs_diag) const
2622 {
2623  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2624 
2625  const int H = height;
2626  const int nnz = J.Capacity();
2627  const bool use_dev = b.UseDevice() || x.UseDevice();
2628 
2629  const auto Ap = Read(A, nnz, use_dev);
2630  const auto Ip = Read(I, height+1, use_dev);
2631  const auto Jp = Read(J, nnz, use_dev);
2632 
2633  const auto bp = b.Read(use_dev);
2634  auto xp = x.Write(use_dev);
2635 
2636  MFEM_FORALL_SWITCH(use_dev, i, H,
2637  {
2638  const int end = Ip[i+1];
2639  for (int j = Ip[i]; true; j++)
2640  {
2641  if (j == end)
2642  {
2643  MFEM_ABORT_KERNEL("Diagonal not found in SparseMatrix::DiagScale");
2644  }
2645  if (Jp[j] == i)
2646  {
2647  const double diag = (use_abs_diag) ? fabs(Ap[j]) : Ap[j];
2648  if (diag == 0.0)
2649  {
2650  MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
2651  }
2652  xp[i] = sc * bp[i] / diag;
2653  break;
2654  }
2655  }
2656  });
2657 }
2658 
2659 template <bool useFabs>
2660 static void JacobiDispatch(const Vector &b, const Vector &x0, Vector &x1,
2661  const Memory<int> &I, const Memory<int> &J,
2662  const Memory<double> &A, const int height,
2663  const double sc)
2664 {
2665  const bool useDevice = b.UseDevice() || x0.UseDevice() || x1.UseDevice();
2666 
2667  const auto bp = b.Read(useDevice);
2668  const auto x0p = x0.Read(useDevice);
2669  auto x1p = x1.Write(useDevice);
2670 
2671  const auto Ip = Read(I, height+1, useDevice);
2672  const auto Jp = Read(J, J.Capacity(), useDevice);
2673  const auto Ap = Read(A, J.Capacity(), useDevice);
2674 
2675  MFEM_FORALL_SWITCH(useDevice, i, height,
2676  {
2677  double resi = bp[i], norm = 0.0;
2678  for (int j = Ip[i]; j < Ip[i+1]; j++)
2679  {
2680  resi -= Ap[j] * x0p[Jp[j]];
2681  if (useFabs)
2682  {
2683  norm += fabs(Ap[j]);
2684  }
2685  else
2686  {
2687  norm += Ap[j];
2688  }
2689  }
2690  if (norm > 0.0)
2691  {
2692  x1p[i] = x0p[i] + sc * resi / norm;
2693  }
2694  else
2695  {
2696  if (useFabs)
2697  {
2698  MFEM_ABORT_KERNEL("L1 norm of row is zero.");
2699  }
2700  else
2701  {
2702  MFEM_ABORT_KERNEL("sum of row is zero.");
2703  }
2704  }
2705  });
2706 }
2707 
2708 void SparseMatrix::Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
2709  double sc) const
2710 {
2711  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2712  JacobiDispatch<true>(b,x0,x1,I,J,A,height,sc);
2713 }
2714 
2715 void SparseMatrix::Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
2716  double sc) const
2717 {
2718  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2719  JacobiDispatch<false>(b,x0,x1,I,J,A,height,sc);
2720 }
2721 
2722 void SparseMatrix::AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
2723  const DenseMatrix &subm, int skip_zeros)
2724 {
2725  int i, j, gi, gj, s, t;
2726  double a;
2727 
2728  if (Finalized())
2729  {
2730  HostReadI();
2731  HostReadJ();
2733  }
2734 
2735  for (i = 0; i < rows.Size(); i++)
2736  {
2737  if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2738  else { s = 1; }
2739  MFEM_ASSERT(gi < height,
2740  "Trying to insert a row " << gi << " outside the matrix height "
2741  << height);
2742  SetColPtr(gi);
2743  for (j = 0; j < cols.Size(); j++)
2744  {
2745  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2746  else { t = s; }
2747  MFEM_ASSERT(gj < width,
2748  "Trying to insert a column " << gj << " outside the matrix width "
2749  << width);
2750  a = subm(i, j);
2751  if (skip_zeros && a == 0.0)
2752  {
2753  // Skip assembly of zero elements if either:
2754  // (i) user specified to skip zeros regardless of symmetry, or
2755  // (ii) symmetry is not broken.
2756  if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2757  {
2758  continue;
2759  }
2760  }
2761  if (t < 0) { a = -a; }
2762  _Add_(gj, a);
2763  }
2764  ClearColPtr();
2765  }
2766 }
2767 
2768 void SparseMatrix::Set(const int i, const int j, const double val)
2769 {
2770  double a = val;
2771  int gi, gj, s, t;
2772 
2773  if ((gi=i) < 0) { gi = -1-gi, s = -1; }
2774  else { s = 1; }
2775  MFEM_ASSERT(gi < height,
2776  "Trying to set a row " << gi << " outside the matrix height "
2777  << height);
2778  if ((gj=j) < 0) { gj = -1-gj, t = -s; }
2779  else { t = s; }
2780  MFEM_ASSERT(gj < width,
2781  "Trying to set a column " << gj << " outside the matrix width "
2782  << width);
2783  if (t < 0) { a = -a; }
2784  _Set_(gi, gj, a);
2785 }
2786 
2787 void SparseMatrix::Add(const int i, const int j, const double val)
2788 {
2789  int gi, gj, s, t;
2790  double a = val;
2791 
2792  if ((gi=i) < 0) { gi = -1-gi, s = -1; }
2793  else { s = 1; }
2794  MFEM_ASSERT(gi < height,
2795  "Trying to insert a row " << gi << " outside the matrix height "
2796  << height);
2797  if ((gj=j) < 0) { gj = -1-gj, t = -s; }
2798  else { t = s; }
2799  MFEM_ASSERT(gj < width,
2800  "Trying to insert a column " << gj << " outside the matrix width "
2801  << width);
2802  if (t < 0) { a = -a; }
2803  _Add_(gi, gj, a);
2804 }
2805 
2806 void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
2807  const DenseMatrix &subm, int skip_zeros)
2808 {
2809  int i, j, gi, gj, s, t;
2810  double a;
2811 
2812  for (i = 0; i < rows.Size(); i++)
2813  {
2814  if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2815  else { s = 1; }
2816  MFEM_ASSERT(gi < height,
2817  "Trying to set a row " << gi << " outside the matrix height "
2818  << height);
2819  SetColPtr(gi);
2820  for (j = 0; j < cols.Size(); j++)
2821  {
2822  a = subm(i, j);
2823  if (skip_zeros && a == 0.0)
2824  {
2825  // Skip assembly of zero elements if either:
2826  // (i) user specified to skip zeros regardless of symmetry, or
2827  // (ii) symmetry is not broken.
2828  if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2829  {
2830  continue;
2831  }
2832  }
2833  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2834  else { t = s; }
2835  MFEM_ASSERT(gj < width,
2836  "Trying to set a column " << gj << " outside the matrix width "
2837  << width);
2838  if (t < 0) { a = -a; }
2839  _Set_(gj, a);
2840  }
2841  ClearColPtr();
2842  }
2843 }
2844 
2846  const Array<int> &cols,
2847  const DenseMatrix &subm,
2848  int skip_zeros)
2849 {
2850  int i, j, gi, gj, s, t;
2851  double a;
2852 
2853  for (i = 0; i < rows.Size(); i++)
2854  {
2855  if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2856  else { s = 1; }
2857  MFEM_ASSERT(gi < height,
2858  "Trying to set a row " << gi << " outside the matrix height "
2859  << height);
2860  SetColPtr(gi);
2861  for (j = 0; j < cols.Size(); j++)
2862  {
2863  a = subm(j, i);
2864  if (skip_zeros && a == 0.0)
2865  {
2866  // Skip assembly of zero elements if either:
2867  // (i) user specified to skip zeros regardless of symmetry, or
2868  // (ii) symmetry is not broken.
2869  if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2870  {
2871  continue;
2872  }
2873  }
2874  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2875  else { t = s; }
2876  MFEM_ASSERT(gj < width,
2877  "Trying to set a column " << gj << " outside the matrix width "
2878  << width);
2879  if (t < 0) { a = -a; }
2880  _Set_(gj, a);
2881  }
2882  ClearColPtr();
2883  }
2884 }
2885 
2886 void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
2887  DenseMatrix &subm) const
2888 {
2889  int i, j, gi, gj, s, t;
2890  double a;
2891 
2892  for (i = 0; i < rows.Size(); i++)
2893  {
2894  if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2895  else { s = 1; }
2896  MFEM_ASSERT(gi < height,
2897  "Trying to read a row " << gi << " outside the matrix height "
2898  << height);
2899  SetColPtr(gi);
2900  for (j = 0; j < cols.Size(); j++)
2901  {
2902  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2903  else { t = s; }
2904  MFEM_ASSERT(gj < width,
2905  "Trying to read a column " << gj << " outside the matrix width "
2906  << width);
2907  a = _Get_(gj);
2908  subm(i, j) = (t < 0) ? (-a) : (a);
2909  }
2910  ClearColPtr();
2911  }
2912 }
2913 
2914 bool SparseMatrix::RowIsEmpty(const int row) const
2915 {
2916  int gi;
2917 
2918  if ((gi=row) < 0)
2919  {
2920  gi = -1-gi;
2921  }
2922  MFEM_ASSERT(gi < height,
2923  "Trying to query a row " << gi << " outside the matrix height "
2924  << height);
2925  if (Rows)
2926  {
2927  return (Rows[gi] == NULL);
2928  }
2929  else
2930  {
2931  return (I[gi] == I[gi+1]);
2932  }
2933 }
2934 
2935 int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
2936 {
2937  RowNode *n;
2938  int j, gi;
2939 
2940  if ((gi=row) < 0) { gi = -1-gi; }
2941  MFEM_ASSERT(gi < height,
2942  "Trying to read a row " << gi << " outside the matrix height "
2943  << height);
2944  if (Rows)
2945  {
2946  for (n = Rows[gi], j = 0; n; n = n->Prev)
2947  {
2948  j++;
2949  }
2950  cols.SetSize(j);
2951  srow.SetSize(j);
2952  for (n = Rows[gi], j = 0; n; n = n->Prev, j++)
2953  {
2954  cols[j] = n->Column;
2955  srow(j) = n->Value;
2956  }
2957  if (row < 0)
2958  {
2959  srow.Neg();
2960  }
2961 
2962  return 0;
2963  }
2964  else
2965  {
2966  j = I[gi];
2967  cols.MakeRef(const_cast<int*>((const int*)J) + j, I[gi+1]-j);
2968  srow.NewDataAndSize(
2969  const_cast<double*>((const double*)A) + j, cols.Size());
2970  MFEM_ASSERT(row >= 0, "Row not valid: " << row << ", height: " << height);
2971  return 1;
2972  }
2973 }
2974 
2975 void SparseMatrix::SetRow(const int row, const Array<int> &cols,
2976  const Vector &srow)
2977 {
2978  int gi, gj, s, t;
2979  double a;
2980 
2981  if ((gi=row) < 0) { gi = -1-gi, s = -1; }
2982  else { s = 1; }
2983  MFEM_ASSERT(gi < height,
2984  "Trying to set a row " << gi << " outside the matrix height "
2985  << height);
2986 
2987  if (!Finalized())
2988  {
2989  SetColPtr(gi);
2990  for (int j = 0; j < cols.Size(); j++)
2991  {
2992  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2993  else { t = s; }
2994  MFEM_ASSERT(gj < width,
2995  "Trying to set a column " << gj << " outside the matrix"
2996  " width " << width);
2997  a = srow(j);
2998  if (t < 0) { a = -a; }
2999  _Set_(gj, a);
3000  }
3001  ClearColPtr();
3002  }
3003  else
3004  {
3005  MFEM_ASSERT(cols.Size() == RowSize(gi), "");
3006  MFEM_ASSERT(cols.Size() == srow.Size(), "");
3007 
3008  for (int i = I[gi], j = 0; j < cols.Size(); j++, i++)
3009  {
3010  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
3011  else { t = s; }
3012  MFEM_ASSERT(gj < width,
3013  "Trying to set a column " << gj << " outside the matrix"
3014  " width " << width);
3015 
3016  J[i] = gj;
3017  A[i] = srow[j] * t;
3018  }
3019  }
3020 }
3021 
3022 void SparseMatrix::AddRow(const int row, const Array<int> &cols,
3023  const Vector &srow)
3024 {
3025  int j, gi, gj, s, t;
3026  double a;
3027 
3028  MFEM_VERIFY(!Finalized(), "Matrix must NOT be finalized.");
3029 
3030  if ((gi=row) < 0) { gi = -1-gi, s = -1; }
3031  else { s = 1; }
3032  MFEM_ASSERT(gi < height,
3033  "Trying to insert a row " << gi << " outside the matrix height "
3034  << height);
3035  SetColPtr(gi);
3036  for (j = 0; j < cols.Size(); j++)
3037  {
3038  if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
3039  else { t = s; }
3040  MFEM_ASSERT(gj < width,
3041  "Trying to insert a column " << gj << " outside the matrix width "
3042  << width);
3043  a = srow(j);
3044  if (a == 0.0)
3045  {
3046  continue;
3047  }
3048  if (t < 0) { a = -a; }
3049  _Add_(gj, a);
3050  }
3051  ClearColPtr();
3052 }
3053 
3054 void SparseMatrix::ScaleRow(const int row, const double scale)
3055 {
3056  int i;
3057 
3058  if ((i=row) < 0)
3059  {
3060  i = -1-i;
3061  }
3062  if (Rows != NULL)
3063  {
3064  RowNode *aux;
3065 
3066  for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
3067  {
3068  aux -> Value *= scale;
3069  }
3070  }
3071  else
3072  {
3073  int j, end = I[i+1];
3074 
3075  for (j = I[i]; j < end; j++)
3076  {
3077  A[j] *= scale;
3078  }
3079  }
3080 }
3081 
3083 {
3084  double scale;
3085  if (Rows != NULL)
3086  {
3087  RowNode *aux;
3088  for (int i=0; i < height; ++i)
3089  {
3090  scale = sl(i);
3091  for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
3092  {
3093  aux -> Value *= scale;
3094  }
3095  }
3096  }
3097  else
3098  {
3099  int j, end;
3100 
3101  for (int i=0; i < height; ++i)
3102  {
3103  end = I[i+1];
3104  scale = sl(i);
3105  for (j = I[i]; j < end; j++)
3106  {
3107  A[j] *= scale;
3108  }
3109  }
3110  }
3111 }
3112 
3114 {
3115  if (Rows != NULL)
3116  {
3117  RowNode *aux;
3118  for (int i=0; i < height; ++i)
3119  {
3120  for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
3121  {
3122  aux -> Value *= sr(aux->Column);
3123  }
3124  }
3125  }
3126  else
3127  {
3128  int j, end;
3129 
3130  for (int i=0; i < height; ++i)
3131  {
3132  end = I[i+1];
3133  for (j = I[i]; j < end; j++)
3134  {
3135  A[j] *= sr(J[j]);
3136  }
3137  }
3138  }
3139 }
3140 
3142 {
3143  MFEM_ASSERT(height == B.height && width == B.width,
3144  "Mismatch of this matrix size and rhs. This height = "
3145  << height << ", width = " << width << ", B.height = "
3146  << B.height << ", B.width = " << B.width);
3147 
3148  for (int i = 0; i < height; i++)
3149  {
3150  SetColPtr(i);
3151  if (B.Rows)
3152  {
3153  for (RowNode *aux = B.Rows[i]; aux != NULL; aux = aux->Prev)
3154  {
3155  _Add_(aux->Column, aux->Value);
3156  }
3157  }
3158  else
3159  {
3160  for (int j = B.I[i]; j < B.I[i+1]; j++)
3161  {
3162  _Add_(B.J[j], B.A[j]);
3163  }
3164  }
3165  ClearColPtr();
3166  }
3167 
3168  return (*this);
3169 }
3170 
3171 void SparseMatrix::Add(const double a, const SparseMatrix &B)
3172 {
3173  for (int i = 0; i < height; i++)
3174  {
3175  B.SetColPtr(i);
3176  if (Rows)
3177  {
3178  for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
3179  {
3180  np->Value += a * B._Get_(np->Column);
3181  }
3182  }
3183  else
3184  {
3185  for (int j = I[i]; j < I[i+1]; j++)
3186  {
3187  A[j] += a * B._Get_(J[j]);
3188  }
3189  }
3190  B.ClearColPtr();
3191  }
3192 }
3193 
3195 {
3196  if (Rows == NULL)
3197  {
3198  const int nnz = J.Capacity();
3199  double *h_A = HostWrite(A, nnz);
3200  for (int i = 0; i < nnz; i++)
3201  {
3202  h_A[i] = a;
3203  }
3204  }
3205  else
3206  {
3207  for (int i = 0; i < height; i++)
3208  {
3209  for (RowNode *node_p = Rows[i]; node_p != NULL;
3210  node_p = node_p -> Prev)
3211  {
3212  node_p -> Value = a;
3213  }
3214  }
3215  }
3216 
3217  return (*this);
3218 }
3219 
3221 {
3222  if (Rows == NULL)
3223  {
3224  for (int i = 0, nnz = I[height]; i < nnz; i++)
3225  {
3226  A[i] *= a;
3227  }
3228  }
3229  else
3230  {
3231  for (int i = 0; i < height; i++)
3232  {
3233  for (RowNode *node_p = Rows[i]; node_p != NULL;
3234  node_p = node_p -> Prev)
3235  {
3236  node_p -> Value *= a;
3237  }
3238  }
3239  }
3240 
3241  return (*this);
3242 }
3243 
3244 void SparseMatrix::Print(std::ostream & os, int width_) const
3245 {
3246  int i, j;
3247 
3248  if (A == NULL)
3249  {
3250  RowNode *nd;
3251  for (i = 0; i < height; i++)
3252  {
3253  os << "[row " << i << "]\n";
3254  for (nd = Rows[i], j = 0; nd != NULL; nd = nd->Prev, j++)
3255  {
3256  os << " (" << nd->Column << "," << nd->Value << ")";
3257  if ( !((j+1) % width_) )
3258  {
3259  os << '\n';
3260  }
3261  }
3262  if (j % width_)
3263  {
3264  os << '\n';
3265  }
3266  }
3267  return;
3268  }
3269 
3270  // HostRead forces synchronization
3271  HostReadI();
3272  HostReadJ();
3273  HostReadData();
3274  for (i = 0; i < height; i++)
3275  {
3276  os << "[row " << i << "]\n";
3277  for (j = I[i]; j < I[i+1]; j++)
3278  {
3279  os << " (" << J[j] << "," << A[j] << ")";
3280  if ( !((j+1-I[i]) % width_) )
3281  {
3282  os << '\n';
3283  }
3284  }
3285  if ((j-I[i]) % width_)
3286  {
3287  os << '\n';
3288  }
3289  }
3290 }
3291 
3292 void SparseMatrix::PrintMatlab(std::ostream & os) const
3293 {
3294  os << "% size " << height << " " << width << "\n";
3295  os << "% Non Zeros " << NumNonZeroElems() << "\n";
3296  int i, j;
3297  ios::fmtflags old_fmt = os.flags();
3298  os.setf(ios::scientific);
3299  std::streamsize old_prec = os.precision(14);
3300 
3301  for (i = 0; i < height; i++)
3302  {
3303  for (j = I[i]; j < I[i+1]; j++)
3304  {
3305  os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
3306  }
3307  }
3308  // Write a zero entry at (m,n) to make sure MATLAB doesn't shrink the matrix
3309  os << height << " " << width << " 0.0\n";
3310  os.precision(old_prec);
3311  os.flags(old_fmt);
3312 }
3313 
3314 void SparseMatrix::PrintMM(std::ostream & os) const
3315 {
3316  int i, j;
3317  ios::fmtflags old_fmt = os.flags();
3318  os.setf(ios::scientific);
3319  std::streamsize old_prec = os.precision(14);
3320 
3321  os << "%%MatrixMarket matrix coordinate real general" << '\n'
3322  << "% Generated by MFEM" << '\n';
3323 
3324  os << height << " " << width << " " << NumNonZeroElems() << '\n';
3325  for (i = 0; i < height; i++)
3326  {
3327  for (j = I[i]; j < I[i+1]; j++)
3328  {
3329  os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
3330  }
3331  }
3332  os.precision(old_prec);
3333  os.flags(old_fmt);
3334 }
3335 
3336 void SparseMatrix::PrintCSR(std::ostream & os) const
3337 {
3338  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
3339 
3340  int i;
3341 
3342  os << height << '\n'; // number of rows
3343 
3344  for (i = 0; i <= height; i++)
3345  {
3346  os << I[i]+1 << '\n';
3347  }
3348 
3349  for (i = 0; i < I[height]; i++)
3350  {
3351  os << J[i]+1 << '\n';
3352  }
3353 
3354  for (i = 0; i < I[height]; i++)
3355  {
3356  os << A[i] << '\n';
3357  }
3358 }
3359 
3360 void SparseMatrix::PrintCSR2(std::ostream & os) const
3361 {
3362  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
3363 
3364  int i;
3365 
3366  os << height << '\n'; // number of rows
3367  os << width << '\n'; // number of columns
3368 
3369  for (i = 0; i <= height; i++)
3370  {
3371  os << I[i] << '\n';
3372  }
3373 
3374  for (i = 0; i < I[height]; i++)
3375  {
3376  os << J[i] << '\n';
3377  }
3378 
3379  for (i = 0; i < I[height]; i++)
3380  {
3381  os << A[i] << '\n';
3382  }
3383 }
3384 
3385 void SparseMatrix::PrintInfo(std::ostream &os) const
3386 {
3387  const double MiB = 1024.*1024;
3388  int nnz = NumNonZeroElems();
3389  double pz = 100./nnz;
3390  int nz = CountSmallElems(0.0);
3391  double max_norm = MaxNorm();
3392  double symm = IsSymmetric();
3393  int nnf = CheckFinite();
3394  int ns12 = CountSmallElems(1e-12*max_norm);
3395  int ns15 = CountSmallElems(1e-15*max_norm);
3396  int ns18 = CountSmallElems(1e-18*max_norm);
3397 
3398  os <<
3399  "SparseMatrix statistics:\n"
3400  " Format : " <<
3401  (Empty() ? "(empty)" : (Finalized() ? "CSR" : "LIL")) << "\n"
3402  " Dimensions : " << height << " x " << width << "\n"
3403  " Number of entries (total) : " << nnz << "\n"
3404  " Number of entries (per row) : " << 1.*nnz/Height() << "\n"
3405  " Number of stored zeros : " << nz*pz << "% (" << nz << ")\n"
3406  " Number of Inf/Nan entries : " << nnf*pz << "% ("<< nnf << ")\n"
3407  " Norm, max |a_ij| : " << max_norm << "\n"
3408  " Symmetry, max |a_ij-a_ji| : " << symm << "\n"
3409  " Number of small entries:\n"
3410  " |a_ij| <= 1e-12*Norm : " << ns12*pz << "% (" << ns12 << ")\n"
3411  " |a_ij| <= 1e-15*Norm : " << ns15*pz << "% (" << ns15 << ")\n"
3412  " |a_ij| <= 1e-18*Norm : " << ns18*pz << "% (" << ns18 << ")\n";
3413  if (Finalized())
3414  {
3415  os << " Memory used by CSR : " <<
3416  (sizeof(int)*(height+1+nnz)+sizeof(double)*nnz)/MiB << " MiB\n";
3417  }
3418  if (Rows != NULL)
3419  {
3420  size_t used_mem = sizeof(RowNode*)*height;
3421 #ifdef MFEM_USE_MEMALLOC
3422  used_mem += NodesMem->MemoryUsage();
3423 #else
3424  for (int i = 0; i < height; i++)
3425  {
3426  for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
3427  {
3428  used_mem += sizeof(RowNode);
3429  }
3430  }
3431 #endif
3432  os << " Memory used by LIL : " << used_mem/MiB << " MiB\n";
3433  }
3434 }
3435 
3437 {
3438  I.Delete();
3439  J.Delete();
3440  A.Delete();
3441 
3442  if (Rows != NULL)
3443  {
3444 #if !defined(MFEM_USE_MEMALLOC)
3445  for (int i = 0; i < height; i++)
3446  {
3447  RowNode *aux, *node_p = Rows[i];
3448  while (node_p != NULL)
3449  {
3450  aux = node_p;
3451  node_p = node_p->Prev;
3452  delete aux;
3453  }
3454  }
3455 #endif
3456  delete [] Rows;
3457  }
3458 
3459  delete [] ColPtrJ;
3460  delete [] ColPtrNode;
3461 #ifdef MFEM_USE_MEMALLOC
3462  delete NodesMem;
3463 #endif
3464  delete At;
3465 
3466  ClearGPUSparse();
3467 }
3468 
3470 {
3471  int awidth = 0;
3472  if (A)
3473  {
3474  const int *start_j = J;
3475  const int *end_j = J + I[height];
3476  for (const int *jptr = start_j; jptr != end_j; ++jptr)
3477  {
3478  awidth = std::max(awidth, *jptr + 1);
3479  }
3480  }
3481  else
3482  {
3483  RowNode *aux;
3484  for (int i = 0; i < height; i++)
3485  {
3486  for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
3487  {
3488  awidth = std::max(awidth, aux->Column + 1);
3489  }
3490  }
3491  }
3492  return awidth;
3493 }
3494 
3495 void SparseMatrixFunction (SparseMatrix & S, double (*f)(double))
3496 {
3497  int n = S.NumNonZeroElems();
3498  double * s = S.GetData();
3499 
3500  for (int i = 0; i < n; i++)
3501  {
3502  s[i] = f(s[i]);
3503  }
3504 }
3505 
3507 {
3508  MFEM_VERIFY(
3509  A.Finalized(),
3510  "Finalize must be called before Transpose. Use TransposeRowMatrix instead");
3511 
3512  int i, j, end;
3513  const int *A_i, *A_j;
3514  int m, n, nnz, *At_i, *At_j;
3515  const double *A_data;
3516  double *At_data;
3517 
3518  m = A.Height(); // number of rows of A
3519  n = A.Width(); // number of columns of A
3520  nnz = A.NumNonZeroElems();
3521  A_i = A.HostReadI();
3522  A_j = A.HostReadJ();
3523  A_data = A.HostReadData();
3524 
3525  At_i = Memory<int>(n+1);
3526  At_j = Memory<int>(nnz);
3527  At_data = Memory<double>(nnz);
3528 
3529  for (i = 0; i <= n; i++)
3530  {
3531  At_i[i] = 0;
3532  }
3533  for (i = 0; i < nnz; i++)
3534  {
3535  At_i[A_j[i]+1]++;
3536  }
3537  for (i = 1; i < n; i++)
3538  {
3539  At_i[i+1] += At_i[i];
3540  }
3541 
3542  for (i = j = 0; i < m; i++)
3543  {
3544  end = A_i[i+1];
3545  for ( ; j < end; j++)
3546  {
3547  At_j[At_i[A_j[j]]] = i;
3548  At_data[At_i[A_j[j]]] = A_data[j];
3549  At_i[A_j[j]]++;
3550  }
3551  }
3552 
3553  for (i = n; i > 0; i--)
3554  {
3555  At_i[i] = At_i[i-1];
3556  }
3557  At_i[0] = 0;
3558 
3559  return new SparseMatrix(At_i, At_j, At_data, n, m);
3560 }
3561 
3563  int useActualWidth)
3564 {
3565  int i, j;
3566  int m, n, nnz, *At_i, *At_j;
3567  double *At_data;
3568  Array<int> Acols;
3569  Vector Avals;
3570 
3571  m = A.Height(); // number of rows of A
3572  if (useActualWidth)
3573  {
3574  n = 0;
3575  int tmp;
3576  for (i = 0; i < m; i++)
3577  {
3578  A.GetRow(i, Acols, Avals);
3579  if (Acols.Size())
3580  {
3581  tmp = Acols.Max();
3582  if (tmp > n)
3583  {
3584  n = tmp;
3585  }
3586  }
3587  }
3588  ++n;
3589  }
3590  else
3591  {
3592  n = A.Width(); // number of columns of A
3593  }
3594  nnz = A.NumNonZeroElems();
3595 
3596  At_i = Memory<int>(n+1);
3597  At_j = Memory<int>(nnz);
3598  At_data = Memory<double>(nnz);
3599 
3600  for (i = 0; i <= n; i++)
3601  {
3602  At_i[i] = 0;
3603  }
3604 
3605  for (i = 0; i < m; i++)
3606  {
3607  A.GetRow(i, Acols, Avals);
3608  for (j = 0; j<Acols.Size(); ++j)
3609  {
3610  At_i[Acols[j]+1]++;
3611  }
3612  }
3613  for (i = 1; i < n; i++)
3614  {
3615  At_i[i+1] += At_i[i];
3616  }
3617 
3618  for (i = 0; i < m; i++)
3619  {
3620  A.GetRow(i, Acols, Avals);
3621  for (j = 0; j<Acols.Size(); ++j)
3622  {
3623  At_j[At_i[Acols[j]]] = i;
3624  At_data[At_i[Acols[j]]] = Avals[j];
3625  At_i[Acols[j]]++;
3626  }
3627  }
3628 
3629  for (i = n; i > 0; i--)
3630  {
3631  At_i[i] = At_i[i-1];
3632  }
3633  At_i[0] = 0;
3634 
3635  return new SparseMatrix(At_i, At_j, At_data, n, m);
3636 }
3637 
3638 
3640  SparseMatrix *OAB)
3641 {
3642  int nrowsA, ncolsA, nrowsB, ncolsB;
3643  const int *A_i, *A_j, *B_i, *B_j;
3644  int *C_i, *C_j, *B_marker;
3645  const double *A_data, *B_data;
3646  double *C_data;
3647  int ia, ib, ic, ja, jb, num_nonzeros;
3648  int row_start, counter;
3649  double a_entry, b_entry;
3650  SparseMatrix *C;
3651 
3652  nrowsA = A.Height();
3653  ncolsA = A.Width();
3654  nrowsB = B.Height();
3655  ncolsB = B.Width();
3656 
3657  MFEM_VERIFY(ncolsA == nrowsB,
3658  "number of columns of A (" << ncolsA
3659  << ") must equal number of rows of B (" << nrowsB << ")");
3660 
3661  A_i = A.HostReadI();
3662  A_j = A.HostReadJ();
3663  A_data = A.HostReadData();
3664  B_i = B.HostReadI();
3665  B_j = B.HostReadJ();
3666  B_data = B.HostReadData();
3667 
3668  B_marker = new int[ncolsB];
3669 
3670  for (ib = 0; ib < ncolsB; ib++)
3671  {
3672  B_marker[ib] = -1;
3673  }
3674 
3675  if (OAB == NULL)
3676  {
3677  C_i = Memory<int>(nrowsA+1);
3678 
3679  C_i[0] = num_nonzeros = 0;
3680  for (ic = 0; ic < nrowsA; ic++)
3681  {
3682  for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3683  {
3684  ja = A_j[ia];
3685  for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
3686  {
3687  jb = B_j[ib];
3688  if (B_marker[jb] != ic)
3689  {
3690  B_marker[jb] = ic;
3691  num_nonzeros++;
3692  }
3693  }
3694  }
3695  C_i[ic+1] = num_nonzeros;
3696  }
3697 
3698  C_j = Memory<int>(num_nonzeros);
3699  C_data = Memory<double>(num_nonzeros);
3700 
3701  C = new SparseMatrix(C_i, C_j, C_data, nrowsA, ncolsB);
3702 
3703  for (ib = 0; ib < ncolsB; ib++)
3704  {
3705  B_marker[ib] = -1;
3706  }
3707  }
3708  else
3709  {
3710  C = OAB;
3711 
3712  MFEM_VERIFY(nrowsA == C->Height() && ncolsB == C->Width(),
3713  "Input matrix sizes do not match output sizes"
3714  << " nrowsA = " << nrowsA
3715  << ", C->Height() = " << C->Height()
3716  << " ncolsB = " << ncolsB
3717  << ", C->Width() = " << C->Width());
3718 
3719  // C_i = C->HostReadI(); // not used
3720  C_j = C->HostWriteJ();
3721  C_data = C->HostWriteData();
3722  }
3723 
3724  counter = 0;
3725  for (ic = 0; ic < nrowsA; ic++)
3726  {
3727  // row_start = C_i[ic];
3728  row_start = counter;
3729  for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3730  {
3731  ja = A_j[ia];
3732  a_entry = A_data[ia];
3733  for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
3734  {
3735  jb = B_j[ib];
3736  b_entry = B_data[ib];
3737  if (B_marker[jb] < row_start)
3738  {
3739  B_marker[jb] = counter;
3740  if (OAB == NULL)
3741  {
3742  C_j[counter] = jb;
3743  }
3744  C_data[counter] = a_entry*b_entry;
3745  counter++;
3746  }
3747  else
3748  {
3749  C_data[B_marker[jb]] += a_entry*b_entry;
3750  }
3751  }
3752  }
3753  }
3754 
3755  MFEM_VERIFY(
3756  OAB == NULL || counter == OAB->NumNonZeroElems(),
3757  "With pre-allocated output matrix, number of non-zeros ("
3758  << OAB->NumNonZeroElems()
3759  << ") did not match number of entries changed from matrix-matrix multiply, "
3760  << counter);
3761 
3762  delete [] B_marker;
3763 
3764  return C;
3765 }
3766 
3768 {
3769  SparseMatrix *At = Transpose(A);
3770  SparseMatrix *AtB = Mult(*At, B);
3771  delete At;
3772  return AtB;
3773 }
3774 
3776  const AbstractSparseMatrix &B)
3777 {
3778  int nrowsA, ncolsA, nrowsB, ncolsB;
3779  int *C_i, *C_j, *B_marker;
3780  double *C_data;
3781  int ia, ib, ic, ja, jb, num_nonzeros;
3782  int row_start, counter;
3783  double a_entry, b_entry;
3784  SparseMatrix *C;
3785 
3786  nrowsA = A.Height();
3787  ncolsA = A.Width();
3788  nrowsB = B.Height();
3789  ncolsB = B.Width();
3790 
3791  MFEM_VERIFY(ncolsA == nrowsB,
3792  "number of columns of A (" << ncolsA
3793  << ") must equal number of rows of B (" << nrowsB << ")");
3794 
3795  B_marker = new int[ncolsB];
3796 
3797  for (ib = 0; ib < ncolsB; ib++)
3798  {
3799  B_marker[ib] = -1;
3800  }
3801 
3802  C_i = Memory<int>(nrowsA+1);
3803 
3804  C_i[0] = num_nonzeros = 0;
3805 
3806  Array<int> colsA, colsB;
3807  Vector dataA, dataB;
3808  for (ic = 0; ic < nrowsA; ic++)
3809  {
3810  A.GetRow(ic, colsA, dataA);
3811  for (ia = 0; ia < colsA.Size(); ia++)
3812  {
3813  ja = colsA[ia];
3814  B.GetRow(ja, colsB, dataB);
3815  for (ib = 0; ib < colsB.Size(); ib++)
3816  {
3817  jb = colsB[ib];
3818  if (B_marker[jb] != ic)
3819  {
3820  B_marker[jb] = ic;
3821  num_nonzeros++;
3822  }
3823  }
3824  }
3825  C_i[ic+1] = num_nonzeros;
3826  }
3827 
3828  C_j = Memory<int>(num_nonzeros);
3829  C_data = Memory<double>(num_nonzeros);
3830 
3831  C = new SparseMatrix(C_i, C_j, C_data, nrowsA, ncolsB);
3832 
3833  for (ib = 0; ib < ncolsB; ib++)
3834  {
3835  B_marker[ib] = -1;
3836  }
3837 
3838  counter = 0;
3839  for (ic = 0; ic < nrowsA; ic++)
3840  {
3841  row_start = counter;
3842  A.GetRow(ic, colsA, dataA);
3843  for (ia = 0; ia < colsA.Size(); ia++)
3844  {
3845  ja = colsA[ia];
3846  a_entry = dataA[ia];
3847  B.GetRow(ja, colsB, dataB);
3848  for (ib = 0; ib < colsB.Size(); ib++)
3849  {
3850  jb = colsB[ib];
3851  b_entry = dataB[ib];
3852  if (B_marker[jb] < row_start)
3853  {
3854  B_marker[jb] = counter;
3855  C_j[counter] = jb;
3856  C_data[counter] = a_entry*b_entry;
3857  counter++;
3858  }
3859  else
3860  {
3861  C_data[B_marker[jb]] += a_entry*b_entry;
3862  }
3863  }
3864  }
3865  }
3866 
3867  delete [] B_marker;
3868 
3869  return C;
3870 }
3871 
3873 {
3874  DenseMatrix *C = new DenseMatrix(A.Height(), B.Width());
3875  Vector columnB, columnC;
3876  for (int j = 0; j < B.Width(); ++j)
3877  {
3878  B.GetColumnReference(j, columnB);
3879  C->GetColumnReference(j, columnC);
3880  A.Mult(columnB, columnC);
3881  }
3882  return C;
3883 }
3884 
3886 {
3887  DenseMatrix R (P, 't'); // R = P^T
3888  DenseMatrix *AP = Mult (A, P);
3889  DenseMatrix *RAP_ = new DenseMatrix(R.Height(), AP->Width());
3890  Mult (R, *AP, *RAP_);
3891  delete AP;
3892  return RAP_;
3893 }
3894 
3896 {
3897  SparseMatrix *R = Transpose(P);
3898  DenseMatrix *RA = Mult(*R, A);
3899  DenseMatrix AtP(*RA, 't');
3900  delete RA;
3901  DenseMatrix *RAtP = Mult(*R, AtP);
3902  delete R;
3903  DenseMatrix * RAP_ = new DenseMatrix(*RAtP, 't');
3904  delete RAtP;
3905  return RAP_;
3906 }
3907 
3909  SparseMatrix *ORAP)
3910 {
3911  SparseMatrix *P = Transpose (R);
3912  SparseMatrix *AP = Mult (A, *P);
3913  delete P;
3914  SparseMatrix *RAP_ = Mult (R, *AP, ORAP);
3915  delete AP;
3916  return RAP_;
3917 }
3918 
3920  const SparseMatrix &P)
3921 {
3922  SparseMatrix * R = Transpose(Rt);
3923  SparseMatrix * RA = Mult(*R,A);
3924  delete R;
3925  SparseMatrix * RAP_ = Mult(*RA, P);
3926  delete RA;
3927  return RAP_;
3928 }
3929 
3931  SparseMatrix *OAtDA)
3932 {
3933  int i, At_nnz, *At_j;
3934  double *At_data;
3935 
3936  SparseMatrix *At = Transpose (A);
3937  At_nnz = At -> NumNonZeroElems();
3938  At_j = At -> GetJ();
3939  At_data = At -> GetData();
3940  for (i = 0; i < At_nnz; i++)
3941  {
3942  At_data[i] *= D(At_j[i]);
3943  }
3944  SparseMatrix *AtDA = Mult (*At, A, OAtDA);
3945  delete At;
3946  return AtDA;
3947 }
3948 
3949 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
3950  const SparseMatrix & B)
3951 {
3952  int nrows = A.Height();
3953  int ncols = A.Width();
3954 
3955  int * C_i = Memory<int>(nrows+1);
3956  int * C_j;
3957  double * C_data;
3958 
3959  const int *A_i = A.GetI();
3960  const int *A_j = A.GetJ();
3961  const double *A_data = A.GetData();
3962 
3963  const int *B_i = B.GetI();
3964  const int *B_j = B.GetJ();
3965  const double *B_data = B.GetData();
3966 
3967  int * marker = new int[ncols];
3968  std::fill(marker, marker+ncols, -1);
3969 
3970  int num_nonzeros = 0, jcol;
3971  C_i[0] = 0;
3972  for (int ic = 0; ic < nrows; ic++)
3973  {
3974  for (int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3975  {
3976  jcol = A_j[ia];
3977  marker[jcol] = ic;
3978  num_nonzeros++;
3979  }
3980  for (int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
3981  {
3982  jcol = B_j[ib];
3983  if (marker[jcol] != ic)
3984  {
3985  marker[jcol] = ic;
3986  num_nonzeros++;
3987  }
3988  }
3989  C_i[ic+1] = num_nonzeros;
3990  }
3991 
3992  C_j = Memory<int>(num_nonzeros);
3993  C_data = Memory<double>(num_nonzeros);
3994 
3995  for (int ia = 0; ia < ncols; ia++)
3996  {
3997  marker[ia] = -1;
3998  }
3999 
4000  int pos = 0;
4001  for (int ic = 0; ic < nrows; ic++)
4002  {
4003  for (int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
4004  {
4005  jcol = A_j[ia];
4006  C_j[pos] = jcol;
4007  C_data[pos] = a*A_data[ia];
4008  marker[jcol] = pos;
4009  pos++;
4010  }
4011  for (int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
4012  {
4013  jcol = B_j[ib];
4014  if (marker[jcol] < C_i[ic])
4015  {
4016  C_j[pos] = jcol;
4017  C_data[pos] = b*B_data[ib];
4018  marker[jcol] = pos;
4019  pos++;
4020  }
4021  else
4022  {
4023  C_data[marker[jcol]] += b*B_data[ib];
4024  }
4025  }
4026  }
4027 
4028  delete[] marker;
4029  return new SparseMatrix(C_i, C_j, C_data, nrows, ncols);
4030 }
4031 
4033 {
4034  return Add(1.,A,1.,B);
4035 }
4036 
4038 {
4039  MFEM_ASSERT(Ai.Size() > 0, "invalid size Ai.Size() = " << Ai.Size());
4040 
4041  SparseMatrix * accumulate = Ai[0];
4042  SparseMatrix * result = accumulate;
4043 
4044  for (int i=1; i < Ai.Size(); ++i)
4045  {
4046  result = Add(*accumulate, *Ai[i]);
4047  if (i != 1)
4048  {
4049  delete accumulate;
4050  }
4051 
4052  accumulate = result;
4053  }
4054 
4055  return result;
4056 }
4057 
4058 /// B += alpha * A
4059 void Add(const SparseMatrix &A,
4060  double alpha, DenseMatrix &B)
4061 {
4062  for (int r = 0; r < B.Height(); r++)
4063  {
4064  const int * colA = A.GetRowColumns(r);
4065  const double * valA = A.GetRowEntries(r);
4066  for (int i=0; i<A.RowSize(r); i++)
4067  {
4068  B(r, colA[i]) += alpha * valA[i];
4069  }
4070  }
4071 }
4072 
4073 /// Produces a block matrix with blocks A_{ij}*B
4075 {
4076  int mA = A.Height(), nA = A.Width();
4077  int mB = B.Height(), nB = B.Width();
4078 
4079  DenseMatrix *C = new DenseMatrix(mA * mB, nA * nB);
4080  *C = 0.0;
4081  for (int i=0; i<mA; i++)
4082  {
4083  for (int j=0; j<nA; j++)
4084  {
4085  C->AddMatrix(A(i,j), B, i * mB, j * nB);
4086  }
4087  }
4088  return C;
4089 }
4090 
4091 /// Produces a block matrix with blocks A_{ij}*B
4093 {
4094  int mA = A.Height(), nA = A.Width();
4095  int mB = B.Height(), nB = B.Width();
4096 
4097  SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
4098 
4099  for (int i=0; i<mA; i++)
4100  {
4101  for (int j=0; j<nA; j++)
4102  {
4103  for (int r=0; r<mB; r++)
4104  {
4105  const int * colB = B.GetRowColumns(r);
4106  const double * valB = B.GetRowEntries(r);
4107 
4108  for (int cj=0; cj<B.RowSize(r); cj++)
4109  {
4110  C->Set(i * mB + r, j * nB + colB[cj], A(i,j) * valB[cj]);
4111  }
4112  }
4113  }
4114  }
4115  C->Finalize();
4116 
4117  return C;
4118 }
4119 
4120 /// Produces a block matrix with blocks A_{ij}*B
4122 {
4123  int mA = A.Height(), nA = A.Width();
4124  int mB = B.Height(), nB = B.Width();
4125 
4126  SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
4127 
4128  for (int r=0; r<mA; r++)
4129  {
4130  const int * colA = A.GetRowColumns(r);
4131  const double * valA = A.GetRowEntries(r);
4132 
4133  for (int aj=0; aj<A.RowSize(r); aj++)
4134  {
4135  for (int i=0; i<mB; i++)
4136  {
4137  for (int j=0; j<nB; j++)
4138  {
4139  C->Set(r * mB + i, colA[aj] * nB + j, valA[aj] * B(i, j));
4140  }
4141  }
4142  }
4143  }
4144  C->Finalize();
4145 
4146  return C;
4147 }
4148 
4149 /// Produces a block matrix with blocks A_{ij}*B
4151 {
4152  int mA = A.Height(), nA = A.Width();
4153  int mB = B.Height(), nB = B.Width();
4154 
4155  SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
4156 
4157  for (int ar=0; ar<mA; ar++)
4158  {
4159  const int * colA = A.GetRowColumns(ar);
4160  const double * valA = A.GetRowEntries(ar);
4161 
4162  for (int aj=0; aj<A.RowSize(ar); aj++)
4163  {
4164  for (int br=0; br<mB; br++)
4165  {
4166  const int * colB = B.GetRowColumns(br);
4167  const double * valB = B.GetRowEntries(br);
4168 
4169  for (int bj=0; bj<B.RowSize(br); bj++)
4170  {
4171  C->Set(ar * mB + br, colA[aj] * nB + colB[bj],
4172  valA[aj] * valB[bj]);
4173  }
4174  }
4175  }
4176  }
4177  C->Finalize();
4178 
4179  return C;
4180 }
4181 
4183 {
4184  mfem::Swap(width, other.width);
4185  mfem::Swap(height, other.height);
4186  mfem::Swap(I, other.I);
4187  mfem::Swap(J, other.J);
4188  mfem::Swap(A, other.A);
4189  mfem::Swap(Rows, other.Rows);
4191  mfem::Swap(ColPtrJ, other.ColPtrJ);
4193  mfem::Swap(At, other.At);
4194 
4195 #ifdef MFEM_USE_MEMALLOC
4196  mfem::Swap(NodesMem, other.NodesMem);
4197 #endif
4198 
4199  mfem::Swap(isSorted, other.isSorted);
4200 }
4201 
4203 {
4204  Destroy();
4205 #ifdef MFEM_USE_CUDA_OR_HIP
4207  {
4208  if (SparseMatrixCount==1)
4209  {
4210  if (handle)
4211  {
4212  MFEM_cu_or_hip(sparseDestroy)(handle);
4213  handle = nullptr;
4214  }
4215  if (dBuffer)
4216  {
4217  MFEM_Cu_or_Hip(MemFree)(dBuffer);
4218  dBuffer = nullptr;
4219  bufferSize = 0;
4220  }
4221  }
4223  }
4224 #endif // MFEM_USE_CUDA_OR_HIP
4225 }
4226 
4227 }
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
RowNode ** ColPtrNode
Definition: sparsemat.hpp:77
Memory< int > I
Array with size (height+1) containing the row offsets.
Definition: sparsemat.hpp:62
Elem * Alloc()
Definition: mem_alloc.hpp:166
const double * HostReadData() const
Definition: sparsemat.hpp:277
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
Definition: sparsemat.cpp:3562
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:492
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1588
void * CuMemFree(void *dptr)
Frees device memory and returns destination ptr.
Definition: cuda.cpp:79
void _Add_(const int col, const double a)
Add a value to an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:585
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:162
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:209
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
Definition: sparsemat.cpp:1741
Memory< T > & GetMemory()
Return a reference to the Memory object used by the Array.
Definition: array.hpp:120
SparseMatrix * At
Transpose of A. Owned. Used to perform MultTranspose() on devices.
Definition: sparsemat.hpp:80
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:669
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
Definition: sparsemat.cpp:3314
static size_t bufferSize
Definition: sparsemat.hpp:101
double GetRowNorml1(int irow) const
For i = irow compute .
Definition: sparsemat.cpp:1237
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:546
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:313
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
double & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:846
void Delete()
Delete the owned pointers and reset the Memory object.
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:452
size_t MemoryUsage() const
Definition: mem_alloc.hpp:212
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
Definition: sparsemat.cpp:3336
double GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
Definition: sparsemat.cpp:2555
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void SetColPtr(const int row) const
Initialize the SparseMatrix for fast access to the entries of the given row which becomes the "curren...
Definition: sparsemat.hpp:792
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Biwise-OR of all HIP backends.
Definition: device.hpp:90
Abstract data type for sparse matrices.
Definition: matrix.hpp:73
double IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
Definition: sparsemat.cpp:1527
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
SparseMatrix & operator+=(const SparseMatrix &B)
Add the sparse matrix &#39;B&#39; to &#39;*this&#39;. This operation will cause an error if &#39;*this&#39; is finalized and ...
Definition: sparsemat.cpp:3141
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:336
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
static void * dBuffer
Definition: sparsemat.hpp:102
void EliminateRowColMultipleRHS(int rc, const Vector &sol, DenseMatrix &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but multiple values for elim...
Definition: sparsemat.cpp:1953
virtual int NumNonZeroElems() const =0
Returns the number of non-zeros in a matrix.
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:740
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Extract all column indices and values from a given row.
Definition: sparsemat.cpp:2935
void ClearGPUSparse()
Clear the cuSPARSE/hipSPARSE descriptors. This must be called after releasing the device memory of A...
Definition: sparsemat.cpp:81
STL namespace.
double MaxNorm() const
Definition: sparsemat.cpp:1611
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: sparsemat.cpp:601
void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: sparsemat.cpp:3244
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:3022
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:405
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:457
cusparseDnVecDescr_t vecY_descr
Definition: sparsemat.hpp:113
double * HostReadWriteData()
Definition: sparsemat.hpp:281
void _Set_(const int col, const double a)
Set an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:588
SparseMatrix * Mult_AtDA(const SparseMatrix &A, const Vector &D, SparseMatrix *OAtDA)
Matrix multiplication A^t D A. All matrices must be finalized.
Definition: sparsemat.cpp:3930
unsigned flags
Bit flags defined from the FlagMask enum.
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:2472
void ScaleRow(const int row, const double scale)
Definition: sparsemat.cpp:3054
RowNodeAlloc * NodesMem
Definition: sparsemat.hpp:84
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
virtual MatrixInverse * Inverse() const
This virtual method is not supported: it always returns NULL.
Definition: sparsemat.cpp:1688
DenseMatrix * ToDenseMatrix() const
Produces a DenseMatrix from a SparseMatrix.
Definition: sparsemat.cpp:703
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: sparsemat.cpp:1569
bool IsFinite(const double &val)
Definition: vector.hpp:481
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2715
void MoveDiagonalFirst()
Move the diagonal entry to the first position in each row, preserving the order of the rest of the co...
Definition: sparsemat.cpp:577
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2321
void Set(const int i, const int j, const double val)
Definition: sparsemat.cpp:2768
void Wrap(T *ptr, int size, bool own)
Wrap an externally allocated host pointer, ptr with the current host memory type returned by MemoryMa...
int NumRows() const
Definition: array.hpp:377
void GetBlocks(Array2D< SparseMatrix *> &blocks) const
Definition: sparsemat.cpp:1437
void ScaleColumns(const Vector &sr)
this = this * diag(sr);
Definition: sparsemat.cpp:3113
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:1214
virtual ~SparseMatrix()
Destroys sparse matrix.
Definition: sparsemat.cpp:4202
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
bool Empty() const
Check if the SparseMatrix is empty.
Definition: sparsemat.hpp:219
bool isSorted
Are the columns sorted already.
Definition: sparsemat.hpp:88
void ClearColPtr() const
Reset the "current row" set by calling SetColPtr(). This method must be called between any two calls ...
Definition: sparsemat.hpp:827
Memory< double > A
Array with size I[height], containing the actual entries of the sparse matrix, as indexed by the I ar...
Definition: sparsemat.hpp:68
SparseMatrix()
Create an empty SparseMatrix.
Definition: sparsemat.hpp:131
Data type sparse matrix.
Definition: sparsemat.hpp:50
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:456
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
Definition: sparsemat.cpp:419
double b
Definition: lissajous.cpp:42
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(), and AbsMultTranspose().
Definition: sparsemat.cpp:974
const int * HostReadJ() const
Definition: sparsemat.hpp:261
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:3292
int CountSmallElems(double tol) const
Count the number of entries with |a_ij| <= tol.
Definition: sparsemat.cpp:1636
void * HipMemFree(void *dptr)
Frees device memory.
Definition: hip.cpp:79
static int SparseMatrixCount
Definition: sparsemat.hpp:100
const double * ReadData(bool on_dev=true) const
Definition: sparsemat.hpp:271
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:1173
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
Definition: sparsemat.cpp:2358
void ScaleRows(const Vector &sl)
this = diag(sl) * this;
Definition: sparsemat.cpp:3082
static cusparseHandle_t handle
Definition: sparsemat.hpp:107
DenseMatrix * OuterProduct(const DenseMatrix &A, const DenseMatrix &B)
Produces a block matrix with blocks A_{ij}*B.
Definition: sparsemat.cpp:4074
void DiagScale(const Vector &b, Vector &x, double sc=1.0, bool use_abs_diag=false) const
x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
Definition: sparsemat.cpp:2620
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition: array.hpp:251
int Capacity() const
Return the size of the allocated memory.
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const =0
Gets the columns indexes and values for row row.
Set the diagonal value to one.
Definition: operator.hpp:50
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
bool RowIsEmpty(const int row) const
Definition: sparsemat.cpp:2914
void Reset()
Reset the memory to be empty, ensuring that Delete() will be a no-op.
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2708
Biwise-OR of all CUDA backends.
Definition: device.hpp:88
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2722
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2787
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2806
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
Biwise-OR of all CPU backends.
Definition: device.hpp:86
void * HipMemAlloc(void **dptr, size_t bytes)
Allocates device memory.
Definition: hip.cpp:34
T * HostWrite(Memory< T > &mem, int size)
Shortcut to Write(const Memory<T> &mem, int size, false)
Definition: device.hpp:343
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
void SparseMatrixFunction(SparseMatrix &S, double(*f)(double))
Applies f() to each element of the matrix (after it is finalized).
Definition: sparsemat.cpp:3495
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:635
const int * ReadI(bool on_dev=true) const
Definition: sparsemat.hpp:239
double _Get_(const int col) const
Read the value of an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:873
void EliminateRowColDiag(int rc, double value)
Perform elimination and set the diagonal entry to the given value.
Definition: sparsemat.cpp:2158
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2845
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:2388
void SetHostPtrOwner(bool own) const
Set/clear the ownership flag for the host pointer. Ownership indicates whether the pointer will be de...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
int ActualWidth() const
Returns the actual Width of the matrix.
Definition: sparsemat.cpp:3469
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1854
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
void Threshold(double tol, bool fix_empty_rows=false)
Remove entries smaller in absolute value than a given tolerance tol. If fix_empty_rows is true...
Definition: sparsemat.cpp:1261
double * HostWriteData()
Definition: sparsemat.hpp:279
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width.
Definition: densemat.cpp:1723
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Definition: device.hpp:258
int * HostReadWriteJ()
Definition: sparsemat.hpp:265
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format...
Definition: sparsemat.hpp:73
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc, bool use_abs_diag=false) const
Definition: sparsemat.cpp:2588
const int * HostReadI() const
Definition: sparsemat.hpp:245
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1693
double a
Definition: lissajous.cpp:41
void BuildTranspose() const
Build and store internally the transpose of this matrix which will be used in the methods AddMultTran...
Definition: sparsemat.cpp:960
const T * HostRead(const Memory< T > &mem, int size)
Shortcut to Read(const Memory<T> &mem, int size, false)
Definition: device.hpp:326
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:464
MemAlloc< RowNode, 1024 > RowNodeAlloc
Definition: sparsemat.hpp:83
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
Definition: sparsemat.cpp:3767
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Keep the diagonal value.
Definition: operator.hpp:51
virtual void EliminateZeroRows(const double threshold=1e-12)
If a row contains only zeros, set its diagonal to 1.
Definition: sparsemat.cpp:2369
int NumCols() const
Definition: array.hpp:378
void ClearOwnerFlags() const
Clear the ownership flags for the host and device pointers, as well as any internal data allocated by...
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
double & operator()(int i, int j)
Returns reference to A[i][j].
Definition: sparsemat.cpp:611
void AbsMult(const Vector &x, Vector &y) const
y = |A| * x, using entry-wise absolute values of matrix A
Definition: sparsemat.cpp:1079
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
Definition: sparsemat.cpp:3385
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
Definition: sparsemat.cpp:1781
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Definition: sparsemat.cpp:297
const int * ReadJ(bool on_dev=true) const
Definition: sparsemat.hpp:255
SparseMatrix & operator=(const SparseMatrix &rhs)
Assignment operator: deep copy.
Definition: sparsemat.cpp:303
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
Definition: sparsemat.cpp:1128
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
SparseMatrix & operator*=(double a)
Definition: sparsemat.cpp:3220
void ResetTranspose() const
Definition: sparsemat.cpp:968
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
Definition: sparsemat.cpp:3360
int MaxRowSize() const
Returns the maximum number of elements among all rows.
Definition: sparsemat.cpp:367
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2975
cusparseSpMatDescr_t matA_descr
Definition: sparsemat.hpp:111
T * Write(bool on_dev=true)
Shortcut for mfem::Write(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:315
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:869
const double alpha
Definition: ex15.cpp:369
double * ReadWriteData(bool on_dev=true)
Definition: sparsemat.hpp:275
int * ReadWriteJ(bool on_dev=true)
Definition: sparsemat.hpp:259
Vector data type.
Definition: vector.hpp:60
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Definition: sparsemat.cpp:1008
cusparseDnVecDescr_t vecX_descr
Definition: sparsemat.hpp:112
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
Definition: sparsemat.cpp:1664
int * HostReadWriteI()
Definition: sparsemat.hpp:249
RefCoord s[3]
SparseMatrix * MultAbstractSparseMatrix(const AbstractSparseMatrix &A, const AbstractSparseMatrix &B)
Matrix product of sparse matrices. A and B do not need to be CSR matrices.
Definition: sparsemat.cpp:3775
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Memory< int > J
Array with size I[height], containing the column indices for all matrix entries, as indexed by the I ...
Definition: sparsemat.hpp:65
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:4182
void * CuMemAlloc(void **dptr, size_t bytes)
Allocates device memory and returns destination ptr.
Definition: cuda.cpp:34
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:468
Set the diagonal value to zero.
Definition: operator.hpp:49
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
void EliminateBC(const Array< int > &ess_dofs, DiagonalPolicy diag_policy)
Eliminate essential (Dirichlet) boundary conditions.
Definition: sparsemat.cpp:2315
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += At * x (default) or y += a * At * x
Definition: sparsemat.cpp:915
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:982
void Neg()
(*this) = -(*this)
Definition: vector.cpp:306
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1056