Initial timestep functionality and jac optimizations

include/cantera/numerics/DenseMatrix.h

@@ -155,6 +155,8 @@ class DenseMatrix : public Array2D
     friend int invert(DenseMatrix& A, int nn);
 };
 
+int factor(DenseMatrix& A);
+int solveFactored(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0);
 
 //! Solve Ax = b. Array b is overwritten on exit with x.
 /*!

include/cantera/numerics/Newton.h

@@ -52,6 +52,10 @@ class Newton
     //! Compute the weighted 2-norm of `step`.
     doublereal weightedNorm(const doublereal* x, const doublereal* step) const;
 
+    int hybridSolve();
+
+    int timestep();
+
     /**
      * Find the solution to F(X) = 0 by damped Newton iteration.
      */
@@ -75,32 +79,46 @@ class Newton
         m_max[n] = upper;
     }
 
-    void setConstant(size_t component, bool constant) {
-        m_constant[component] = constant;
+    void setConstants(vector_int constantComponents) {
+        m_constantComponents = constantComponents;
     }
 
     void evalJacobian(doublereal* x, doublereal* xdot);
 
+    void getSolution(double* x) {
+        for (size_t i = 0; i < m_nv; i++) {
+            x[i] = m_x[i];
+        }
+    }
+
 protected:
     FuncEval* m_residfunc;
 
-    DenseMatrix m_jacobian;
+    //! number of variables
+    size_t m_nv;
+
+    //! solution converged if [weightedNorm(sol, step) < m_convergenceThreshold]
+    doublereal m_convergenceThreshold;
+
+    DenseMatrix m_jacobian, m_jacFactored;
     int m_jacAge, m_jacMaxAge;
     doublereal m_jacRtol, m_jacAtol;
 
 
-    //! Work arrays of size #m_nv used in solve().
+    //! 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_rtol_ss, m_rtol_ts;
     vector_fp m_atol_ss, m_atol_ts;
 
-    //! number of variables
-    size_t m_nv;
+    vector_fp m_xlast, m_xsave;
+
+    //! the indexes of any constant variables
+    vector_int m_constantComponents;
 
-    //! constant variables
-    std::vector<bool> m_constant;
+    //! current timestep reciprocal
+    doublereal m_rdt = 0;
 };
 }
 

src/numerics/DenseMatrix.cpp

@@ -140,6 +140,83 @@ vector_int& DenseMatrix::ipiv()
     return m_ipiv;
 }
 
+int factor(DenseMatrix& A)
+{
+    int info = 0;
+    #if CT_USE_LAPACK
+        ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
+                  A.nRows(), &A.ipiv()[0], info);
+        if (info > 0) {
+            if (A.m_printLevel) {
+                writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d   U(i,i) is exactly zero. The factorization has"
+                          " been completed, but the factor U is exactly singular, and division by zero will occur if "
+                          "it is used to solve a system of equations.\n", info);
+            }
+            if (!A.m_useReturnErrorCode) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                                    "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has"
+                                    " been completed, but the factor U is exactly singular, and division by zero will occur if "
+                                    "it is used to solve a system of equations.", info);
+            }
+            return info;
+        } else if (info < 0) {
+            if (A.m_printLevel) {
+                writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
+            }
+
+            throw CanteraError("solve(DenseMatrix& A, double* b)",
+                               "DGETRF returned INFO = {}. The argument i has an illegal value", info);
+        }
+    #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
+    #endif
+    return info;
+}
+
+int solveFactored(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
+{
+    if (A.nColumns() != A.nRows()) {
+        if (A.m_printLevel) {
+            writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
+        }
+        throw CanteraError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix");
+    }
+
+    int info = 0;
+    if (ldb == 0) {
+        ldb = A.nColumns();
+    }
+    #if CT_USE_LAPACK
+        ct_dgetrs(ctlapack::NoTranspose, A.nRows(), nrhs, A.ptrColumn(0),
+                  A.nRows(), &A.ipiv()[0], b, ldb, info);
+        if (info != 0) {
+            if (A.m_printLevel) {
+                writelog("solve(DenseMatrix& A, double* b): DGETRS returned INFO = {}\n", info);
+            }
+            if (info < 0 || !A.m_useReturnErrorCode) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                                   "DGETRS returned INFO = {}", info);
+            }
+        }
+    #else
+        for (size_t i = 0; i < nrhs; i++) {
+            MappedVector bm(b + ldb*i, A.nColumns());
+            bm = lu.solve(bm);
+        }
+    #endif
+    return info;
+}
+
 int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
 {
     if (A.nColumns() != A.nRows()) {