From b9ac39bf68980acaa992056888548c575d58e476 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 11 Nov 2016 19:07:15 -0500 Subject: [PATCH] [Numerics] Failures in BandMatrix raise exceptions Store the 'info' flag returned by the underlying LAPACK function so that callers can catch the exception and retrieve this information if desired. Eliminate automatic 'bandmatrix.csv' output. --- include/cantera/numerics/BandMatrix.h | 5 +++ src/numerics/BandMatrix.cpp | 62 +++++++++++++-------------- src/oneD/MultiNewton.cpp | 40 +++++++++-------- 3 files changed, 55 insertions(+), 52 deletions(-) diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h index a7e8a4f5c..1c04dfd03 100644 --- a/include/cantera/numerics/BandMatrix.h +++ b/include/cantera/numerics/BandMatrix.h @@ -280,6 +280,9 @@ public: */ virtual size_t checkColumns(doublereal& valueSmall) const; + //! LAPACK "info" flag after last factor/solve operation + int info() const { return m_info; }; + protected: //! Matrix data vector_fp data; @@ -311,6 +314,8 @@ protected: //! Extra dp work array needed - size = 3n vector_fp work_; + + int m_info; }; //! Utility routine to print out the matrix diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index 4830d257b..c2ec04e7d 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -30,7 +30,8 @@ BandMatrix::BandMatrix() : m_n(0), m_kl(0), m_ku(0), - m_zero(0.0) + m_zero(0.0), + m_info(0) { } @@ -38,7 +39,8 @@ BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) : m_n(n), m_kl(kl), m_ku(ku), - m_zero(0.0) + m_zero(0.0), + m_info(0) { data.resize(n*(2*kl + ku + 1)); ludata.resize(n*(2*kl + ku + 1)); @@ -59,7 +61,8 @@ BandMatrix::BandMatrix(const BandMatrix& y) : m_n(0), m_kl(0), m_ku(0), - m_zero(0.0) + m_zero(0.0), + m_info(y.m_info) { m_n = y.m_n; m_kl = y.m_kl; @@ -95,6 +98,7 @@ BandMatrix& BandMatrix::operator=(const BandMatrix& y) m_colPtrs[j] = &data[ldab * j]; m_lu_col_ptrs[j] = &ludata[ldab * j]; } + m_info = y.m_info; return *this; } @@ -235,27 +239,23 @@ void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) con int BandMatrix::factor() { - int info=0; ludata = data; #if CT_USE_LAPACK ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(), - ludata.data(), ldim(), ipiv().data(), info); + ludata.data(), ldim(), ipiv().data(), m_info); #else long int nu = static_cast(nSuperDiagonals()); long int nl = static_cast(nSubDiagonals()); long int smu = nu + nl; - info = bandGBTRF(m_lu_col_ptrs.data(), static_cast(nColumns()), - nu, nl, smu, m_ipiv.data()); + m_info = bandGBTRF(m_lu_col_ptrs.data(), static_cast(nColumns()), + nu, nl, smu, m_ipiv.data()); #endif - // if info = 0, LU decomp succeeded. - if (info == 0) { - m_factored = true; - } else { - m_factored = false; - ofstream fout("bandmatrix.csv"); - fout << *this << endl; + if (m_info != 0) { + throw Cantera::CanteraError("BandMatrix::factor", + "Factorization failed with DGBTRF error code {}.", m_info); } - return info; + m_factored = true; + return m_info; } int BandMatrix::solve(const doublereal* const b, doublereal* const x) @@ -266,34 +266,30 @@ int BandMatrix::solve(const doublereal* const b, doublereal* const x) int BandMatrix::solve(doublereal* b, size_t nrhs, size_t ldb) { - int info = 0; if (!m_factored) { - info = factor(); + factor(); } if (ldb == 0) { ldb = nColumns(); } - if (info == 0) { #if CT_USE_LAPACK - ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(), - nSuperDiagonals(), nrhs, ludata.data(), ldim(), - ipiv().data(), b, ldb, info); + ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(), + nSuperDiagonals(), nrhs, ludata.data(), ldim(), + ipiv().data(), b, ldb, m_info); #else - long int nu = static_cast(nSuperDiagonals()); - long int nl = static_cast(nSubDiagonals()); - long int smu = nu + nl; - double** a = m_lu_col_ptrs.data(); - bandGBTRS(a, static_cast(nColumns()), smu, nl, m_ipiv.data(), - b); + long int nu = static_cast(nSuperDiagonals()); + long int nl = static_cast(nSubDiagonals()); + long int smu = nu + nl; + double** a = m_lu_col_ptrs.data(); + bandGBTRS(a, static_cast(nColumns()), smu, nl, m_ipiv.data(), b); + m_info = 0; #endif - } - // error handling - if (info != 0) { - ofstream fout("bandmatrix.csv"); - fout << *this << endl; + if (m_info != 0) { + throw Cantera::CanteraError("BandMatrix::solve", + "Linear solve failed with DGBTRS error code {}.", m_info); } - return info; + return m_info; } vector_fp::iterator BandMatrix::begin() diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index fbb07c0ab..ca97e8820 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -167,27 +167,29 @@ void MultiNewton::step(doublereal* x, doublereal* step, step[n] = -step[n]; } - size_t iok = jac.solve(step, step); - // if iok is non-zero, then solve failed - if (iok != 0) { - iok--; - size_t nd = r.nDomains(); - size_t n; - for (n = nd-1; n != npos; n--) { - if (iok >= r.start(n)) { - break; + try { + jac.solve(step, step); + } catch (CanteraError&) { + int iok = jac.info() - 1; + if (iok >= 0) { + size_t nd = r.nDomains(); + size_t n; + for (n = nd-1; n != npos; n--) { + if (iok >= static_cast(r.start(n))) { + break; + } } + Domain1D& dom = r.domain(n); + size_t offset = iok - r.start(n); + size_t pt = offset/dom.nComponents(); + size_t comp = offset - pt*dom.nComponents(); + throw CanteraError("MultiNewton::step", + "Jacobian is singular for domain {}, component {} at point {}\n" + "(Matrix row {})", + dom.id(), dom.componentName(comp), pt, iok); + } else { + throw; } - Domain1D& dom = r.domain(n); - size_t offset = iok - r.start(n); - size_t pt = offset/dom.nComponents(); - size_t comp = offset - pt*dom.nComponents(); - throw CanteraError("MultiNewton::step", - "Jacobian is singular for domain {}, component {} at point {}\n" - "(Matrix row {}) \nsee file bandmatrix.csv\n", - dom.id(), dom.componentName(comp), pt, iok); - } else if (int(iok) < 0) { - throw CanteraError("MultiNewton::step", "iok = {}", iok); } }