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