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