From 961ca06cd511e2c757348f1497c84e93474a390f Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sun, 10 Apr 2016 21:25:55 -0400
Subject: [PATCH] Eliminate the use of SquareMatrix in favor of DenseMatrix

---
 include/cantera/kinetics/solveSP.h         |  6 ++--
 include/cantera/transport/MultiTransport.h |  3 +-
 src/kinetics/solveSP.cpp                   | 37 ++--------------------
 src/transport/MultiTransport.cpp           | 24 ++------------
 4 files changed, 10 insertions(+), 60 deletions(-)

diff --git a/include/cantera/kinetics/solveSP.h b/include/cantera/kinetics/solveSP.h
index 9870f2857..b2b0ec3f4 100644
--- a/include/cantera/kinetics/solveSP.h
+++ b/include/cantera/kinetics/solveSP.h
@@ -13,7 +13,7 @@
 #define SOLVESP_H
 
 #include "cantera/kinetics/InterfaceKinetics.h"
-#include "cantera/numerics/SquareMatrix.h"
+#include "cantera/numerics/DenseMatrix.h"
 
 //! @defgroup solvesp_methods Surface Problem Solver Methods
 //! @{
@@ -301,7 +301,7 @@ private:
      *  @param do_time Calculate a time dependent residual
      *  @param deltaT  Delta time for time dependent problem.
      */
-    void resjac_eval(SquareMatrix& jac, doublereal* resid,
+    void resjac_eval(DenseMatrix& jac, doublereal* resid,
                      doublereal* CSolnSP,
                      const doublereal* CSolnSPOld, const bool do_time,
                      const doublereal deltaT);
@@ -504,7 +504,7 @@ private:
 
     //! Jacobian. m_neq by m_neq computed Jacobian matrix for the local
     //! Newton's method.
-    SquareMatrix m_Jac;
+    DenseMatrix m_Jac;
 
 public:
     int m_ioflag;
diff --git a/include/cantera/transport/MultiTransport.h b/include/cantera/transport/MultiTransport.h
index 186d3d2e2..dc6caed96 100644
--- a/include/cantera/transport/MultiTransport.h
+++ b/include/cantera/transport/MultiTransport.h
@@ -10,7 +10,6 @@
 
 // Cantera includes
 #include "GasTransport.h"
-#include "cantera/numerics/SquareMatrix.h"
 
 namespace Cantera
 {
@@ -151,7 +150,7 @@ protected:
 
     // L matrix quantities
     DenseMatrix m_Lmatrix;
-    SquareMatrix m_aa;
+    DenseMatrix m_aa;
     vector_fp m_a;
     vector_fp m_b;
 
diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp
index 965d01c5f..ddb24682e 100644
--- a/src/kinetics/solveSP.cpp
+++ b/src/kinetics/solveSP.cpp
@@ -223,39 +223,8 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
         // Find the weighted norm of the residual
         double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
 
-        // Solve Linear system.  The solution is in resid[]
-        int info = m_Jac.factor();
-        if (info==0) {
-            m_Jac.solve(&m_resid[0]);
-        } else {
-            // Force convergence if residual is small to avoid "nan" results
-            // from the linear solve.
-            if (m_ioflag) {
-                writelogf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
-                          resid_norm, info);
-            }
-            for (size_t jcol = 0; jcol < m_neq; jcol++) {
-                m_resid[jcol] = 0.0;
-            }
-
-            // print out some helpful info
-            if (m_ioflag > 1) {
-                writelog("-----\n");
-                writelogf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
-                          iter,t_real, 1.0/inv_t);
-                writelog("solveSurfProb: init guess, current concentration,"
-                         "and prod rate:\n");
-                for (size_t jcol = 0; jcol < m_neq; jcol++) {
-                    writelog("\t%d  %g %g %g\n", jcol,
-                           m_CSolnSPInit[jcol], m_CSolnSP[jcol],
-                           m_netProductionRatesSave[m_kinSpecIndex[jcol]]);
-                }
-                writelog("-----\n");
-            }
-            if (do_time) {
-                t_real += time_scale;
-            }
-        }
+        // Solve Linear system.  The solution is in m_resid
+        solve(m_Jac, m_resid.data());
 
         // Calculate the Damping factor needed to keep all unknowns between 0
         // and 1, and not allow too large a change (factor of 2) in any unknown.
@@ -476,7 +445,7 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln,
     }
 }
 
-void solveSP::resjac_eval(SquareMatrix& jac,
+void solveSP::resjac_eval(DenseMatrix& jac,
                           doublereal resid[], doublereal CSoln[],
                           const doublereal CSolnOld[], const bool do_time,
                           const doublereal deltaT)
diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp
index 87465dfcf..604cbf046 100644
--- a/src/transport/MultiTransport.cpp
+++ b/src/transport/MultiTransport.cpp
@@ -226,17 +226,8 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_
         fluxes[jmax + n*ldf] = 0.0;
     }
 
-    // use LAPACK to solve the equations
-    int info = m_aa.factor();
-    if (info) {
-        throw CanteraError("MultiTransport::getSpeciesFluxes",
-                           "Error factorizing matrix.");
-    }
-    info = m_aa.solve(fluxes, ndim, ldf);
-    if (info) {
-        throw CanteraError("MultiTransport::getSpeciesFluxes",
-                           "Error solving linear system.");
-    }
+    // solve the equations
+    solve(m_aa, fluxes, ndim, ldf);
 
     size_t offset;
     doublereal pp = pressure_ig();
@@ -332,16 +323,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s
     fluxes[jmax] = 0.0;
 
     // Solve the equations
-    int info = m_aa.factor();
-    if (info) {
-        throw CanteraError("MultiTransport::getMassFluxes",
-                           "Error in factorization.  Info = {}", info);
-    }
-    info = m_aa.solve(fluxes);
-    if (info) {
-        throw CanteraError("MultiTransport::getMassFluxes",
-                           "Error in linear solve. Info = {}", info);
-    }
+    solve(m_aa, fluxes);
 
     doublereal pp = pressure_ig();
     // multiply diffusion velocities by rho * Y_k to create