19 #include "../general/table.hpp"
20 #include "../general/sort_pairs.hpp"
32 Rows(new RowNode *[nrows]),
40 for (
int i = 0; i < nrows; i++)
45 #ifdef MFEM_USE_MEMALLOC
62 #ifdef MFEM_USE_MEMALLOC
68 bool ownij,
bool owna,
bool issorted)
80 #ifdef MFEM_USE_MEMALLOC
87 A =
new double[ I[
height] ];
99 RowNode *row = Rows[i];
100 for ( ; row != NULL; row = row->Prev)
101 if (row->Value != 0.0)
114 for (
int i=0; i <
height; ++i)
116 rowSize = I[i+1]-I[i];
117 out = (out > rowSize) ? out : rowSize;
122 for (
int i=0; i <
height; ++i)
125 out = (out > rowSize) ? out : rowSize;
134 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
141 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
148 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
155 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
162 if (newWidth ==
width)
167 else if ( newWidth == -1)
173 else if (newWidth >
width)
178 delete [] ColPtrNode;
179 ColPtrNode =
static_cast<RowNode **
>(NULL);
184 ColPtrJ =
static_cast<int *
>(NULL);
192 "The new width needs to be bigger or equal to the actual width");
200 MFEM_VERIFY(
Finalized(),
"Matrix is not Finalized!");
208 for (
int j = 0, i = 0; i <
height; i++)
212 for (
int k = 0; k < row.
Size(); k++)
218 for (
int k = 0; k < row.Size(); k++, j++)
241 MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
242 "Trying to access element outside of the matrix. "
243 <<
"height = " <<
height <<
", "
244 <<
"width = " <<
width <<
", "
245 <<
"i = " << i <<
", "
248 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
251 for (k = I[i]; k < end; k++)
257 MFEM_ABORT(
"Did not find i = " << i <<
", j = " << j <<
" in matrix.");
264 static const double zero = 0.0;
266 MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
267 "Trying to access element outside of the matrix. "
268 <<
"height = " <<
height <<
", "
269 <<
"width = " <<
width <<
", "
270 <<
"i = " << i <<
", "
273 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
275 for (k = I[i]; k < end; k++)
287 "Matrix must be square, not height = " <<
height <<
", width = " <<
width);
288 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
293 for (
int i = 0; i <
height; i++)
297 for (j = I[i]; j < end; j++)
321 "Input vector size (" << x.
Size() <<
") must match matrix width (" <<
width
324 "Output vector size (" << y.
Size() <<
") must match matrix height (" <<
height
328 double *Ap = A, *yp = y.
GetData();
329 const double *xp = x.
GetData();
334 for (i = 0; i <
height; i++)
336 RowNode *row = Rows[i];
338 for ( ; row != NULL; row = row->Prev)
340 b += row->Value * xp[row->Column];
348 int *Jp = J, *Ip = I;
352 #ifndef MFEM_USE_OPENMP
353 for (i = j = 0; i <
height; i++)
356 for (end = Ip[i+1]; j < end; j++)
358 d += Ap[j] * xp[Jp[j]];
363 #pragma omp parallel for private(j,end)
364 for (i = 0; i <
height; i++)
367 for (j = Ip[i], end = Ip[i+1]; j < end; j++)
369 d += Ap[j] * xp[Jp[j]];
376 for (i = j = 0; i <
height; i++)
379 for (end = Ip[i+1]; j < end; j++)
381 d += Ap[j] * xp[Jp[j]];
394 const double a)
const
397 "Input vector size (" << x.
Size() <<
") must match matrix height (" <<
height
400 "Output vector size (" << y.
Size() <<
") must match matrix width (" <<
width
409 for (i = 0; i <
height; i++)
411 RowNode *row = Rows[i];
413 for ( ; row != NULL; row = row->Prev)
415 yp[row->Column] += row->Value * b;
421 for (i = 0; i <
height; i++)
423 double xi = a * x(i);
425 for (j = I[i]; j < end; j++)
435 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
437 for (
int i = 0; i < rows.
Size(); i++)
442 for (
int j = I[r]; j < end; j++)
453 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
455 for (
int i = 0; i < rows.
Size(); i++)
460 for (
int j = I[r]; j < end; j++)
462 val += A[j] * x(J[j]);
471 for (
int i = 0; i <
height; i++)
475 for (
int j = I[i], end = I[i+1]; j < end; j++)
480 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
482 a += np->Value * x(np->Column);
492 for (
int i = 0; i <
height; i++)
496 for (
int j = I[i], end = I[i+1]; j < end; j++)
501 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
511 MFEM_VERIFY(irow <
height,
512 "row " << irow <<
" not in matrix with height " <<
height);
516 for (
int j = I[irow], end = I[irow+1]; j < end; j++)
521 for (RowNode *np = Rows[irow]; np != NULL; np = np->Prev)
523 a += fabs(np->Value);
539 delete [] ColPtrNode;
544 for (i = 1; i <=
height; i++)
547 for (aux = Rows[i-1]; aux != NULL; aux = aux->Prev)
548 if (!skip_zeros || aux->Value != 0.0)
560 for (j = i = 0; i <
height; i++)
563 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
565 if (!skip_zeros || aux->Value != 0.0)
570 if ( lastCol > J[j] )
581 #ifdef MFEM_USE_MEMALLOC
585 for (i = 0; i <
height; i++)
587 RowNode *node_p = Rows[i];
588 while (node_p != NULL)
591 node_p = node_p->Prev;
603 MFEM_VERIFY(!
Finalized(),
"Matrix must NOT be finalized.");
606 int nr = (
height + br - 1)/br, nc = (
width + bc - 1)/bc;
608 for (
int j = 0; j < bc; j++)
609 for (
int i = 0; i < br; i++)
611 int *bI =
new int[nr + 1];
612 for (
int k = 0; k <= nr; k++)
619 for (
int gr = 0; gr <
height; gr++)
621 int bi = gr/nr, i = gr%nr + 1;
622 for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
623 if (n_p->Value != 0.0)
625 blocks(bi,n_p->Column/nc)->I[i]++;
629 for (
int j = 0; j < bc; j++)
630 for (
int i = 0; i < br; i++)
634 for (
int k = 1; k <= nr; k++)
636 rs = b.I[k], b.I[k] = nnz, nnz += rs;
639 b.A =
new double[nnz];
642 for (
int gr = 0; gr <
height; gr++)
644 int bi = gr/nr, i = gr%nr + 1;
645 for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
646 if (n_p->Value != 0.0)
649 b.J[b.I[i]] = n_p->Column % nc;
650 b.A[b.I[i]] = n_p->Value;
658 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
664 for (i = 1; i <
height; i++)
665 for (j = I[i]; j < I[i+1]; j++)
668 a = fabs ( A[j] - (*
this)(J[j],i) );
680 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
683 for (i = 1; i <
height; i++)
684 for (j = I[i]; j < I[i+1]; j++)
687 A[j] += (*this)(J[j],i);
689 (*this)(J[j],i) = A[j];
703 for (
int i = 0; i <
height; i++)
704 for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
720 for (
int j = 0; j < nnz; j++)
722 m = fmax(m, fabs(A[j]));
727 for (
int i = 0; i <
height; i++)
728 for (RowNode *n_p = Rows[i]; n_p != NULL; n_p = n_p->Prev)
730 m = fmax(m, fabs(n_p->Value));
745 for (i = 0; i < nz; i++)
746 if (fabs(Ap[i]) < tol)
755 for (i = 0; i <
height; i++)
756 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
757 if (fabs(aux -> Value) < tol)
775 MFEM_ASSERT(row < height && row >= 0,
776 "Row " << row <<
" not in matrix of height " <<
height);
778 MFEM_VERIFY(!
Finalized(),
"Matrix must NOT be finalized.");
780 for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
782 rhs(aux->Column) -= sol * aux->Value;
791 MFEM_ASSERT(row < height && row >= 0,
792 "Row " << row <<
" not in matrix of height " <<
height);
794 "if setOneDiagonal, must be rectangular matrix, not height = "
798 for (
int i=I[row]; i < I[row+1]; ++i)
803 for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
810 SearchRow(row, row) = 1.;
818 MFEM_VERIFY(!
Finalized(),
"Matrix must NOT be finalized.");
820 for (
int i = 0; i <
height; i++)
821 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
822 if (aux -> Column == col)
832 for (
int i = 0; i <
height; i++)
833 for (
int jpos = I[i]; jpos != I[i+1]; ++jpos)
838 (*b)(i) -= A[jpos] * (*x)( J[jpos] );
846 for (
int i = 0; i <
height; i++)
847 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
848 if (cols[aux -> Column])
852 (*b)(i) -= aux -> Value * (*x)(aux -> Column);
864 MFEM_ASSERT(rc < height && rc >= 0,
865 "Row " << rc <<
" not in matrix of height " <<
height);
868 for (
int j = I[rc]; j < I[rc+1]; j++)
869 if ((col = J[j]) == rc)
872 rhs(rc) = A[j] * sol;
882 for (
int k = I[col]; 1; k++)
885 mfem_error(
"SparseMatrix::EliminateRowCol () #2");
889 rhs(col) -= sol * A[k];
895 for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
896 if ((col = aux->Column) == rc)
899 rhs(rc) = aux->Value * sol;
909 for (RowNode *node = Rows[col]; 1; node = node->Prev)
912 mfem_error(
"SparseMatrix::EliminateRowCol () #3");
914 else if (node->Column == rc)
916 rhs(col) -= sol * node->Value;
927 int num_rhs = rhs.
Width();
929 MFEM_ASSERT(rc < height && rc >= 0,
930 "Row " << rc <<
" not in matrix of height " <<
height);
931 MFEM_ASSERT(sol.
Size() == num_rhs,
"solution size (" << sol.
Size()
932 <<
") must match rhs width (" << num_rhs <<
")");
935 for (
int j = I[rc]; j < I[rc+1]; j++)
936 if ((col = J[j]) == rc)
939 for (
int r = 0; r < num_rhs; r++)
941 rhs(rc,r) = A[j] * sol(r);
947 for (
int r = 0; r < num_rhs; r++)
955 for (
int k = I[col]; 1; k++)
958 mfem_error(
"SparseMatrix::EliminateRowColMultipleRHS() #3");
962 for (
int r = 0; r < num_rhs; r++)
964 rhs(col,r) -= sol(r) * A[k];
971 for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
972 if ((col = aux->Column) == rc)
975 for (
int r = 0; r < num_rhs; r++)
977 rhs(rc,r) = aux->Value * sol(r);
983 for (
int r = 0; r < num_rhs; r++)
991 for (RowNode *node = Rows[col]; 1; node = node->Prev)
994 mfem_error(
"SparseMatrix::EliminateRowColMultipleRHS() #4");
996 else if (node->Column == rc)
998 for (
int r = 0; r < num_rhs; r++)
1000 rhs(col,r) -= sol(r) * node->Value;
1012 MFEM_ASSERT(rc < height && rc >= 0,
1013 "Row " << rc <<
" not in matrix of height " <<
height);
1017 for (
int j = I[rc]; j < I[rc+1]; j++)
1018 if ((col = J[j]) == rc)
1028 for (
int k = I[col]; 1; k++)
1031 mfem_error(
"SparseMatrix::EliminateRowCol() #2");
1033 else if (J[k] == rc)
1042 RowNode *aux, *node;
1044 for (aux = Rows[rc]; aux != NULL; aux = aux->Prev)
1046 if ((col = aux->Column) == rc)
1056 for (node = Rows[col]; 1; node = node->Prev)
1059 mfem_error(
"SparseMatrix::EliminateRowCol() #3");
1061 else if (node->Column == rc)
1078 for (nd = Rows[rc]; nd != NULL; nd = nd->Prev)
1080 if ((col = nd->Column) == rc)
1084 Ae.
Add(rc, rc, nd->Value - 1.0);
1090 Ae.
Add(rc, col, nd->Value);
1092 for (nd2 = Rows[col]; 1; nd2 = nd2->Prev)
1098 else if (nd2->Column == rc)
1100 Ae.
Add(col, rc, nd2->Value);
1110 for (
int j = I[rc]; j < I[rc+1]; j++)
1111 if ((col = J[j]) == rc)
1115 Ae.
Add(rc, rc, A[j] - 1.0);
1121 Ae.
Add(rc, col, A[j]);
1123 for (
int k = I[col];
true; k++)
1128 else if (J[k] == rc)
1130 Ae.
Add(col, rc, A[k]);
1140 for (
int i = 0; i <
height; i++)
1141 if (I[i+1] == I[i]+1 && fabs(A[I[i]]) < 1e-16)
1152 for (i = 0; i <
height; i++)
1155 for (j = I[i]; j < I[i+1]; j++)
1161 for (j = I[i]; j < I[i+1]; j++)
1177 double sum, *yp = y.
GetData();
1178 const double *xp = x.
GetData();
1182 RowNode *diag_p, *n_p, **R = Rows;
1184 for (i = 0; i < s; i++)
1188 for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
1189 if ((c = n_p->Column) == i)
1195 sum += n_p->Value * yp[c];
1198 if (diag_p != NULL && diag_p->Value != 0.0)
1200 yp[i] = (xp[i] - sum) / diag_p->Value;
1202 else if (xp[i] == sum)
1208 mfem_error(
"SparseMatrix::Gauss_Seidel_forw()");
1214 int j, end, d, *Ip = I, *Jp = J;
1218 for (i = 0; i < s; i++)
1223 for ( ; j < end; j++)
1224 if ((c = Jp[j]) == i)
1230 sum += Ap[j] * yp[c];
1233 if (d >= 0 && Ap[d] != 0.0)
1235 yp[i] = (xp[i] - sum) / Ap[d];
1237 else if (xp[i] == sum)
1243 mfem_error(
"SparseMatrix::Gauss_Seidel_forw(...) #2");
1252 double sum, *yp = y.
GetData();
1253 const double *xp = x.
GetData();
1257 RowNode *diag_p, *n_p, **R = Rows;
1259 for (i =
height-1; i >= 0; i--)
1263 for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
1264 if ((c = n_p->Column) == i)
1270 sum += n_p->Value * yp[c];
1273 if (diag_p != NULL && diag_p->Value != 0.0)
1275 yp[i] = (xp[i] - sum) / diag_p->Value;
1277 else if (xp[i] == sum)
1283 mfem_error(
"SparseMatrix::Gauss_Seidel_back()");
1289 int j, beg, d, *Ip = I, *Jp = J;
1293 for (i =
height-1; i >= 0; i--)
1298 for ( ; j >= beg; j--)
1299 if ((c = Jp[j]) == i)
1305 sum += Ap[j] * yp[c];
1308 if (d >= 0 && Ap[d] != 0.0)
1310 yp[i] = (xp[i] - sum) / Ap[d];
1312 else if (xp[i] == sum)
1318 mfem_error(
"SparseMatrix::Gauss_Seidel_back(...) #2");
1326 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
1329 for (
int i = 0; i <
height; i++)
1333 for (
int j = I[i]; j < I[i+1]; j++)
1341 if (d >= 0 && A[d] != 0.0)
1343 double a = 1.8 * fabs(A[d]) / norm;
1351 mfem_error(
"SparseMatrix::GetJacobiScaling() #2");
1360 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
1362 for (
int i = 0; i <
height; i++)
1366 for (
int j = I[i]; j < I[i+1]; j++)
1374 sum -= A[j] * x0(J[j]);
1377 if (d >= 0 && A[d] != 0.0)
1379 x1(i) = sc * (sum / A[d]) + (1.0 - sc) * x0(i);
1390 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
1392 bool scale = (sc != 1.0);
1393 for (
int i = 0, j = 0; i <
height; i++)
1398 MFEM_VERIFY(j != end,
"Couldn't find diagonal in row. i = " << i
1400 <<
", I[i+1] = " << end );
1403 MFEM_VERIFY(std::abs(A[j]) > 0.0,
"Diagonal " << j <<
" must be nonzero");
1406 x(i) = sc * b(i) / A[j];
1423 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
1425 for (
int i = 0; i <
height; i++)
1427 double resi = b(i), norm = 0.0;
1428 for (
int j = I[i]; j < I[i+1]; j++)
1430 resi -= A[j] * x0(J[j]);
1435 x1(i) = x0(i) + sc * resi / norm;
1439 MFEM_ABORT(
"L1 norm of row " << i <<
" is zero.");
1447 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
1449 for (
int i = 0; i <
height; i++)
1451 double resi = b(i), sum = 0.0;
1452 for (
int j = I[i]; j < I[i+1]; j++)
1454 resi -= A[j] * x0(J[j]);
1459 x1(i) = x0(i) + sc * resi / sum;
1463 MFEM_ABORT(
"sum of row " << i <<
" is zero.");
1471 int i, j, gi, gj, s, t;
1474 for (i = 0; i < rows.
Size(); i++)
1476 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
1479 "Trying to insert a row " << gi <<
" outside the matrix height "
1482 for (j = 0; j < cols.
Size(); j++)
1484 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1486 MFEM_ASSERT(gj <
width,
1487 "Trying to insert a column " << gj <<
" outside the matrix width "
1490 if (skip_zeros && a == 0.0)
1494 if (&rows != &cols || subm(j, i) == 0.0)
1499 if (t < 0) { a = -a; }
1511 if ((gi=i) < 0) { gi = -1-gi, s = -1; }
1514 "Trying to insert a row " << gi <<
" outside the matrix height "
1516 if ((gj=j) < 0) { gj = -1-gj, t = -s; }
1518 MFEM_ASSERT(gj <
width,
1519 "Trying to insert a column " << gj <<
" outside the matrix width "
1521 if (t < 0) { a = -a; }
1530 if ((gi=i) < 0) { gi = -1-gi, s = -1; }
1533 "Trying to insert a row " << gi <<
" outside the matrix height "
1535 if ((gj=j) < 0) { gj = -1-gj, t = -s; }
1537 MFEM_ASSERT(gj <
width,
1538 "Trying to insert a column " << gj <<
" outside the matrix width "
1540 if (t < 0) { a = -a; }
1547 int i, j, gi, gj, s, t;
1550 for (i = 0; i < rows.
Size(); i++)
1552 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
1555 "Trying to insert a row " << gi <<
" outside the matrix height "
1558 for (j = 0; j < cols.
Size(); j++)
1561 if (skip_zeros && a == 0.0)
1565 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1567 MFEM_ASSERT(gj <
width,
1568 "Trying to insert a column " << gj <<
" outside the matrix width "
1570 if (t < 0) { a = -a; }
1582 int i, j, gi, gj, s, t;
1585 for (i = 0; i < rows.
Size(); i++)
1587 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
1590 "Trying to insert a row " << gi <<
" outside the matrix height "
1593 for (j = 0; j < cols.
Size(); j++)
1596 if (skip_zeros && a == 0.0)
1600 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1602 MFEM_ASSERT(gj <
width,
1603 "Trying to insert a column " << gj <<
" outside the matrix width "
1605 if (t < 0) { a = -a; }
1615 int i, j, gi, gj, s, t;
1618 for (i = 0; i < rows.
Size(); i++)
1620 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
1623 "Trying to insert a row " << gi <<
" outside the matrix height "
1626 for (j = 0; j < cols.
Size(); j++)
1628 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1630 MFEM_ASSERT(gj <
width,
1631 "Trying to insert a column " << gj <<
" outside the matrix width "
1634 subm(i, j) = (t < 0) ? (-a) : (a);
1649 "Trying to insert a row " << gi <<
" outside the matrix height "
1653 return (Rows[gi] == NULL);
1657 return (I[gi] == I[gi+1]);
1666 if ((gi=row) < 0) { gi = -1-gi; }
1668 "Trying to insert a row " << gi <<
" outside the matrix height "
1672 for (n = Rows[gi], j = 0; n; n = n->Prev)
1678 for (n = Rows[gi], j = 0; n; n = n->Prev, j++)
1680 cols[j] = n->Column;
1693 cols.
MakeRef(J + j, I[gi+1]-j);
1695 MFEM_ASSERT(row >= 0,
"Row not valid: " << row );
1703 int j, gi, gj, s, t;
1706 MFEM_VERIFY(!
Finalized(),
"Matrix must NOT be finalized.");
1708 if ((gi=row) < 0) { gi = -1-gi, s = -1; }
1711 "Trying to insert a row " << gi <<
" outside the matrix height "
1714 for (j = 0; j < cols.
Size(); j++)
1716 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1718 MFEM_ASSERT(gj <
width,
1719 "Trying to insert a column " << gj <<
" outside the matrix width "
1722 if (t < 0) { a = -a; }
1731 int j, gi, gj, s, t;
1734 MFEM_VERIFY(!
Finalized(),
"Matrix must NOT be finalized.");
1736 if ((gi=row) < 0) { gi = -1-gi, s = -1; }
1739 "Trying to insert a row " << gi <<
" outside the matrix height "
1742 for (j = 0; j < cols.
Size(); j++)
1744 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
1746 MFEM_ASSERT(gj <
width,
1747 "Trying to insert a column " << gj <<
" outside the matrix width "
1754 if (t < 0) { a = -a; }
1772 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
1774 aux -> Value *= scale;
1779 int j, end = I[i+1];
1781 for (j = I[i]; j < end; j++)
1794 for (
int i=0; i <
height; ++i)
1797 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
1799 aux -> Value *= scale;
1807 for (
int i=0; i <
height; ++i)
1811 for (j = I[i]; j < end; j++)
1824 for (
int i=0; i <
height; ++i)
1826 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
1828 aux -> Value *= sr(aux->Column);
1836 for (
int i=0; i <
height; ++i)
1839 for (j = I[i]; j < end; j++)
1850 "Mismatch of this matrix size and rhs. This height = "
1851 <<
height <<
", width = " <<
width <<
", B.height = "
1854 for (
int i = 0; i <
height; i++)
1859 for (RowNode *aux = B.Rows[i]; aux != NULL; aux = aux->Prev)
1861 _Add_(aux->Column, aux->Value);
1866 for (
int j = B.I[i]; j < B.I[i+1]; j++)
1868 _Add_(B.J[j], B.A[j]);
1879 for (
int i = 0; i <
height; i++)
1884 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
1886 np->Value += a * B._Get_(np->Column);
1891 for (
int j = I[i]; j < I[i+1]; j++)
1893 A[j] += a * B._Get_(J[j]);
1903 for (
int i = 0, nnz = I[
height]; i < nnz; i++)
1908 for (
int i = 0; i <
height; i++)
1909 for (RowNode *node_p = Rows[i]; node_p != NULL;
1910 node_p = node_p -> Prev)
1912 node_p -> Value = a;
1921 for (
int i = 0, nnz = I[
height]; i < nnz; i++)
1926 for (
int i = 0; i <
height; i++)
1927 for (RowNode *node_p = Rows[i]; node_p != NULL;
1928 node_p = node_p -> Prev)
1930 node_p -> Value *= a;
1943 for (i = 0; i <
height; i++)
1945 out <<
"[row " << i <<
"]\n";
1946 for (nd = Rows[i], j = 0; nd != NULL; nd = nd->Prev, j++)
1948 out <<
" (" << nd->Column <<
"," << nd->Value <<
")";
1949 if ( !((j+1) % _width) )
1962 for (i = 0; i <
height; i++)
1964 out <<
"[row " << i <<
"]\n";
1965 for (j = I[i]; j < I[i+1]; j++)
1967 out <<
" (" << J[j] <<
"," << A[j] <<
")";
1968 if ( !((j+1-I[i]) % _width) )
1973 if ((j-I[i]) % _width)
1982 out <<
"% size " <<
height <<
" " <<
width <<
"\n";
1985 ios::fmtflags old_fmt = out.flags();
1986 out.setf(ios::scientific);
1987 std::streamsize old_prec = out.precision(14);
1989 for (i = 0; i <
height; i++)
1990 for (j = I[i]; j < I[i+1]; j++)
1992 out << i+1 <<
" " << J[j]+1 <<
" " << A[j] <<
'\n';
1994 out.precision(old_prec);
2001 ios::fmtflags old_fmt = out.flags();
2002 out.setf(ios::scientific);
2003 std::streamsize old_prec = out.precision(14);
2005 out <<
"%%MatrixMarket matrix coordinate real general" <<
'\n'
2006 <<
"% Generated by MFEM" <<
'\n';
2009 for (i = 0; i <
height; i++)
2010 for (j = I[i]; j < I[i+1]; j++)
2012 out << i+1 <<
" " << J[j]+1 <<
" " << A[j] <<
'\n';
2014 out.precision(old_prec);
2020 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
2026 for (i = 0; i <=
height; i++)
2028 out << I[i]+1 <<
'\n';
2031 for (i = 0; i < I[
height]; i++)
2033 out << J[i]+1 <<
'\n';
2036 for (i = 0; i < I[
height]; i++)
2038 out << A[i] <<
'\n';
2044 MFEM_VERIFY(
Finalized(),
"Matrix must be finalized.");
2049 out <<
width <<
'\n';
2051 for (i = 0; i <=
height; i++)
2053 out << I[i] <<
'\n';
2056 for (i = 0; i < I[
height]; i++)
2058 out << J[i] <<
'\n';
2061 for (i = 0; i < I[
height]; i++)
2063 out << A[i] <<
'\n';
2069 if ( I != NULL && ownGraph )
2073 if ( J != NULL && ownGraph )
2077 if (A != NULL && ownData )
2084 #if !defined(MFEM_USE_MEMALLOC)
2085 for (
int i = 0; i <
height; i++)
2087 RowNode *aux, *node_p = Rows[i];
2088 while (node_p != NULL)
2091 node_p = node_p->Prev;
2099 if (ColPtrJ != NULL)
2103 if (ColPtrNode != NULL)
2105 delete [] ColPtrNode;
2107 #ifdef MFEM_USE_MEMALLOC
2108 if (NodesMem != NULL)
2121 int * end_j(J + I[
height]);
2122 for (
int * jptr(start_j); jptr != end_j; ++jptr)
2124 awidth = (*jptr > awidth) ? *jptr : awidth;
2130 for (
int i = 0; i <
height; i++)
2131 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
2133 awidth =(aux->Column > awidth) ? aux->Column : awidth;
2147 for (
int i = 0; i < n; i++)
2157 "Finalize must be called before Transpose. Use TransposeRowMatrix instead");
2160 int m, n, nnz, *A_i, *A_j, *At_i, *At_j;
2161 double *A_data, *At_data;
2170 At_i =
new int[n+1];
2171 At_j =
new int[nnz];
2172 At_data =
new double[nnz];
2174 for (i = 0; i <= n; i++)
2178 for (i = 0; i < nnz; i++)
2182 for (i = 1; i < n; i++)
2184 At_i[i+1] += At_i[i];
2187 for (i = j = 0; i < m; i++)
2190 for ( ; j < end; j++)
2192 At_j[At_i[A_j[j]]] = i;
2193 At_data[At_i[A_j[j]]] = A_data[j];
2198 for (i = n; i > 0; i--)
2200 At_i[i] = At_i[i-1];
2211 int m, n, nnz, *At_i, *At_j;
2221 for (i = 0; i < m; i++)
2223 A.
GetRow(i, Acols, Avals);
2241 At_i =
new int[n+1];
2242 At_j =
new int[nnz];
2243 At_data =
new double[nnz];
2245 for (i = 0; i <= n; i++)
2250 for (i = 0; i < m; i++)
2252 A.
GetRow(i, Acols, Avals);
2253 for (j = 0; j<Acols.
Size(); ++j)
2258 for (i = 1; i < n; i++)
2260 At_i[i+1] += At_i[i];
2263 for (i = 0; i < m; i++)
2265 A.
GetRow(i, Acols, Avals);
2266 for (j = 0; j<Acols.
Size(); ++j)
2268 At_j[At_i[Acols[j]]] = i;
2269 At_data[At_i[Acols[j]]] = Avals[j];
2274 for (i = n; i > 0; i--)
2276 At_i[i] = At_i[i-1];
2287 int nrowsA, ncolsA, nrowsB, ncolsB;
2288 int *A_i, *A_j, *B_i, *B_j, *C_i, *C_j, *B_marker;
2289 double *A_data, *B_data, *C_data;
2290 int ia, ib, ic, ja, jb, num_nonzeros;
2291 int row_start, counter;
2292 double a_entry, b_entry;
2300 MFEM_VERIFY(ncolsA == nrowsB,
2301 "number of columns of A (" << ncolsA
2302 <<
") must equal number of rows of B (" << nrowsB <<
")");
2311 B_marker =
new int[ncolsB];
2313 for (ib = 0; ib < ncolsB; ib++)
2320 C_i =
new int[nrowsA+1];
2322 C_i[0] = num_nonzeros = 0;
2323 for (ic = 0; ic < nrowsA; ic++)
2325 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
2328 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
2331 if (B_marker[jb] != ic)
2338 C_i[ic+1] = num_nonzeros;
2341 C_j =
new int[num_nonzeros];
2342 C_data =
new double[num_nonzeros];
2344 C =
new SparseMatrix (C_i, C_j, C_data, nrowsA, ncolsB);
2346 for (ib = 0; ib < ncolsB; ib++)
2355 MFEM_VERIFY(nrowsA == C -> Height() && ncolsB == C -> Width(),
2356 "Input matrix sizes do not match output sizes"
2357 <<
" nrowsA = " << nrowsA
2358 <<
", C->Height() = " << C->
Height()
2359 <<
" ncolsB = " << ncolsB
2360 <<
", C->Width() = " << C->
Width());
2364 C_data = C -> GetData();
2368 for (ic = 0; ic < nrowsA; ic++)
2371 row_start = counter;
2372 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
2375 a_entry = A_data[ia];
2376 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
2379 b_entry = B_data[ib];
2380 if (B_marker[jb] < row_start)
2382 B_marker[jb] = counter;
2387 C_data[counter] = a_entry*b_entry;
2392 C_data[B_marker[jb]] += a_entry*b_entry;
2400 "With pre-allocated output matrix, number of non-zeros ("
2402 <<
") did not match number of entries changed from matrix-matrix multiply, "
2413 int nrowsA, ncolsA, nrowsB, ncolsB;
2414 int *C_i, *C_j, *B_marker;
2416 int ia, ib, ic, ja, jb, num_nonzeros;
2417 int row_start, counter;
2418 double a_entry, b_entry;
2426 MFEM_VERIFY(ncolsA == nrowsB,
2427 "number of columns of A (" << ncolsA
2428 <<
") must equal number of rows of B (" << nrowsB <<
")");
2430 B_marker =
new int[ncolsB];
2432 for (ib = 0; ib < ncolsB; ib++)
2437 C_i =
new int[nrowsA+1];
2439 C_i[0] = num_nonzeros = 0;
2443 for (ic = 0; ic < nrowsA; ic++)
2445 A.
GetRow(ic, colsA, dataA);
2446 for (ia = 0; ia < colsA.
Size(); ia++)
2449 B.
GetRow(ja, colsB, dataB);
2450 for (ib = 0; ib < colsB.
Size(); ib++)
2453 if (B_marker[jb] != ic)
2460 C_i[ic+1] = num_nonzeros;
2463 C_j =
new int[num_nonzeros];
2464 C_data =
new double[num_nonzeros];
2466 C =
new SparseMatrix(C_i, C_j, C_data, nrowsA, ncolsB);
2468 for (ib = 0; ib < ncolsB; ib++)
2474 for (ic = 0; ic < nrowsA; ic++)
2476 row_start = counter;
2477 A.
GetRow(ic, colsA, dataA);
2478 for (ia = 0; ia < colsA.
Size(); ia++)
2481 a_entry = dataA[ia];
2482 B.
GetRow(ja, colsB, dataB);
2483 for (ib = 0; ib < colsB.
Size(); ib++)
2486 b_entry = dataB[ib];
2487 if (B_marker[jb] < row_start)
2489 B_marker[jb] = counter;
2491 C_data[counter] = a_entry*b_entry;
2496 C_data[B_marker[jb]] += a_entry*b_entry;
2532 int i, At_nnz, *At_j;
2536 At_nnz = At -> NumNonZeroElems();
2537 At_j = At -> GetJ();
2538 At_data = At -> GetData();
2539 for (i = 0; i < At_nnz; i++)
2541 At_data[i] *= D(At_j[i]);
2552 int ncols = A.
Width();
2554 int * C_i =
new int[nrows+1];
2558 int * A_i = A.
GetI();
2559 int * A_j = A.
GetJ();
2560 double * A_data = A.
GetData();
2562 int * B_i = B.
GetI();
2563 int * B_j = B.
GetJ();
2564 double * B_data = B.
GetData();
2566 int * marker =
new int[ncols];
2567 std::fill(marker, marker+ncols, -1);
2569 int num_nonzeros = 0, jcol;
2571 for (
int ic = 0; ic < nrows; ic++)
2573 for (
int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
2579 for (
int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
2582 if (marker[jcol] != ic)
2588 C_i[ic+1] = num_nonzeros;
2591 C_j =
new int[num_nonzeros];
2592 C_data =
new double[num_nonzeros];
2594 for (
int ia = 0; ia < ncols; ia++)
2600 for (
int ic = 0; ic < nrows; ic++)
2602 for (
int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
2606 C_data[pos] = a*A_data[ia];
2610 for (
int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
2613 if (marker[jcol] < C_i[ic])
2616 C_data[pos] = b*B_data[ib];
2622 C_data[marker[jcol]] += b*B_data[ib];
2628 return new SparseMatrix(C_i, C_j, C_data, nrows, ncols);
2633 return Add(1.,A,1.,B);
2638 MFEM_ASSERT(Ai.
Size() > 0,
"invalid size Ai.Size() = " << Ai.
Size());
2643 for (
int i=1; i < Ai.
Size(); ++i)
2645 result =
Add(*accumulate, *Ai[i]);
2651 accumulate = result;
2664 Swap(A.Rows, B.Rows);
2665 Swap(A.current_row, B.current_row);
2666 Swap(A.ColPtrJ, B.ColPtrJ);
2668 #ifdef MFEM_USE_MEMALLOC
2669 Swap(A.NodesMem, B.NodesMem);
double GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
int Size() const
Logical size of the array.
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
int RowSize(const int i) const
Returns the number of elements in row i.
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
void Add(const int i, const int j, const double a)
void NewDataAndSize(double *d, int s)
void Print(std::ostream &out=std::cout, int width_=4) const
Prints matrix to stream out.
void DiagScale(const Vector &b, Vector &x, double sc=1.0) const
void SetSize(int s)
Resizes the vector if the new size is different.
void PrintMM(std::ostream &out=std::cout) const
Prints matrix in Matrix Market sparse format.
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Abstract data type for sparse matrices.
template void SortPairs< int, double >(Pair< int, double > *, int)
void EliminateRowCol(int rc, const double sol, Vector &rhs, int d=0)
void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
int Size() const
Returns the size of the vector.
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += At * x (default) or y += a * At * x
Abstract data type for matrix inverse.
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)
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
void SortColumnIndices()
Sort the column indices corresponding to each row.
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
bool RowIsEmpty(const int row) const
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
void ScaleRow(const int row, const double scale)
int ActualWidth()
Returns the actual Width of the matrix.
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void EliminateCol(int col)
void ScaleColumns(const Vector &sr)
virtual ~SparseMatrix()
Destroys sparse matrix.
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
void EliminateCols(Array< int > &cols, Vector *x=NULL, Vector *b=NULL)
Eliminate all columns 'i' for which cols[i] != 0.
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
double IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
void ScaleRows(const Vector &sl)
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const =0
virtual void Finalize(int skip_zeros=1)
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
double * GetData() const
Return element data.
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
int * GetI() const
Return the array I.
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void SparseMatrixFunction(SparseMatrix &S, double(*f)(double))
Applies f() to each element of the matrix (after it is finalized).
virtual MatrixInverse * Inverse() const
Returns a pointer to approximation of the matrix inverse.
void Swap(Array< T > &, Array< T > &)
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
SparseMatrix & operator=(double a)
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void mfem_error(const char *msg)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void Set(const int i, const int j, const double a)
void PrintMatlab(std::ostream &out=std::cout) const
Prints matrix in matlab format.
void Gauss_Seidel_back(const Vector &x, Vector &y) const
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
SparseMatrix(int nrows, int ncols=0)
Creates sparse matrix.
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
void GetDiag(Vector &d) const
Returns the Diagonal of A.
SparseMatrix * Mult_AtDA(const SparseMatrix &A, const Vector &D, SparseMatrix *OAtDA)
double GetRowNorml1(int irow) const
For i = irow compute .
double & operator()(int i, int j)
Returns reference to A[i][j].
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
SparseMatrix & operator*=(double a)
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
void MakeRef(T *, int)
Make this Array a reference to a poiter.
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm)
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
int MaxRowSize() const
Returns the maximum number of elements among all rows.
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
SparseMatrix * MultAbstractSparseMatrix(const AbstractSparseMatrix &A, const AbstractSparseMatrix &B)
Matrix product of sparse matrices. A and B do not need to be CSR matrices.
SparseMatrix & operator+=(SparseMatrix &B)
int * GetJ() const
Return the array J.
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
int CountSmallElems(double tol) const
Count the number of entries with |a_ij| < tol.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
void GetRowSums(Vector &x) const
For all i compute .
void Neg()
(*this) = -(*this)