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