Solver rewrite

include/cantera/numerics/Newton.h

@@ -31,24 +31,6 @@ class Newton
     //! at `x`, but the Jacobian is not recomputed.
     void step(doublereal* x, doublereal* step, int loglevel);
 
-    /**
-     * Return the factor by which the undamped Newton step 'step0'
-     * must be multiplied in order to keep all solution components in
-     * all domains between their specified lower and upper bounds.
-     */
-    doublereal boundStep(const doublereal* x0, const doublereal* step0, int loglevel);
-
-    /**
-     * On entry, step0 must contain an undamped Newton step for the solution x0.
-     * This method attempts to find a damping coefficient such that the next
-     * undamped step would have a norm smaller than that of step0. If
-     * successful, the new solution after taking the damped step is returned in
-     * x1, and the undamped step at x1 is returned in step1.
-     */
-    int dampStep(const doublereal* x0, const doublereal* step0,
-                 doublereal* x1, doublereal* step1, doublereal& s1,
-                 int loglevel, bool writetitle);
-
     //! Compute the weighted 2-norm of `step`.
     doublereal weightedNorm(const doublereal* x, const doublereal* step) const;
 
@@ -75,8 +57,8 @@ class Newton
 
     //! Set upper and lower bounds on the nth component
     void setBounds(size_t n, doublereal lower, doublereal upper) {
-        m_min[n] = lower;
-        m_max[n] = upper;
+        m_lower_bounds[n] = lower;
+        m_upper_bounds[n] = upper;
     }
 
     void setConstants(vector_int constantComponents) {
@@ -98,17 +80,17 @@ class Newton
     size_t m_nv;
 
     //! solution converged if [weightedNorm(sol, step) < m_convergenceThreshold]
-    doublereal m_convergenceThreshold;
+    doublereal m_converge_tol;
 
     DenseMatrix m_jacobian, m_jacFactored;
-    int m_jacAge, m_jacMaxAge;
+    size_t m_jacAge, m_jacMaxAge;
     doublereal m_jacRtol, m_jacAtol;
 
 
     //! work arrays of size #m_nv used in solve().
     vector_fp m_x, m_x1, m_stp, m_stp1;
 
-    vector_fp m_max, m_min;
+    vector_fp m_upper_bounds, m_lower_bounds;
     vector_fp m_rtol_ss, m_rtol_ts;
     vector_fp m_atol_ss, m_atol_ts;
 
@@ -120,6 +102,18 @@ class Newton
     //! current timestep reciprocal
     doublereal m_rdt = 0;
 };
+
+// //! Returns the weighted Root Mean Square Deviation given a vector of residuals and
+// //  vectors of the corresponding weights and absolute tolerances for each component.
+// double weightedRMS(vector_fp residuals, vector_fp weights, vector_fp atols) {
+//     size_t n = residuals.size();
+//     double square = 0.0;
+//     for (size_t i = 0; i < n; i++) {
+//         square += pow(residuals[i]/(weights[i] + atols[i]), 2);
+//     }
+//     return sqrt(square/n);
+// }
+
 }
 
 #endif

src/numerics/DenseMatrix.cpp

@@ -209,6 +209,17 @@ int solveFactored(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
             }
         }
     #else
+        MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns());
+        #ifdef NDEBUG
+            auto lu = Am.partialPivLu();
+        #else
+            auto lu = Am.fullPivLu();
+            if (lu.nonzeroPivots() < static_cast<long int>(A.nColumns())) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                    "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}",
+                    lu.nonzeroPivots(), A.nColumns());
+            }
+        #endif
         for (size_t i = 0; i < nrhs; i++) {
             MappedVector bm(b + ldb*i, A.nColumns());
             bm = lu.solve(bm);