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