Skip to content
Snippets Groups Projects
Commit f19172aa authored by Tobias Leibner's avatar Tobias Leibner
Browse files

[algorithms.cholesky] also throw on NaN input

parent 8c660bda
No related branches found
No related tags found
No related merge requests found
...@@ -18,7 +18,7 @@ status = 1a3bcab04b011a5d6e44f9983cae6ff89fa695e8 bin (heads/master) ...@@ -18,7 +18,7 @@ status = 1a3bcab04b011a5d6e44f9983cae6ff89fa695e8 bin (heads/master)
b8cf8efc32e46511990b691afee55da3b7af15c8 dune-xt-data (heads/master) b8cf8efc32e46511990b691afee55da3b7af15c8 dune-xt-data (heads/master)
0c9df39934e23b950f357c912b4c800b86de683f dune-xt-functions (heads/dailywork_tleibner) 0c9df39934e23b950f357c912b4c800b86de683f dune-xt-functions (heads/dailywork_tleibner)
dd30fcd7d4485eb2a8158d5ddf01333f58502c40 dune-xt-grid (heads/master) dd30fcd7d4485eb2a8158d5ddf01333f58502c40 dune-xt-grid (heads/master)
+9762d4f7bfcaa3424fc4aeaf362fa2b7710898ea dune-xt-la (heads/dailywork-tleibner) +8c660bdad165a189b44c3f6e92e73e07cff1d774 dune-xt-la (heads/dailywork-tleibner)
09d0378f616b94d68bcdd9fc6114813181849ec0 scripts (heads/master) 09d0378f616b94d68bcdd9fc6114813181849ec0 scripts (heads/master)
commit = 02dfba6e3d4fa4659a18ce0f63b5956e1e557f16 commit = 02dfba6e3d4fa4659a18ce0f63b5956e1e557f16
...@@ -115,7 +115,7 @@ commit = dd30fcd7d4485eb2a8158d5ddf01333f58502c40 ...@@ -115,7 +115,7 @@ commit = dd30fcd7d4485eb2a8158d5ddf01333f58502c40
[submodule.dune-xt-la] [submodule.dune-xt-la]
remote = git@github.com:dune-community/dune-xt-la.git remote = git@github.com:dune-community/dune-xt-la.git
status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master) status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master)
commit = 9762d4f7bfcaa3424fc4aeaf362fa2b7710898ea commit = 8c660bdad165a189b44c3f6e92e73e07cff1d774
[submodule.scripts] [submodule.scripts]
remote = https://github.com/wwu-numerik/scripts.git remote = https://github.com/wwu-numerik/scripts.git
......
...@@ -37,7 +37,7 @@ void tridiagonal_ldlt(FirstVectorType& diag, SecondVectorType& subdiag) ...@@ -37,7 +37,7 @@ void tridiagonal_ldlt(FirstVectorType& diag, SecondVectorType& subdiag)
typedef Common::VectorAbstraction<SecondVectorType> V2; typedef Common::VectorAbstraction<SecondVectorType> V2;
V1::set_entry(diag, 0, V1::get_entry(diag, 0)); V1::set_entry(diag, 0, V1::get_entry(diag, 0));
for (size_t ii = 0; ii < diag.size() - 1; ++ii) { for (size_t ii = 0; ii < diag.size() - 1; ++ii) {
if (V1::get_entry(diag, ii) <= 0) if (!(V1::get_entry(diag, ii) > 0)) // use !(.. > 0) instead of (.. <= 0) to also throw on NaNs
DUNE_THROW(MathError, "LDL^T factorization failed!"); DUNE_THROW(MathError, "LDL^T factorization failed!");
V2::set_entry(subdiag, ii, V2::get_entry(subdiag, ii) / V1::get_entry(diag, ii)); V2::set_entry(subdiag, ii, V2::get_entry(subdiag, ii) / V1::get_entry(diag, ii));
V1::set_entry( V1::set_entry(
...@@ -104,7 +104,7 @@ void cholesky_rowwise(MatrixType& A) ...@@ -104,7 +104,7 @@ void cholesky_rowwise(MatrixType& A)
auto L_ii = M::get_entry(A, ii, ii); auto L_ii = M::get_entry(A, ii, ii);
for (size_t kk = 0; kk < ii; ++kk) for (size_t kk = 0; kk < ii; ++kk)
L_ii -= std::pow(M::get_entry(L, ii, kk), 2); L_ii -= std::pow(M::get_entry(L, ii, kk), 2);
if (L_ii <= 0) if (!(L_ii > 0)) // use !(.. > 0) instead of (.. <= 0) to also throw on NaNs
DUNE_THROW(MathError, "Cholesky factorization failed!"); DUNE_THROW(MathError, "Cholesky factorization failed!");
M::set_entry(L, ii, ii, std::sqrt(L_ii)); M::set_entry(L, ii, ii, std::sqrt(L_ii));
} // ii } // ii
...@@ -120,7 +120,7 @@ void cholesky_colwise(MatrixType& A) ...@@ -120,7 +120,7 @@ void cholesky_colwise(MatrixType& A)
auto L_jj = M::get_entry(A, jj, jj); auto L_jj = M::get_entry(A, jj, jj);
for (size_t kk = 0; kk < jj; ++kk) for (size_t kk = 0; kk < jj; ++kk)
L_jj -= std::pow(M::get_entry(L, jj, kk), 2); L_jj -= std::pow(M::get_entry(L, jj, kk), 2);
if (L_jj <= 0) if (!(L_jj > 0)) // use !(.. > 0) instead of (.. <= 0) to also throw on NaNs
DUNE_THROW(MathError, "Cholesky factorization failed!"); DUNE_THROW(MathError, "Cholesky factorization failed!");
L_jj = std::sqrt(L_jj); L_jj = std::sqrt(L_jj);
M::set_entry(L, jj, jj, L_jj); M::set_entry(L, jj, jj, L_jj);
...@@ -167,7 +167,7 @@ cholesky_csr(MatrixType& A) ...@@ -167,7 +167,7 @@ cholesky_csr(MatrixType& A)
auto kk = row_pointers[ii]; auto kk = row_pointers[ii];
while (kk < row_pointers[ii + 1] && size_t(column_indices[kk]) < ii) while (kk < row_pointers[ii + 1] && size_t(column_indices[kk]) < ii)
L_ii -= std::pow(entries[kk++], 2); L_ii -= std::pow(entries[kk++], 2);
if (L_ii <= 0) if (!(L_ii > 0)) // use !(.. > 0) instead of (.. <= 0) to also throw on NaNs
DUNE_THROW(MathError, "Cholesky factorization failed!"); DUNE_THROW(MathError, "Cholesky factorization failed!");
M::set_entry(L, ii, ii, std::sqrt(L_ii)); M::set_entry(L, ii, ii, std::sqrt(L_ii));
} // ii } // ii
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment