31#if defined(_MSC_VER) && (_MSC_VER < 1800)
33#define copysign _copysign
49 MFEM_ASSERT(m.data,
"invalid source matrix");
51 std::memcpy(data, m.data,
sizeof(
real_t)*hw);
57 MFEM_ASSERT(s >= 0,
"invalid DenseMatrix size: " << s);
67 MFEM_ASSERT(m >= 0 && n >= 0,
68 "invalid DenseMatrix size: " << m <<
" x " << n);
69 const int capacity = m*n;
78 :
Matrix(mat.width, mat.height)
80 MFEM_CONTRACT_VAR(ch);
86 for (
int i = 0; i <
height; i++)
88 for (
int j = 0; j <
width; j++)
90 (*this)(i,j) = mat(j,i);
98 MFEM_ASSERT(h >= 0 && w >= 0,
99 "invalid DenseMatrix size: " << h <<
" x " << w);
133 MFEM_ASSERT(
height == y.
Size(),
"incompatible dimensions");
141 MFEM_ASSERT(
width == x.
Size(),
"incompatible dimensions");
150 "incompatible dimensions");
160 "incompatible dimensions");
164 for (
int i = 0; i < hw; i++)
166 a += data[i] * m.data[i];
176 for (
int col = 0; col <
width; col++)
179 for (
int row = 0; row <
height; row++)
181 y_col += x[row]*d_col[row];
190 MFEM_ASSERT(
width == y.
Size(),
"incompatible dimensions");
198 MFEM_ASSERT(
height == x.
Size(),
"incompatible dimensions");
207 "incompatible dimensions");
222 "incompatible dimensions");
226 for (
int col = 0; col <
width; col++)
229 for (
int row = 0; row <
height; row++)
231 yp[row] += x_col*d_col[row];
246 "incompatible dimensions");
248 const real_t *d_col = data;
249 for (
int col = 0; col <
width; col++)
252 for (
int row = 0; row <
height; row++)
254 y_col += x[row]*d_col[row];
264 "incompatible dimensions");
271 for (
int col = 0; col <
width; col++)
273 const real_t x_col =
a*xp[col];
274 for (
int row = 0; row <
height; row++)
276 yp[row] += x_col*d_col[row];
286 "incompatible dimensions");
288 const real_t *d_col = data;
289 for (
int col = 0; col <
width; col++)
292 for (
int row = 0; row <
height; row++)
294 y_col += x[row]*d_col[row];
305 for (
int i = 0; i <
height; i++)
308 for (
int j = 0; j <
width; j++)
310 Axi += (*this)(i,j) * x[j];
322 for (
int j = 0; j <
width; ++j)
324 for (
int i = 0; i <
height; ++i)
326 *(it_data++) *= s(i);
335 for (
int j = 0; j <
width; ++j)
337 for (
int i = 0; i <
height; ++i)
339 *(it_data++) /= s(i);
349 for (
int j = 0; j <
width; ++j)
352 for (
int i = 0; i <
height; ++i)
363 for (
int j = 0; j <
width; ++j)
365 const real_t sj = 1./s(j);
366 for (
int i = 0; i <
height; ++i)
378 mfem_error(
"DenseMatrix::SymmetricScaling: dimension mismatch");
384 for (
real_t * end_s = it_s +
width; it_s != end_s; ++it_s)
386 *(it_ss++) = sqrt(*it_s);
390 for (
int j = 0; j <
width; ++j)
392 for (
int i = 0; i <
height; ++i)
394 *(it_data++) *= ss[i]*ss[j];
406 mfem_error(
"DenseMatrix::InvSymmetricScaling: dimension mismatch");
412 for (
real_t * end_s = it_s +
width; it_s != end_s; ++it_s)
414 *(it_ss++) = 1./sqrt(*it_s);
418 for (
int j = 0; j <
width; ++j)
420 for (
int i = 0; i <
height; ++i)
422 *(it_data++) *= ss[i]*ss[j];
434 mfem_error(
"DenseMatrix::Trace() : not a square matrix!");
440 for (
int i = 0; i <
width; i++)
456 "The matrix must be square and "
457 <<
"of size less than or equal to 2."
458 <<
" Height() = " <<
Height()
459 <<
", Width() = " <<
Width());
465 data[0] = std::exp(data[0]);
476 const real_t e = (
a - d)*(
a - d) + 4*
b*c;
477 const real_t f = std::exp((
a + d)/2.0);
478 const real_t g = std::sqrt(std::abs(e)) / 2.0;
482 data[0] = 1.0 + (
a - d)/2.0;
483 data[3] = 1.0 - (
a - d)/2.0;
487 data[0] = std::cosh(g) + (
a - d)/2 * std::sinh(g) / g;
488 data[1] =
b * std::sinh(g) / g;
489 data[2] = c * std::sinh(g) / g;
490 data[3] = std::cosh(g) - (
a - d)/2 * std::sinh(g) / g;
494 data[0] = std::cos(g) + (
a - d)/2 * std::sin(g) / g;
495 data[1] =
b * std::sin(g) / g;
496 data[2] = c * std::sin(g) / g;
497 data[3] = std::cos(g) - (
a - d)/2 * std::sin(g) / g;
499 for (
int i = 0; i < 4; i++)
507 MFEM_ABORT(
"3x3 matrices are not currently supported");
511 MFEM_ABORT(
"Only 1x1 and 2x2 matrices are currently supported");
519 "The matrix must be square and "
520 <<
"sized larger than zero to compute the determinant."
521 <<
" Height() = " <<
Height()
522 <<
", Width() = " <<
Width());
530 return data[0] * data[3] - data[1] * data[2];
536 d[0] * (d[4] * d[8] - d[5] * d[7]) +
537 d[3] * (d[2] * d[7] - d[1] * d[8]) +
538 d[6] * (d[1] * d[5] - d[2] * d[4]);
544 d[ 0] * (d[ 5] * (d[10] * d[15] - d[11] * d[14]) -
545 d[ 9] * (d[ 6] * d[15] - d[ 7] * d[14]) +
546 d[13] * (d[ 6] * d[11] - d[ 7] * d[10])
548 d[ 4] * (d[ 1] * (d[10] * d[15] - d[11] * d[14]) -
549 d[ 9] * (d[ 2] * d[15] - d[ 3] * d[14]) +
550 d[13] * (d[ 2] * d[11] - d[ 3] * d[10])
552 d[ 8] * (d[ 1] * (d[ 6] * d[15] - d[ 7] * d[14]) -
553 d[ 5] * (d[ 2] * d[15] - d[ 3] * d[14]) +
554 d[13] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
556 d[12] * (d[ 1] * (d[ 6] * d[11] - d[ 7] * d[10]) -
557 d[ 5] * (d[ 2] * d[11] - d[ 3] * d[10]) +
558 d[ 9] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
567 return lu_factors.
Det();
582 return sqrt(data[0] * data[0] + data[1] * data[1]);
586 return sqrt(data[0] * data[0] + data[1] * data[1] + data[2] * data[2]);
591 real_t E = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
592 real_t G = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
593 real_t F = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
594 return sqrt(E * G - F * F);
596 mfem_error(
"DenseMatrix::Weight(): mismatched or unsupported dimensions");
603 for (
int i = 0; i < s; i++)
605 data[i] =
alpha*A[i];
611 for (
int j = 0; j <
Width(); j++)
613 for (
int i = 0; i <
Height(); i++)
615 (*this)(i,j) += c * A(i,j);
623 for (
int i = 0; i < s; i++)
632 for (
int i = 0; i < s; i++)
642 for (
int i = 0; i < s; i++)
654 for (
int i = 0; i < hw; i++)
671 "incompatible matrix sizes.");
677 for (
int j = 0; j <
width; j++)
679 for (
int i = 0; i <
height; i++)
681 (*this)(i, j) -= m(i, j);
691 for (
int i = 0; i < s; i++)
701 for (
int i = 0; i < hw; i++)
712 mfem_error(
"DenseMatrix::Invert(): dimension mismatch");
716#ifdef MFEM_USE_LAPACK
717 int *ipiv =
new int[
width];
726 mfem_error(
"DenseMatrix::Invert() : Error in DGETRF");
729 MFEM_LAPACK_PREFIX(getri_)(&
width, data, &
width, ipiv, &qwork, &lwork, &info);
734 MFEM_LAPACK_PREFIX(getri_)(&
width, data, &
width, ipiv, work, &lwork, &info);
738 mfem_error(
"DenseMatrix::Invert() : Error in DGETRI");
744 int c, i, j, n =
Width();
748 for (c = 0; c < n; c++)
750 a = fabs((*
this)(c, c));
752 for (j = c + 1; j < n; j++)
754 b = fabs((*
this)(j, c));
763 mfem_error(
"DenseMatrix::Invert() : singular matrix");
766 for (j = 0; j < n; j++)
771 a = (*this)(c, c) = 1.0 / (*
this)(c, c);
772 for (j = 0; j < c; j++)
776 for (j++; j < n; j++)
780 for (i = 0; i < c; i++)
782 (*this)(i, c) =
a * (
b = -(*
this)(i, c));
783 for (j = 0; j < c; j++)
785 (*this)(i, j) +=
b * (*
this)(c, j);
787 for (j++; j < n; j++)
789 (*this)(i, j) +=
b * (*
this)(c, j);
792 for (i++; i < n; i++)
794 (*this)(i, c) =
a * (
b = -(*
this)(i, c));
795 for (j = 0; j < c; j++)
797 (*this)(i, j) +=
b * (*
this)(c, j);
799 for (j++; j < n; j++)
801 (*this)(i, j) +=
b * (*
this)(c, j);
806 for (c = n - 1; c >= 0; c--)
809 for (i = 0; i < n; i++)
823 mfem_error(
"DenseMatrix::SquareRootInverse() matrix not square.");
833 for (
int v = 0; v <
Height() ; v++) { (*this)(v,v) = 1.0; }
835 for (
int j = 0; j < 10; j++)
837 for (
int i = 0; i < 10; i++)
852 for (
int v = 0; v <
Height() ; v++) { tmp2(v,v) -= 1.0; }
853 if (tmp2.FNorm() < 1e-10) {
break; }
856 if (tmp2.FNorm() > 1e-10)
858 mfem_error(
"DenseMatrix::SquareRootInverse not converged");
864 for (
int j = 0; j <
Width(); j++)
867 for (
int i = 0; i <
Height(); i++)
869 v[j] += (*this)(i,j)*(*
this)(i,j);
879 real_t norm = 0.0, abs_entry;
881 for (
int i = 0; i < hw; i++)
883 abs_entry = fabs(d[i]);
884 if (norm < abs_entry)
896 real_t max_norm = 0.0, entry, fnorm2;
898 for (i = 0; i < hw; i++)
900 entry = fabs(data[i]);
901 if (entry > max_norm)
909 scale_factor = scaled_fnorm2 = 0.0;
914 for (i = 0; i < hw; i++)
916 entry = data[i] / max_norm;
917 fnorm2 += entry * entry;
920 scale_factor = max_norm;
921 scaled_fnorm2 = fnorm2;
926#ifdef MFEM_USE_LAPACK
944 int *ISUPPZ =
new int[2*N];
962 int hw =
a.Height() *
a.Width();
965 for (
int i = 0; i < hw; i++)
970 MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
971 &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, &QWORK,
972 &LWORK, &QIWORK, &LIWORK, &INFO);
978 IWORK =
new int[LIWORK];
980 MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
981 &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, WORK,
982 &LWORK, IWORK, &LIWORK, &INFO);
986 mfem::err <<
"dsyevr_Eigensystem(...): DSYEVR error code: "
994 mfem::err <<
"dsyevr_Eigensystem(...):\n"
995 <<
" DSYEVR did not find all eigenvalues "
996 << M <<
"/" << N << endl;
1001 mfem_error(
"dsyevr_Eigensystem(...): inf/nan values in W");
1005 mfem_error(
"dsyevr_Eigensystem(...): inf/nan values in Z");
1008 for (IL = 0; IL < N; IL++)
1009 for (IU = 0; IU <= IL; IU++)
1012 for (M = 0; M < N; M++)
1014 VL += Z[M+IL*N] * Z[M+IU*N];
1031 <<
" Z^t Z - I deviation = " << VU
1032 <<
"\n W[max] = " << W[N-1] <<
", W[min] = "
1033 << W[0] <<
", N = " << N << endl;
1040 <<
" Z^t Z - I deviation = " << VU
1041 <<
"\n W[max] = " << W[N-1] <<
", W[min] = "
1042 << W[0] <<
", N = " << N << endl;
1046 mfem_error(
"dsyevr_Eigensystem(...): ERROR: ...");
1049 for (IL = 0; IL < N; IL++)
1050 for (IU = 0; IU < N; IU++)
1053 for (M = 0; M < N; M++)
1055 VL += Z[IL+M*N] * W[M] * Z[IU+M*N];
1057 VL = fabs(VL-data[IL+N*IU]);
1066 <<
" max matrix deviation = " << VU
1067 <<
"\n W[max] = " << W[N-1] <<
", W[min] = "
1068 << W[0] <<
", N = " << N << endl;
1072 mfem_error(
"dsyevr_Eigensystem(...): ERROR: ...");
1081 MFEM_CONTRACT_VAR(
a);
1082 MFEM_CONTRACT_VAR(ev);
1083 MFEM_CONTRACT_VAR(evect);
1089#ifdef MFEM_USE_LAPACK
1115 int hw =
a.Height() *
a.Width();
1117 for (
int i = 0; i < hw; i++)
1122 MFEM_LAPACK_PREFIX(syev_)(&JOBZ, &UPLO, &N, A, &LDA, W, &QWORK, &LWORK, &INFO);
1124 LWORK = (int) QWORK;
1125 WORK =
new real_t[LWORK];
1127 MFEM_LAPACK_PREFIX(syev_)(&JOBZ, &UPLO, &N, A, &LDA, W, WORK, &LWORK, &INFO);
1131 mfem::err <<
"dsyev_Eigensystem: DSYEV error code: " << INFO << endl;
1136 if (evect == NULL) {
delete [] A; }
1138 MFEM_CONTRACT_VAR(
a);
1139 MFEM_CONTRACT_VAR(ev);
1140 MFEM_CONTRACT_VAR(evect);
1144void DenseMatrix::Eigensystem(Vector &ev, DenseMatrix *evect)
1146#ifdef MFEM_USE_LAPACK
1154 MFEM_CONTRACT_VAR(ev);
1155 MFEM_CONTRACT_VAR(evect);
1156 mfem_error(
"DenseMatrix::Eigensystem: Compiled without LAPACK");
1164#ifdef MFEM_USE_LAPACK
1193 int hw =
a.Height() *
a.Width();
1196 for (
int i = 0; i < hw; i++)
1202 MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W,
1203 &QWORK, &LWORK, &INFO);
1205 LWORK = (int) QWORK;
1206 WORK =
new real_t[LWORK];
1208 MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK,
1213 mfem::err <<
"dsygv_Eigensystem: DSYGV error code: " << INFO << endl;
1219 if (evect == NULL) {
delete [] A; }
1221 MFEM_CONTRACT_VAR(
a);
1222 MFEM_CONTRACT_VAR(
b);
1223 MFEM_CONTRACT_VAR(ev);
1224 MFEM_CONTRACT_VAR(evect);
1228void DenseMatrix::Eigensystem(DenseMatrix &
b, Vector &ev,
1231#ifdef MFEM_USE_LAPACK
1236 MFEM_CONTRACT_VAR(
b);
1237 MFEM_CONTRACT_VAR(ev);
1238 MFEM_CONTRACT_VAR(evect);
1239 mfem_error(
"DenseMatrix::Eigensystem(generalized): Compiled without LAPACK");
1245#ifdef MFEM_USE_LAPACK
1251 real_t *
a = copy_of_this.data;
1261 MFEM_LAPACK_PREFIX(gesvd_)(&jobu, &jobvt, &m, &n,
a, &m, s,
u, &m, vt, &n,
1262 &qwork, &lwork, &info);
1264 lwork = (int) qwork;
1265 work =
new real_t[lwork];
1267 MFEM_LAPACK_PREFIX(gesvd_)(&jobu, &jobvt, &m, &n,
a, &m, s,
u, &m, vt, &n,
1268 work, &lwork, &info);
1273 mfem::err <<
"DenseMatrix::SingularValues : info = " << info << endl;
1277 MFEM_CONTRACT_VAR(sv);
1279 mfem_error(
"DenseMatrix::SingularValues: Compiled without LAPACK");
1289 for (
int i=0; i < sv.
Size(); ++i)
1301 "The matrix must be square and sized 1, 2, or 3 to compute the"
1303 <<
" Height() = " <<
Height()
1304 <<
", Width() = " <<
Width());
1351 const real_t* rp = data + r;
1354 for (
int i = 0; i < n; i++)
1369 for (
int i = 0; i < m; i++)
1383 for (
int i = 0; i <
height; ++i)
1385 d(i) = (*this)(i,i);
1399 for (
int j = 0; j <
width; ++j)
1400 for (
int i = 0; i <
height; ++i)
1402 l(i) += fabs((*
this)(i,j));
1409 for (
int i = 0; i <
height; i++)
1412 for (
int j = 0; j <
width; j++)
1425 for (
int i = 0; i < N; i++)
1429 for (
int i = 0; i < n; i++)
1440 for (i = 0; i < N; i++)
1444 for (i = 0; i < n; i++)
1446 data[i*(n+1)] = diag[i];
1457 for (i = 0; i <
Height(); i++)
1458 for (j = i+1; j <
Width(); j++)
1461 (*this)(i,j) = (*
this)(j,i);
1476 for (
int i = 0; i <
Height(); i++)
1477 for (
int j = 0; j <
Width(); j++)
1479 (*this)(i,j) = A(j,i);
1488 mfem_error(
"DenseMatrix::Symmetrize() : not a square matrix!");
1496 for (
int i = 0; i <
Height(); i++)
1499 for (
int j = 0; j <
Width(); j++)
1502 (*this)(i, j) = 0.0;
1516 mfem_error(
"DenseMatrix::GradToCurl(...): dimension mismatch");
1522 for (
int i = 0; i < n; i++)
1539 for (
int i = 0; i < n; i++)
1569 MFEM_VERIFY(
Width() == 2,
1570 "DenseMatrix::GradToVectorCurl2D(...): dimension must be 2")
1574 for (
int i = 0; i < n; i++)
1576 curl(i,0) = (*this)(i,1);
1577 curl(i,1) = -(*this)(i,0);
1583 MFEM_ASSERT(
Width()*
Height() == div.
Size(),
"incompatible Vector 'div'!");
1590 for (
int i = 0; i < n; i++)
1600 for (
int j = 0; j <
Width(); j++)
1602 for (
int i = row1; i <= row2; i++)
1604 (*this)(i-row1,j) = A(i,j);
1613 for (
int j = col1; j <= col2; j++)
1615 for (
int i = 0; i <
Height(); i++)
1617 (*this)(i,j-col1) = A(i,j);
1626 for (
int j = 0; j < n; j++)
1628 for (
int i = 0; i < m; i++)
1630 (*this)(i,j) = A(Aro+i,Aco+j);
1639 for (
int j = 0; j < A.
Width(); j++)
1641 for (
int i = 0; i < A.
Height(); i++)
1643 (*this)(row_offset+i,col_offset+j) = *(v++);
1652 for (
int i = 0; i < A.
Width(); i++)
1654 for (
int j = 0; j < A.
Height(); j++)
1656 (*this)(row_offset+i,col_offset+j) = *(v++);
1662 int row_offset,
int col_offset)
1664 MFEM_VERIFY(row_offset+m <= this->
Height() && col_offset+n <= this->
Width(),
1665 "this DenseMatrix is too small to accommodate the submatrix. "
1666 <<
"row_offset = " << row_offset
1668 <<
", this->Height() = " << this->
Height()
1669 <<
", col_offset = " << col_offset
1671 <<
", this->Width() = " << this->
Width()
1673 MFEM_VERIFY(Aro+m <= A.
Height() && Aco+n <= A.
Width(),
1674 "The A DenseMatrix is too small to accommodate the submatrix. "
1677 <<
", A.Height() = " << A.
Height()
1678 <<
", Aco = " << Aco
1680 <<
", A.Width() = " << A.
Width()
1683 for (
int j = 0; j < n; j++)
1685 for (
int i = 0; i < m; i++)
1687 (*this)(row_offset+i,col_offset+j) = A(Aro+i,Aco+j);
1694 for (
int i = 0; i < n; i++)
1696 for (
int j = i+1; j < n; j++)
1698 (*this)(row_offset+i,col_offset+j) =
1699 (*
this)(row_offset+j,col_offset+i) = 0.0;
1703 for (
int i = 0; i < n; i++)
1705 (*this)(row_offset+i,col_offset+i) = c;
1712 for (
int i = 0; i < n; i++)
1714 for (
int j = i+1; j < n; j++)
1716 (*this)(row_offset+i,col_offset+j) =
1717 (*
this)(row_offset+j,col_offset+i) = 0.0;
1721 for (
int i = 0; i < n; i++)
1723 (*this)(row_offset+i,col_offset+i) = diag[i];
1731 int i, j, i_off = 0, j_off = 0;
1733 for (j = 0; j < A.
Width(); j++)
1740 for (i = 0; i < A.
Height(); i++)
1747 (*this)(i-i_off,j-j_off) = A(i,j);
1763 if (co+aw >
Width() || ro+ah > h)
1765 mfem_error(
"DenseMatrix::AddMatrix(...) 1 : dimension mismatch");
1769 p = data + ro + co * h;
1772 for (
int c = 0; c < aw; c++)
1774 for (
int r = 0; r < ah; r++)
1793 if (co+aw >
Width() || ro+ah > h)
1795 mfem_error(
"DenseMatrix::AddMatrix(...) 2 : dimension mismatch");
1799 p = data + ro + co * h;
1802 for (
int c = 0; c < aw; c++)
1804 for (
int r = 0; r < ah; r++)
1816 int idx_max = idx.
Max();
1817 MFEM_VERIFY(idx.
Min() >=0 && idx_max < this->
height && idx_max < this->
width,
1818 "DenseMatrix::GetSubMatrix: Index out of bounds");
1823 for (
int i = 0; i<k; i++)
1826 for (
int j = 0; j<k; j++)
1829 adata[i+j*k] = this->data[ii+jj*
height];
1837 int k = idx_i.
Size();
1838 int l = idx_j.
Size();
1840 MFEM_VERIFY(idx_i.
Min() >=0 && idx_i.
Max() < this->height,
1841 "DenseMatrix::GetSubMatrix: Row index out of bounds");
1842 MFEM_VERIFY(idx_j.
Min() >=0 && idx_j.
Max() < this->width,
1843 "DenseMatrix::GetSubMatrix: Col index out of bounds");
1849 for (
int i = 0; i<k; i++)
1852 for (
int j = 0; j<l; j++)
1855 adata[i+j*k] = this->data[ii+jj*
height];
1862 MFEM_VERIFY(iend >= ibeg,
"DenseMatrix::GetSubMatrix: Inconsistent range");
1863 MFEM_VERIFY(ibeg >=0,
1864 "DenseMatrix::GetSubMatrix: Negative index");
1865 MFEM_VERIFY(iend <= this->
height && iend <= this->
width,
1866 "DenseMatrix::GetSubMatrix: Index bigger than upper bound");
1868 int k = iend - ibeg;
1873 for (
int i = 0; i<k; i++)
1876 for (
int j = 0; j<k; j++)
1879 adata[i+j*k] = this->data[ii+jj*
height];
1887 MFEM_VERIFY(iend >= ibeg,
1888 "DenseMatrix::GetSubMatrix: Inconsistent row range");
1889 MFEM_VERIFY(jend >= jbeg,
1890 "DenseMatrix::GetSubMatrix: Inconsistent col range");
1891 MFEM_VERIFY(ibeg >=0,
1892 "DenseMatrix::GetSubMatrix: Negative row index");
1893 MFEM_VERIFY(jbeg >=0,
1894 "DenseMatrix::GetSubMatrix: Negative row index");
1895 MFEM_VERIFY(iend <= this->
height,
1896 "DenseMatrix::GetSubMatrix: Index bigger than row upper bound");
1897 MFEM_VERIFY(jend <= this->
width,
1898 "DenseMatrix::GetSubMatrix: Index bigger than col upper bound");
1900 int k = iend - ibeg;
1901 int l = jend - jbeg;
1906 for (
int i = 0; i<k; i++)
1909 for (
int j = 0; j<l; j++)
1912 adata[i+j*k] = this->data[ii+jj*
height];
1921 "DenseMatrix::SetSubMatrix:Inconsistent matrix dimensions");
1923 int idx_max = idx.
Max();
1925 MFEM_VERIFY(idx.
Min() >=0,
1926 "DenseMatrix::SetSubMatrix: Negative index");
1927 MFEM_VERIFY(idx_max < this->
height,
1928 "DenseMatrix::SetSubMatrix: Index bigger than row upper bound");
1929 MFEM_VERIFY(idx_max < this->
width,
1930 "DenseMatrix::SetSubMatrix: Index bigger than col upper bound");
1935 for (
int i = 0; i<k; i++)
1938 for (
int j = 0; j<k; j++)
1941 this->data[ii+jj*
height] = adata[i+j*k];
1949 int k = idx_i.
Size();
1950 int l = idx_j.
Size();
1952 "DenseMatrix::SetSubMatrix:Inconsistent matrix dimensions");
1953 MFEM_VERIFY(idx_i.
Min() >=0,
1954 "DenseMatrix::SetSubMatrix: Negative row index");
1955 MFEM_VERIFY(idx_j.
Min() >=0,
1956 "DenseMatrix::SetSubMatrix: Negative col index");
1957 MFEM_VERIFY(idx_i.
Max() < this->height,
1958 "DenseMatrix::SetSubMatrix: Index bigger than row upper bound");
1959 MFEM_VERIFY(idx_j.
Max() < this->width,
1960 "DenseMatrix::SetSubMatrix: Index bigger than col upper bound");
1965 for (
int i = 0; i<k; i++)
1968 for (
int j = 0; j<l; j++)
1971 this->data[ii+jj*
height] = adata[i+j*k];
1980 MFEM_VERIFY(A.
Width() == k,
"DenseMatrix::SetSubmatrix: A is not square");
1981 MFEM_VERIFY(ibeg >=0,
1982 "DenseMatrix::SetSubmatrix: Negative index");
1983 MFEM_VERIFY(ibeg + k <= this->
height,
1984 "DenseMatrix::SetSubmatrix: index bigger than row upper bound");
1985 MFEM_VERIFY(ibeg + k <= this->
width,
1986 "DenseMatrix::SetSubmatrix: index bigger than col upper bound");
1991 for (
int i = 0; i<k; i++)
1994 for (
int j = 0; j<k; j++)
1997 this->data[ii+jj*
height] = adata[i+j*k];
2007 MFEM_VERIFY(ibeg>=0,
2008 "DenseMatrix::SetSubmatrix: Negative row index");
2009 MFEM_VERIFY(jbeg>=0,
2010 "DenseMatrix::SetSubmatrix: Negative col index");
2011 MFEM_VERIFY(ibeg + k <= this->
height,
2012 "DenseMatrix::SetSubmatrix: Index bigger than row upper bound");
2013 MFEM_VERIFY(jbeg + l <= this->
width,
2014 "DenseMatrix::SetSubmatrix: Index bigger than col upper bound");
2019 for (
int i = 0; i<k; i++)
2022 for (
int j = 0; j<l; j++)
2025 this->data[ii+jj*
height] = adata[i+j*k];
2034 "DenseMatrix::AddSubMatrix:Inconsistent matrix dimensions");
2036 int idx_max = idx.
Max();
2038 MFEM_VERIFY(idx.
Min() >=0,
"DenseMatrix::AddSubMatrix: Negative index");
2039 MFEM_VERIFY(idx_max < this->
height,
2040 "DenseMatrix::AddSubMatrix: Index bigger than row upper bound");
2041 MFEM_VERIFY(idx_max < this->
width,
2042 "DenseMatrix::AddSubMatrix: Index bigger than col upper bound");
2047 for (
int i = 0; i<k; i++)
2050 for (
int j = 0; j<k; j++)
2053 this->data[ii+jj*
height] += adata[i+j*k];
2061 int k = idx_i.
Size();
2062 int l = idx_j.
Size();
2064 "DenseMatrix::AddSubMatrix:Inconsistent matrix dimensions");
2066 MFEM_VERIFY(idx_i.
Min() >=0,
2067 "DenseMatrix::AddSubMatrix: Negative row index");
2068 MFEM_VERIFY(idx_j.
Min() >=0,
2069 "DenseMatrix::AddSubMatrix: Negative col index");
2070 MFEM_VERIFY(idx_i.
Max() < this->height,
2071 "DenseMatrix::AddSubMatrix: Index bigger than row upper bound");
2072 MFEM_VERIFY(idx_j.
Max() < this->width,
2073 "DenseMatrix::AddSubMatrix: Index bigger than col upper bound");
2078 for (
int i = 0; i<k; i++)
2081 for (
int j = 0; j<l; j++)
2084 this->data[ii+jj*
height] += adata[i+j*k];
2092 MFEM_VERIFY(A.
Width() == k,
"DenseMatrix::AddSubmatrix: A is not square");
2094 MFEM_VERIFY(ibeg>=0,
2095 "DenseMatrix::AddSubmatrix: Negative index");
2096 MFEM_VERIFY(ibeg + k <= this->
Height(),
2097 "DenseMatrix::AddSubmatrix: Index bigger than row upper bound");
2098 MFEM_VERIFY(ibeg + k <= this->
Width(),
2099 "DenseMatrix::AddSubmatrix: Index bigger than col upper bound");
2104 for (
int i = 0; i<k; i++)
2107 for (
int j = 0; j<k; j++)
2110 this->data[ii+jj*
height] += adata[i+j*k];
2120 MFEM_VERIFY(ibeg>=0,
2121 "DenseMatrix::AddSubmatrix: Negative row index");
2122 MFEM_VERIFY(jbeg>=0,
2123 "DenseMatrix::AddSubmatrix: Negative col index");
2124 MFEM_VERIFY(ibeg + k <= this->
height,
2125 "DenseMatrix::AddSubmatrix: Index bigger than row upper bound");
2126 MFEM_VERIFY(jbeg + l <= this->
width,
2127 "DenseMatrix::AddSubmatrix: Index bigger than col upper bound");
2132 for (
int i = 0; i<k; i++)
2135 for (
int j = 0; j<l; j++)
2138 this->data[ii+jj*
height] += adata[i+j*k];
2148 for (
int i = 0; i < n; i++)
2150 vdata[i] += data[i];
2159 for (
int i = 0; i < n; i++)
2172 mfem_error(
"DenseMatrix::AdjustDofDirection(...): dimension mismatch");
2177 for (
int i = 0; i < n-1; i++)
2179 const int s = (dof[i] < 0) ? (-1) : (1);
2180 for (
int j = i+1; j < n; j++)
2182 const int t = (dof[j] < 0) ? (-s) : (s);
2185 (*this)(i,j) = -(*
this)(i,j);
2186 (*this)(j,i) = -(*
this)(j,i);
2194 for (
int j = 0; j <
Width(); j++)
2196 (*this)(row, j) = value;
2202 for (
int i = 0; i <
Height(); i++)
2204 (*this)(i, col) = value;
2210 MFEM_ASSERT(row !=
nullptr,
"supplied row pointer is null");
2211 for (
int j = 0; j <
Width(); j++)
2213 (*this)(r, j) = row[j];
2225 MFEM_ASSERT(col !=
nullptr,
"supplied column pointer is null");
2226 for (
int i = 0; i <
Height(); i++)
2228 (*this)(i, c) = col[i];
2240 for (
int col = 0; col <
Width(); col++)
2242 for (
int row = 0; row <
Height(); row++)
2244 if (std::abs(
operator()(row,col)) <= eps)
2255 ios::fmtflags old_flags = os.flags();
2257 os << setiosflags(ios::scientific | ios::showpos);
2258 for (
int i = 0; i <
height; i++)
2260 os <<
"[row " << i <<
"]\n";
2261 for (
int j = 0; j <
width; j++)
2264 if (j+1 ==
width || (j+1) % width_ == 0)
2275 os.flags(old_flags);
2281 ios::fmtflags old_flags = os.flags();
2283 os << setiosflags(ios::scientific | ios::showpos);
2284 for (
int i = 0; i <
height; i++)
2286 for (
int j = 0; j <
width; j++)
2294 os.flags(old_flags);
2299 ios::fmtflags old_fmt = os.flags();
2300 os.setf(ios::scientific);
2301 std::streamsize old_prec = os.precision(14);
2303 os <<
"(* Read file into Mathematica using: "
2304 <<
"myMat = Get[\"this_file_name\"] *)\n";
2307 for (
int i = 0; i <
height; i++)
2310 for (
int j = 0; j <
width; j++)
2312 os <<
"Internal`StringToMReal[\"" << (*this)(i,j) <<
"\"]";
2313 if (j <
width - 1) { os <<
','; }
2317 if (i <
height - 1) { os <<
','; }
2322 os.precision(old_prec);
2329 ios::fmtflags old_flags = os.flags();
2331 os << setiosflags(ios::scientific | ios::showpos);
2332 for (
int j = 0; j <
width; j++)
2334 os <<
"[col " << j <<
"]\n";
2335 for (
int i = 0; i <
height; i++)
2338 if (i+1 ==
height || (i+1) % width_ == 0)
2349 os.flags(old_flags);
2358 for (
int i = 0; i <
width; i++)
2363 <<
", cond_F = " <<
FNorm()*copy.FNorm() << endl;
2404 MFEM_VERIFY(A.
IsSquare(),
"A must be a square matrix!");
2405 MFEM_ASSERT(A.
NumCols() > 0,
"supplied matrix, A, is empty!");
2406 MFEM_ASSERT(X !=
nullptr,
"supplied vector, X, is null!");
2415 if (std::abs(det) <= TOL) {
return false; }
2423 if (std::abs(det) <= TOL) {
return false; }
2425 real_t invdet = 1. / det;
2430 X[0] = ( A(1,1)*b0 - A(0,1)*b1) * invdet;
2431 X[1] = (-A(1,0)*b0 + A(0,0)*b1) * invdet;
2440 if (!lu.
Factor(N,TOL)) {
return false; }
2452 MFEM_ASSERT(
a.Height() ==
b.Height() &&
a.Width() == c.
Width() &&
2453 b.Width() == c.
Height(),
"incompatible dimensions");
2455#ifdef MFEM_USE_LAPACK
2456 static char transa =
'N', transb =
'N';
2458 int m =
b.Height(), n = c.
Width(), k =
b.Width();
2460 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha,
b.Data(), &m,
2463 const int ah =
a.Height();
2464 const int aw =
a.Width();
2465 const int bw =
b.Width();
2476 MFEM_ASSERT(
a.Height() ==
b.Height() &&
a.Width() == c.
Width() &&
2477 b.Width() == c.
Height(),
"incompatible dimensions");
2479#ifdef MFEM_USE_LAPACK
2480 static char transa =
'N', transb =
'N';
2482 int m =
b.Height(), n = c.
Width(), k =
b.Width();
2484 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha,
b.Data(), &m,
2487 const int ah =
a.Height();
2488 const int aw =
a.Width();
2489 const int bw =
b.Width();
2493 for (
int j = 0; j < aw; j++)
2495 for (
int k = 0; k < bw; k++)
2497 for (
int i = 0; i < ah; i++)
2499 ad[i+j*ah] +=
alpha * bd[i+k*ah] * cd[k+j*bw];
2508 MFEM_ASSERT(
a.Height() ==
b.Height() &&
a.Width() == c.
Width() &&
2509 b.Width() == c.
Height(),
"incompatible dimensions");
2511#ifdef MFEM_USE_LAPACK
2512 static char transa =
'N', transb =
'N';
2514 int m =
b.Height(), n = c.
Width(), k =
b.Width();
2516 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha,
b.Data(), &m,
2519 const int ah =
a.Height();
2520 const int aw =
a.Width();
2521 const int bw =
b.Width();
2525 for (
int j = 0; j < aw; j++)
2527 for (
int k = 0; k < bw; k++)
2529 for (
int i = 0; i < ah; i++)
2531 ad[i+j*ah] += bd[i+k*ah] * cd[k+j*bw];
2541 if (
a.Width() >
a.Height() ||
a.Width() < 1 ||
a.Height() > 3)
2543 mfem_error(
"CalcAdjugate(...): unsupported dimensions");
2545 if (
a.Width() != adja.
Height() ||
a.Height() != adja.
Width())
2547 mfem_error(
"CalcAdjugate(...): dimension mismatch");
2551 if (
a.Width() <
a.Height())
2560 if (
a.Height() == 3)
2569 e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
2570 g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
2571 f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
2573 ad[0] = d[0]*g - d[3]*
f;
2574 ad[1] = d[3]*e - d[0]*
f;
2575 ad[2] = d[1]*g - d[4]*
f;
2576 ad[3] = d[4]*e - d[1]*
f;
2577 ad[4] = d[2]*g - d[5]*
f;
2578 ad[5] = d[5]*e - d[2]*
f;
2587 else if (
a.Width() == 2)
2590 adja(0,1) = -
a(0,1);
2591 adja(1,0) = -
a(1,0);
2596 adja(0,0) =
a(1,1)*
a(2,2)-
a(1,2)*
a(2,1);
2597 adja(0,1) =
a(0,2)*
a(2,1)-
a(0,1)*
a(2,2);
2598 adja(0,2) =
a(0,1)*
a(1,2)-
a(0,2)*
a(1,1);
2600 adja(1,0) =
a(1,2)*
a(2,0)-
a(1,0)*
a(2,2);
2601 adja(1,1) =
a(0,0)*
a(2,2)-
a(0,2)*
a(2,0);
2602 adja(1,2) =
a(0,2)*
a(1,0)-
a(0,0)*
a(1,2);
2604 adja(2,0) =
a(1,0)*
a(2,1)-
a(1,1)*
a(2,0);
2605 adja(2,1) =
a(0,1)*
a(2,0)-
a(0,0)*
a(2,1);
2606 adja(2,2) =
a(0,0)*
a(1,1)-
a(0,1)*
a(1,0);
2613 if (
a.Height() !=
a.Width() || adjat.
Height() != adjat.
Width() ||
2614 a.Width() != adjat.
Width() ||
a.Width() < 1 ||
a.Width() > 3)
2616 mfem_error(
"CalcAdjugateTranspose(...): dimension mismatch");
2623 else if (
a.Width() == 2)
2625 adjat(0,0) =
a(1,1);
2626 adjat(1,0) = -
a(0,1);
2627 adjat(0,1) = -
a(1,0);
2628 adjat(1,1) =
a(0,0);
2632 adjat(0,0) =
a(1,1)*
a(2,2)-
a(1,2)*
a(2,1);
2633 adjat(1,0) =
a(0,2)*
a(2,1)-
a(0,1)*
a(2,2);
2634 adjat(2,0) =
a(0,1)*
a(1,2)-
a(0,2)*
a(1,1);
2636 adjat(0,1) =
a(1,2)*
a(2,0)-
a(1,0)*
a(2,2);
2637 adjat(1,1) =
a(0,0)*
a(2,2)-
a(0,2)*
a(2,0);
2638 adjat(2,1) =
a(0,2)*
a(1,0)-
a(0,0)*
a(1,2);
2640 adjat(0,2) =
a(1,0)*
a(2,1)-
a(1,1)*
a(2,0);
2641 adjat(1,2) =
a(0,1)*
a(2,0)-
a(0,0)*
a(2,1);
2642 adjat(2,2) =
a(0,0)*
a(1,1)-
a(0,1)*
a(1,0);
2648 MFEM_ASSERT(
a.Width() <=
a.Height() &&
a.Width() >= 1 &&
a.Height() <= 3,
"");
2649 MFEM_ASSERT(inva.
Height() ==
a.Width(),
"incorrect dimensions");
2650 MFEM_ASSERT(inva.
Width() ==
a.Height(),
"incorrect dimensions");
2652 if (
a.Width() <
a.Height())
2656 if (
a.Height() == 2)
2676 MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(
a.FNorm()/
a.Width(),
a.Width()),
2677 "singular matrix!");
2683 inva(0,0) = 1.0 /
a.
Det();
2697 if ( (
a.Width() !=
a.Height()) || ( (
a.Height()!= 1) && (
a.Height()!= 2)
2698 && (
a.Height()!= 3) ) )
2700 mfem_error(
"CalcInverseTranspose(...): dimension mismatch");
2709 inva(0,0) = 1.0 /
a(0,0);
2712 inva(0,0) =
a(1,1) * t ;
2713 inva(1,0) = -
a(0,1) * t ;
2714 inva(0,1) = -
a(1,0) * t ;
2715 inva(1,1) =
a(0,0) * t ;
2718 inva(0,0) = (
a(1,1)*
a(2,2)-
a(1,2)*
a(2,1))*t;
2719 inva(1,0) = (
a(0,2)*
a(2,1)-
a(0,1)*
a(2,2))*t;
2720 inva(2,0) = (
a(0,1)*
a(1,2)-
a(0,2)*
a(1,1))*t;
2722 inva(0,1) = (
a(1,2)*
a(2,0)-
a(1,0)*
a(2,2))*t;
2723 inva(1,1) = (
a(0,0)*
a(2,2)-
a(0,2)*
a(2,0))*t;
2724 inva(2,1) = (
a(0,2)*
a(1,0)-
a(0,0)*
a(1,2))*t;
2726 inva(0,2) = (
a(1,0)*
a(2,1)-
a(1,1)*
a(2,0))*t;
2727 inva(1,2) = (
a(0,1)*
a(2,0)-
a(0,0)*
a(2,1))*t;
2728 inva(2,2) = (
a(0,0)*
a(1,1)-
a(0,1)*
a(1,0))*t;
2738 "Matrix must be 3x2 or 2x1, "
2739 <<
"and the Vector must be sized with the rows. "
2740 <<
" J.Height() = " << J.
Height()
2741 <<
", J.Width() = " << J.
Width()
2742 <<
", n.Size() = " << n.
Size()
2753 n(0) = d[1]*d[5] - d[2]*d[4];
2754 n(1) = d[2]*d[3] - d[0]*d[5];
2755 n(2) = d[0]*d[4] - d[1]*d[3];
2761 const int height =
a.Height();
2762 const int width =
a.Width();
2763 for (
int i = 0; i < height; i++)
2765 for (
int j = 0; j <= i; j++)
2768 for (
int k = 0; k < width; k++)
2770 temp +=
a(i,k) *
a(j,k);
2772 aat(j,i) = aat(i,j) = temp;
2779 for (
int i = 0; i < A.
Height(); i++)
2781 for (
int j = 0; j < i; j++)
2784 for (
int k = 0; k < A.
Width(); k++)
2786 t += D(k) * A(i, k) * A(j, k);
2794 for (
int i = 0; i < A.
Height(); i++)
2797 for (
int k = 0; k < A.
Width(); k++)
2799 t += D(k) * A(i, k) * A(i, k);
2807 for (
int i = 0; i < A.
Height(); i++)
2809 for (
int j = 0; j <= i; j++)
2812 for (
int k = 0; k < A.
Width(); k++)
2814 t += D(k) * A(i, k) * A(j, k);
2816 ADAt(j, i) = ADAt(i, j) = t;
2827 mfem_error(
"MultABt(...): dimension mismatch");
2831#ifdef MFEM_USE_LAPACK
2832 static char transa =
'N', transb =
'T';
2836 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &m,
2839 const int ah = A.
Height();
2840 const int bh = B.
Height();
2841 const int aw = A.
Width();
2848 const int ah = A.
Height();
2849 const int bh = B.
Height();
2850 const int aw = A.
Width();
2855 for (
int j = 0; j < bh; j++)
2856 for (
int i = 0; i < ah; i++)
2859 const real_t *ap = ad + i;
2860 const real_t *bp = bd + j;
2861 for (
int k = 0; k < aw; k++)
2873 for (i = 0; i < A.
Height(); i++)
2874 for (j = 0; j < B.
Height(); j++)
2877 for (k = 0; k < A.
Width(); k++)
2879 d += A(i, k) * B(j, k);
2893 mfem_error(
"MultADBt(...): dimension mismatch");
2897 const int ah = A.
Height();
2898 const int bh = B.
Height();
2899 const int aw = A.
Width();
2905 for (
int i = 0, s = ah*bh; i < s; i++)
2909 for (
int k = 0; k < aw; k++)
2912 for (
int j = 0; j < bh; j++)
2914 const real_t dk_bjk = dd[k] * bd[j];
2915 for (
int i = 0; i < ah; i++)
2917 cp[i] += ad[i] * dk_bjk;
2932 mfem_error(
"AddMultABt(...): dimension mismatch");
2936#ifdef MFEM_USE_LAPACK
2937 static char transa =
'N', transb =
'T';
2941 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &m,
2944 const int ah = A.
Height();
2945 const int bh = B.
Height();
2946 const int aw = A.
Width();
2951 for (
int k = 0; k < aw; k++)
2954 for (
int j = 0; j < bh; j++)
2956 const real_t bjk = bd[j];
2957 for (
int i = 0; i < ah; i++)
2959 cp[i] += ad[i] * bjk;
2970 for (i = 0; i < A.
Height(); i++)
2971 for (j = 0; j < B.
Height(); j++)
2974 for (k = 0; k < A.
Width(); k++)
2976 d += A(i, k) * B(j, k);
2990 mfem_error(
"AddMultADBt(...): dimension mismatch");
2994 const int ah = A.
Height();
2995 const int bh = B.
Height();
2996 const int aw = A.
Width();
3002 for (
int k = 0; k < aw; k++)
3005 for (
int j = 0; j < bh; j++)
3007 const real_t dk_bjk = dd[k] * bd[j];
3008 for (
int i = 0; i < ah; i++)
3010 cp[i] += ad[i] * dk_bjk;
3026 mfem_error(
"AddMult_a_ABt(...): dimension mismatch");
3030#ifdef MFEM_USE_LAPACK
3031 static char transa =
'N', transb =
'T';
3036 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &m,
3039 const int ah = A.
Height();
3040 const int bh = B.
Height();
3041 const int aw = A.
Width();
3046 for (
int k = 0; k < aw; k++)
3049 for (
int j = 0; j < bh; j++)
3052 for (
int i = 0; i < ah; i++)
3054 cp[i] += ad[i] * bjk;
3065 for (i = 0; i < A.
Height(); i++)
3066 for (j = 0; j < B.
Height(); j++)
3069 for (k = 0; k < A.
Width(); k++)
3071 d += A(i, k) * B(j, k);
3084 mfem_error(
"MultAtB(...): dimension mismatch");
3088#ifdef MFEM_USE_LAPACK
3089 static char transa =
'T', transb =
'N';
3093 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &k,
3096 const int ah = A.
Height();
3097 const int aw = A.
Width();
3098 const int bw = B.
Width();
3103 for (
int j = 0; j < bw; j++)
3106 for (
int i = 0; i < aw; i++)
3109 for (
int k = 0; k < ah; k++)
3122 for (i = 0; i < A.
Width(); i++)
3123 for (j = 0; j < B.
Width(); j++)
3126 for (k = 0; k < A.
Height(); k++)
3128 d += A(k, i) * B(k, j);
3141#ifdef MFEM_USE_LAPACK
3142 static char transa =
'T', transb =
'N';
3146 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &k,
3149 const int ah = A.
Height();
3150 const int aw = A.
Width();
3151 const int bw = B.
Width();
3156 for (
int j = 0; j < bw; j++)
3159 for (
int i = 0; i < aw; i++)
3162 for (
int k = 0; k < ah; k++)
3180#ifdef MFEM_USE_LAPACK
3181 static char transa =
'T', transb =
'N';
3186 MFEM_LAPACK_PREFIX(gemm_)(&transa, &transb, &m, &n, &k, &
alpha, A.
Data(), &k,
3189 const int ah = A.
Height();
3190 const int aw = A.
Width();
3191 const int bw = B.
Width();
3196 for (
int j = 0; j < bw; j++)
3199 for (
int i = 0; i < aw; i++)
3202 for (
int k = 0; k < ah; k++)
3218 for (
int i = 0; i < A.
Height(); i++)
3220 for (
int j = 0; j < i; j++)
3223 for (
int k = 0; k < A.
Width(); k++)
3225 d += A(i,k) * A(j,k);
3227 AAt(i, j) += (d *=
a);
3231 for (
int k = 0; k < A.
Width(); k++)
3233 d += A(i,k) * A(i,k);
3241 for (
int i = 0; i < A.
Height(); i++)
3243 for (
int j = 0; j <= i; j++)
3246 for (
int k = 0; k < A.
Width(); k++)
3248 d += A(i,k) * A(j,k);
3250 AAt(i, j) = AAt(j, i) =
a * d;
3257 for (
int i = 0; i < v.
Size(); i++)
3259 for (
int j = 0; j <= i; j++)
3261 vvt(i,j) = vvt(j,i) = v(i) * v(j);
3271 mfem_error(
"MultVWt(...): dimension mismatch");
3275 for (
int i = 0; i < v.
Size(); i++)
3278 for (
int j = 0; j < w.
Size(); j++)
3280 VWt(i, j) = vi * w(j);
3287 const int m = v.
Size(), n = w.
Size();
3292 mfem_error(
"AddMultVWt(...): dimension mismatch");
3296 for (
int i = 0; i < m; i++)
3299 for (
int j = 0; j < n; j++)
3301 VWt(i, j) += vi * w(j);
3308 const int n = v.
Size();
3313 mfem_error(
"AddMultVVt(...): dimension mismatch");
3317 for (
int i = 0; i < n; i++)
3320 for (
int j = 0; j < i; j++)
3322 const real_t vivj = vi * v(j);
3326 VVt(i, i) += vi * vi;
3333 const int m = v.
Size(), n = w.
Size();
3338 mfem_error(
"AddMult_a_VWt(...): dimension mismatch");
3342 for (
int j = 0; j < n; j++)
3345 for (
int i = 0; i < m; i++)
3347 VWt(i, j) += v(i) * awj;
3355 "incompatible dimensions!");
3357 const int n = v.
Size();
3358 for (
int i = 0; i < n; i++)
3361 for (
int j = 0; j < i; j++)
3363 const real_t avivj = avi * v(j);
3367 VVt(i, i) += avi * v(i);
3390#ifdef MFEM_USE_LAPACK
3392 if (m) { MFEM_LAPACK_PREFIX(getrf_)(&m, &m,
data, &m,
ipiv, &info); }
3397 for (
int i = 0; i < m; i++)
3402 real_t a = std::abs(data_ptr[piv+i*m]);
3403 for (
int j = i+1; j < m; j++)
3405 const real_t b = std::abs(data_ptr[j+i*m]);
3416 for (
int j = 0; j < m; j++)
3423 if (abs(data_ptr[i + i*m]) <= TOL)
3428 const real_t a_ii_inv = 1.0 / data_ptr[i+i*m];
3429 for (
int j = i+1; j < m; j++)
3431 data_ptr[j+i*m] *= a_ii_inv;
3433 for (
int k = i+1; k < m; k++)
3435 const real_t a_ik = data_ptr[i+k*m];
3436 for (
int j = i+1; j < m; j++)
3438 data_ptr[j+k*m] -= a_ik * data_ptr[j+i*m];
3450 for (
int i=0; i<m; i++)
3454 det *= -
data[m * i + i];
3458 det *=
data[m * i + i];
3467 for (
int k = 0; k < n; k++)
3470 for (
int i = 0; i < m; i++)
3473 for (
int j = i+1; j < m; j++)
3475 x_i += x[j] *
data[i+j*m];
3480 for (
int i = m-1; i >= 0; i--)
3483 for (
int j = 0; j < i; j++)
3485 x_i += x[j] *
data[i+j*m];
3490 for (
int i = m-1; i >= 0; i--)
3501 for (
int k = 0; k < n; k++)
3511 for (
int k = 0; k < n; k++)
3520#ifdef MFEM_USE_LAPACK
3525 MFEM_LAPACK_PREFIX(getrs_)(&
trans, &m, &n,
data, &m,
ipiv, X, &m, &info);
3527 MFEM_VERIFY(!info,
"LAPACK: error in DGETRS");
3538#ifdef MFEM_USE_LAPACK
3539 char n_ch =
'N', side =
'R', u_ch =
'U', l_ch =
'L';
3543 MFEM_LAPACK_PREFIX(trsm_)(&side,&u_ch,&n_ch,&n_ch,&n,&m,&
alpha,
data,&m,X,&n);
3544 MFEM_LAPACK_PREFIX(trsm_)(&side,&l_ch,&n_ch,&u_ch,&n,&m,&
alpha,
data,&m,X,&n);
3550 for (
int k = 0; k < n; k++)
3552 for (
int j = 0; j < m; j++)
3555 for (
int i = j+1; i < m; i++)
3557 x[i*n] -=
data[j + i*m] * x_j;
3565 for (
int k = 0; k < n; k++)
3567 for (
int j = m-1; j >= 0; j--)
3569 const real_t x_j = x[j*n];
3570 for (
int i = 0; i < j; i++)
3572 x[i*n] -=
data[j + i*m] * x_j;
3580 for (
int k = 0; k < n; k++)
3582 for (
int i = m-1; i >= 0; --i)
3595 for (
int k = 0; k < m; k++)
3597 const real_t minus_x_k = -( x[k] = 1.0/
data[k+k*m] );
3598 for (
int i = 0; i < k; i++)
3600 x[i] =
data[i+k*m] * minus_x_k;
3602 for (
int j = k-1; j >= 0; j--)
3605 for (
int i = 0; i < j; i++)
3607 x[i] -=
data[i+j*m] * x_j;
3615 for (
int j = 0; j < k; j++)
3618 for (
int i = 0; i <= j; i++)
3620 X[i+j*m] += X[i+k*m] * minus_L_kj;
3622 for (
int i = j+1; i < m; i++)
3624 X[i+j*m] = X[i+k*m] * minus_L_kj;
3628 for (
int k = m-2; k >= 0; k--)
3630 for (
int j = 0; j < k; j++)
3633 for (
int i = 0; i < m; i++)
3635 X[i+j*m] -= X[i+k*m] * L_kj;
3640 for (
int k = m-1; k >= 0; k--)
3645 for (
int i = 0; i < m; i++)
3671 SubMult(m, n, r, L21, B1, B2);
3678 SubMult(n, m, r, U12, X2, Y1);
3686#ifdef MFEM_USE_LAPACK
3689 MFEM_VERIFY(
data,
"Matrix data not set");
3690 if (m) { MFEM_LAPACK_PREFIX(potrf_)(&uplo, &m,
data, &m, &info); }
3694 for (
int j = 0; j<m; j++)
3697 for (
int k = 0; k<j; k++)
3702 MFEM_VERIFY(
data[j+j*m] -
a > 0.,
3703 "CholeskyFactors::Factor: The matrix is not SPD");
3705 data[j+j*m] = std::sqrt(
data[j+j*m] -
a);
3707 if (
data[j + j*m] <= TOL)
3712 for (
int i = j+1; i<m; i++)
3715 for (
int k = 0; k<j; k++)
3729 for (
int i=0; i<m; i++)
3731 det *=
data[i + i*m];
3740 for (
int k = 0; k < n; k++)
3742 for (
int j = m-1; j >= 0; j--)
3745 for (
int i = 0; i < j; i++)
3747 x_j += x[i] *
data[j+i*m];
3758 for (
int k = 0; k < n; k++)
3760 for (
int i = 0; i < m; i++)
3763 for (
int j = i+1; j < m; j++)
3765 x_i += x[j] *
data[j+i*m];
3776#ifdef MFEM_USE_LAPACK
3782 MFEM_LAPACK_PREFIX(trtrs_)(&uplo, &
trans, &diag, &m, &n,
data, &m, X, &m,
3784 MFEM_VERIFY(!info,
"CholeskyFactors:LSolve:: info");
3788 for (
int k = 0; k < n; k++)
3791 for (
int j = 0; j < m; j++)
3794 for (
int i = j+1; i < m; i++)
3796 x[i] -=
data[i+j*m] * x_j;
3806#ifdef MFEM_USE_LAPACK
3813 MFEM_LAPACK_PREFIX(trtrs_)(&uplo, &
trans, &diag, &m, &n,
data, &m, X, &m,
3815 MFEM_VERIFY(!info,
"CholeskyFactors:USolve:: info");
3820 for (
int k = 0; k < n; k++)
3822 for (
int j = m-1; j >= 0; j--)
3825 for (
int i = 0; i < j; i++)
3827 x[i] -=
data[j+i*m] * x_j;
3837#ifdef MFEM_USE_LAPACK
3840 MFEM_LAPACK_PREFIX(potrs_)(&uplo, &m, &n,
data, &m, X, &m, &info);
3841 MFEM_VERIFY(!info,
"CholeskyFactors:Solve:: info");
3851#ifdef MFEM_USE_LAPACK
3861 MFEM_LAPACK_PREFIX(trsm_)(&side,&uplo,&transt,&diag,&n,&m,&
alpha,
data,&m,X,&n);
3862 MFEM_LAPACK_PREFIX(trsm_)(&side,&uplo,&
trans,&diag,&n,&m,&
alpha,
data,&m,X,&n);
3867 for (
int k = 0; k < n; k++)
3869 for (
int j = 0; j < m; j++)
3872 for (
int i = j+1; i < m; i++)
3874 x[i*n] -=
data[i + j*m] * x_j;
3881 for (
int k = 0; k < n; k++)
3883 for (
int j = m-1; j >= 0; j--)
3886 for (
int i = 0; i < j; i++)
3888 x[i*n] -=
data[j + i*m] * x_j;
3899#ifdef MFEM_USE_LAPACK
3901 for (
int i = 0; i<m; i++)
3903 for (
int j = i; j<m; j++)
3905 X[j+i*m] =
data[j+i*m];
3910 MFEM_LAPACK_PREFIX(potri_)(&uplo, &m, X, &m, &info);
3911 MFEM_VERIFY(!info,
"CholeskyFactors:GetInverseMatrix:: info");
3913 for (
int i = 0; i<m; i++)
3915 for (
int j = i+1; j<m; j++)
3917 X[i+j*m] = X[j+i*m];
3922 for (
int k = 0; k<m; k++)
3924 X[k+k*m] = 1./
data[k+k*m];
3925 for (
int i = k+1; i < m; i++)
3928 for (
int j=k; j<i; j++)
3930 s -=
data[i+j*m] * X[j+k*m]/
data[i+i*m];
3935 for (
int i = 0; i < m; i++)
3937 for (
int j = i; j < m; j++)
3940 for (
int k=j; k<m; k++)
3942 s += X[k+i*m] * X[k+j*m];
3944 X[i+j*m] = X[j+i*m] = s;
3951void DenseMatrixInverse::Init(
int m)
3959 factors =
new LUFactors();
3966 dynamic_cast<LUFactors *
>(factors)->ipiv =
new int[m];
3975 MFEM_ASSERT(
height ==
width,
"not a square matrix");
3984 MFEM_ASSERT(
height ==
width,
"not a square matrix");
3991 MFEM_ASSERT(a,
"DenseMatrix is not given");
3992 const real_t *adata = a->data;
3994 for (
int i = 0; i < s; i++)
3996 factors->
data[i] = adata[i];
4009 MFEM_VERIFY(mat.
height == mat.
width,
"DenseMatrix is not square!");
4013 if (own_data) {
delete [] factors->
data; }
4019 if (own_data) {
delete [] lu->
ipiv; }
4031 MFEM_VERIFY(
p != NULL,
"Operator is not a DenseMatrix!");
4037 for (
int row = 0; row <
height; row++)
4060 for (
int i = 0; i <
width; i++)
4071 delete [] factors->
data;
4074 delete []
dynamic_cast<LUFactors *
>(factors)->ipiv;
4080#ifdef MFEM_USE_LAPACK
4094 MFEM_LAPACK_PREFIX(syev_)(&jobz, &uplo, &n, EVect.
Data(), &n, EVal.
GetData(),
4095 &qwork, &lwork, &info);
4097 lwork = (int) qwork;
4098 work =
new real_t[lwork];
4103 : mat(other.mat), EVal(other.EVal), EVect(other.EVect), ev(NULL, other.n),
4108 lwork = other.lwork;
4110 work =
new real_t[lwork];
4116 if (mat.
Width() != n)
4118 mfem_error(
"DenseMatrixEigensystem::Eval(): dimension mismatch");
4123 MFEM_LAPACK_PREFIX(syev_)(&jobz, &uplo, &n, EVect.
Data(), &n, EVal.
GetData(),
4124 work, &lwork, &info);
4128 mfem::err <<
"DenseMatrixEigensystem::Eval(): DSYEV error code: "
4142 bool left_eigen_vectors,
4143 bool right_eigen_vectors)
4146 MFEM_VERIFY(A.
Height() == A.
Width(),
"A has to be a square matrix");
4147 MFEM_VERIFY(B.
Height() == B.
Width(),
"B has to be a square matrix");
4149 MFEM_VERIFY(B.
Height() == n,
"A and B dimension mismatch");
4155 if (left_eigen_vectors)
4160 if (right_eigen_vectors)
4173 int nl = max(1,Vl.
Height());
4174 int nr = max(1,Vr.
Height());
4176 MFEM_LAPACK_PREFIX(ggev_)(&jobvl,&jobvr,&n,A_copy.
Data(),&n,B_copy.
Data(),&n,
4177 alphar, alphai, beta, Vl.
Data(), &nl, Vr.
Data(),
4178 &nr, &qwork, &lwork, &info);
4180 lwork = (int) qwork;
4181 work =
new real_t[lwork];
4186 int nl = max(1,Vl.
Height());
4187 int nr = max(1,Vr.
Height());
4191 MFEM_LAPACK_PREFIX(ggev_)(&jobvl,&jobvr,&n,A_copy.
Data(),&n,B_copy.
Data(),&n,
4192 alphar, alphai, beta, Vl.
Data(), &nl, Vr.
Data(),
4193 &nr, work, &lwork, &info);
4196 mfem::err <<
"DenseMatrixGeneralizedEigensystem::Eval(): DGGEV error code: "
4202 for (
int i = 0; i<n; i++)
4206 evalues_r(i) = alphar[i]/
beta[i];
4207 evalues_i(i) = alphai[i]/
beta[i];
4226 bool left_singular_vectors,
4227 bool right_singular_vectors)
4231 jobu = (left_singular_vectors)?
'S' :
'N';
4232 jobvt = (right_singular_vectors)?
'S' :
'N';
4237 bool left_singular_vectors,
4238 bool right_singular_vectors)
4242 jobu = (left_singular_vectors)?
'S' :
'N';
4243 jobvt = (right_singular_vectors)?
'S' :
'N';
4248 char left_singular_vectors,
4249 char right_singular_vectors)
4253 jobu = left_singular_vectors;
4254 jobvt = right_singular_vectors;
4259 char left_singular_vectors,
4260 char right_singular_vectors)
4264 jobu = left_singular_vectors;
4265 jobvt = right_singular_vectors;
4269void DenseMatrixSVD::Init()
4274 MFEM_LAPACK_PREFIX(gesvd_)(&jobu, &jobvt, &m, &n, NULL, &m, sv.
GetData(),
4275 NULL, &m, NULL, &n, &qwork, &lwork, &info);
4276 lwork = (int) qwork;
4277 work =
new real_t[lwork];
4288 real_t * datau =
nullptr;
4289 real_t * datavt =
nullptr;
4295 else if (jobu ==
'S')
4305 else if (jobvt ==
'S')
4311 MFEM_LAPACK_PREFIX(gesvd_)(&jobu, &jobvt, &m, &n, Mc.
Data(), &m, sv.
GetData(),
4312 datau, &m, datavt, &n, work, &lwork, &info);
4316 mfem::err <<
"DenseMatrixSVD::Eval() : info = " << info << endl;
4333 const int *I = elem_dof.
GetI(), *J = elem_dof.
GetJ(), *dofs;
4341 for (
int i = 0; i < ne; i++)
4344 for (
int col = 0; col < n; col++)
4346 x_col = xp[dofs[col]];
4347 for (
int row = 0; row < n; row++)
4349 yp[dofs[row]] += x_col*d_col[row];
4358 for (
int i = 0; i < ne; i++)
4361 x_col = xp[dofs[0]];
4362 for (
int row = 0; row < n; row++)
4364 ye(row) = x_col*d_col[row];
4367 for (
int col = 1; col < n; col++)
4369 x_col = xp[dofs[col]];
4370 for (
int row = 0; row < n; row++)
4372 ye(row) += x_col*d_col[row];
4376 for (
int row = 0; row < n; row++)
4378 yp[dofs[row]] += ye(row);
4387 for (
int i=0; i<s; i++)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
static void LUFactor(DenseTensor &A, Array< int > &P)
Replaces the block diagonal matrix with its LU factors. The pivots are stored in P.
static void LUSolve(const DenseTensor &A, const Array< int > &P, Vector &x)
Replaces with , given the LU factors A and pivots P of the block-diagonal matrix .
void UMult(int m, int n, real_t *X) const
void Solve(int m, int n, real_t *X) const override
void RightSolve(int m, int n, real_t *X) const
void USolve(int m, int n, real_t *X) const
void LSolve(int m, int n, real_t *X) const
bool Factor(int m, real_t TOL=0.0) override
Compute the Cholesky factorization of the current matrix.
void GetInverseMatrix(int m, real_t *X) const override
Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1}.
void LMult(int m, int n, real_t *X) const
real_t Det(int m) const override
~DenseMatrixEigensystem()
DenseMatrixEigensystem(DenseMatrix &m)
~DenseMatrixGeneralizedEigensystem()
DenseMatrixGeneralizedEigensystem(DenseMatrix &a, DenseMatrix &b, bool left_eigen_vectors=false, bool right_eigen_vectors=false)
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
DenseMatrixInverse(bool spd_=false)
Default constructor.
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
real_t Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
MFEM_DEPRECATED DenseMatrixSVD(DenseMatrix &M, bool left_singular_vectors=false, bool right_singular_vectors=false)
Constructor for the DenseMatrixSVD.
void Eval(DenseMatrix &M)
Evaluate the SVD.
Data type dense matrix using column-major storage.
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
void CopyMNDiag(real_t c, int n, int row_offset, int col_offset)
Copy c on the diagonal of size n to *this at row_offset, col_offset.
void AddMult_a(real_t a, const Vector &x, Vector &y) const
y += a * A.x
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void AddMultTranspose_a(real_t a, const Vector &x, Vector &y) const
y += a * A^t x
void Set(real_t alpha, const real_t *A)
Set the matrix to alpha * A, assuming that A has the same dimensions as the matrix and uses column-ma...
void TestInversion()
Invert and print the numerical conditioning of the inversion.
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void Transpose()
(*this) = (*this)^t
void AddToVector(int offset, Vector &v) const
Add the matrix 'data' to the Vector 'v' at the given 'offset'.
void Threshold(real_t eps)
Replace small entries, abs(a_ij) <= eps, with zero.
const real_t * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
void CalcEigenvalues(real_t *lambda, real_t *vec) const
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
void SetRow(int r, const real_t *row)
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
real_t & operator()(int i, int j)
Returns reference to a_{ij}.
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
real_t * Data() const
Returns the matrix data array.
void GetSubMatrix(const Array< int > &idx, DenseMatrix &A) const
real_t * GetData() const
Returns the matrix data array.
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
y += a * A.x
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual ~DenseMatrix()
Destroys dense matrix.
friend class DenseMatrixInverse
void SetCol(int c, const real_t *col)
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
void Invert()
Replaces the current matrix with its inverse.
DenseMatrix & operator+=(const real_t *m)
void Neg()
(*this) = -(*this)
real_t operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
Copy matrix A^t to the location in *this at row_offset, col_offset.
void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
y += a * A^t x
void AdjustDofDirection(Array< int > &dofs)
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
void SingularValues(Vector &sv) const
real_t FNorm() const
Compute the Frobenius norm of the matrix.
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
void SetSubMatrix(const Array< int > &idx, const DenseMatrix &A)
Set (*this)(idx[i],idx[j]) = A(i,j)
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
void AddSubMatrix(const Array< int > &idx, const DenseMatrix &A)
(*this)(idx[i],idx[j]) += A(i,j)
real_t Trace() const
Trace of a square matrix.
void Diag(real_t c, int n)
Creates n x n diagonal matrix with diagonal elements c.
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
virtual void PrintMathematica(std::ostream &out=mfem::out) const
DenseMatrix & operator*=(real_t c)
void Swap(DenseMatrix &other)
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0<=i.
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
real_t & Elem(int i, int j) override
Returns reference to a_{ij}.
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
void InvSymmetricScaling(const Vector &s)
InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
int Rank(real_t tol) const
void Add(const real_t c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
void PrintMatlab(std::ostream &out=mfem::out) const override
Prints operator in Matlab format.
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
DenseMatrix & operator-=(const DenseMatrix &m)
real_t CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
void GradToCurl(DenseMatrix &curl)
void Print(std::ostream &out=mfem::out, int width_=4) const override
Prints matrix to stream out.
void GetColumn(int c, Vector &col) const
void GradToVectorCurl2D(DenseMatrix &curl)
MatrixInverse * Inverse() const override
Returns a pointer to the inverse matrix.
DenseMatrix & operator=(real_t c)
Sets the matrix elements equal to constant c.
real_t MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
void Norm2(real_t *v) const
Take the 2-norm of the columns of A and store in v.
void GetFromVector(int offset, const Vector &v)
Get the matrix 'data' from the Vector 'v' at the given 'offset'.
void GetRow(int r, Vector &row) const
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
void GradToDiv(Vector &div)
Rank 3 tensor (array of matrices)
void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const
DenseTensor & operator=(real_t c)
Sets the tensor elements equal to constant c.
void Swap(DenseTensor &t)
virtual void GetInverseMatrix(int m, real_t *X) const
virtual void Solve(int m, int n, real_t *X) const
virtual bool Factor(int m, real_t TOL=0.0)
A class to initialize the size of a Tensor.
void LSolve(int m, int n, real_t *X) const
static void SubMult(int m, int n, int r, const real_t *A21, const real_t *X1, real_t *X2)
bool Factor(int m, real_t TOL=0.0) override
Compute the LU factorization of the current matrix.
void Mult(int m, int n, real_t *X) const
void USolve(int m, int n, real_t *X) const
void BlockFactor(int m, int n, real_t *A12, real_t *A21, real_t *A22) const
real_t Det(int m) const override
void BlockForwSolve(int m, int n, int r, const real_t *L21, real_t *B1, real_t *B2) const
void Solve(int m, int n, real_t *X) const override
void RightSolve(int m, int n, real_t *X) const
void BlockBackSolve(int m, int n, int r, const real_t *U12, const real_t *X2, real_t *Y1) const
static const int ipiv_base
void GetInverseMatrix(int m, real_t *X) const override
Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
Abstract data type for matrix inverse.
Abstract data type matrix.
bool IsSquare() const
Returns whether the matrix is a square matrix.
int Capacity() const
Return the size of the allocated memory.
void Delete()
Delete the owned pointers and reset the Memory object.
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
int width
Dimension of the input / number of columns in the matrix.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int height
Dimension of the output / number of rows in the matrix.
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
void trans(const Vector &u, Vector &x)
MFEM_HOST_DEVICE void CalcInverse(const T *data, T *inv_data)
Return the inverse of a matrix with given size and data into the matrix with data inv_data.
MFEM_HOST_DEVICE real_t CalcSingularvalue< 2 >(const real_t *data, const int i)
Return the i'th singular value of the matrix of size 2 with given data.
MFEM_HOST_DEVICE void CalcEigenvalues< 2 >(const real_t *data, real_t *lambda, real_t *vec)
MFEM_HOST_DEVICE void Add(const int height, const int width, const TALPHA alpha, const TA *Adata, const TB *Bdata, TC *Cdata)
Compute C = A + alpha*B, where the matrices A, B and C are of size height x width with data Adata,...
MFEM_HOST_DEVICE void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix vector multiplication: y = A x, where the matrix A is of size height x width with given data,...
MFEM_HOST_DEVICE void CalcLeftInverse< 2, 1 >(const real_t *d, real_t *left_inv)
MFEM_HOST_DEVICE void LSolve(const real_t *data, const int m, const int *ipiv, real_t *x)
Assuming L.U = P.A factored matrix of size (m x m), compute X <- L^{-1} P X, for a vector X of length...
MFEM_HOST_DEVICE void MultABt(const int Aheight, const int Awidth, const int Bheight, const TA *Adata, const TB *Bdata, TC *ABtdata)
Multiply a matrix of size Aheight x Awidth and data Adata with the transpose of a matrix of size Bhei...
MFEM_HOST_DEVICE void USolve(const real_t *data, const int m, real_t *x)
Assuming L.U = P.A factored matrix of size (m x m), compute X <- U^{-1} X, for a vector X of length m...
MFEM_HOST_DEVICE void Symmetrize(const int size, T *data)
Symmetrize a square matrix with given size and data: A -> (A+A^T)/2.
MFEM_HOST_DEVICE void CalcLeftInverse< 3, 2 >(const real_t *d, real_t *left_inv)
MFEM_HOST_DEVICE void BlockFactor(const real_t *data, int m, const int *ipiv, int n, real_t *A12, real_t *A21, real_t *A22)
MFEM_HOST_DEVICE void CalcLeftInverse< 3, 1 >(const real_t *d, real_t *left_inv)
MFEM_HOST_DEVICE real_t CalcSingularvalue< 3 >(const real_t *data, const int i)
Return the i'th singular value of the matrix of size 3 with given data.
MFEM_HOST_DEVICE void CalcEigenvalues< 3 >(const real_t *data, real_t *lambda, real_t *vec)
MFEM_HOST_DEVICE void SubMult(const int m, const int n, const int r, const real_t *A21, const real_t *X1, real_t *X2)
Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices X1, and X2 of size (m x r) and (...
void CalcOrtho(const DenseMatrix &J, Vector &n)
void Swap(Array< T > &, Array< T > &)
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
void AddMult_a_ABt(real_t a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
real_t u(const Vector &xvec)
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void BatchLUSolve(const DenseTensor &Mlu, const Array< int > &P, Vector &X)
Solve batch linear systems. Calls BatchedLinAlg::LUSolve.
void AddMult_a(real_t alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
const T * HostRead(const Memory< T > &mem, int size)
Shortcut to Read(const Memory<T> &mem, int size, false)
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void dsygv_Eigensystem(DenseMatrix &a, DenseMatrix &b, Vector &ev, DenseMatrix *evect)
void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
void AddMult_a_VWt(const real_t a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
void AddMult_a_VVt(const real_t a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
void BatchLUFactor(DenseTensor &Mlu, Array< int > &P, const real_t TOL)
Compute the LU factorization of a batch of matrices. Calls BatchedLinAlg::LUFactor.
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
int CheckFinite(const real_t *v, const int n)
void AddMult_a_AtB(real_t a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
AtB += a * A^t * B.
void dsyev_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
bool LinearSolve(DenseMatrix &A, real_t *X, real_t TOL)
Solves the dense linear system, A * X = B for X
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
void AddMultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
void AddMultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
AtB += A^t * B.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
void Mult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
real_t p(const Vector &x, real_t t)