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