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