From 68768e56e7834e6eb552a35026045d3ed35415b2 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Wed, 6 Apr 2011 16:26:32 +0000
Subject: [PATCH] Another iteration to get the QR solve up to speed

---
 Cantera/src/numerics/SquareMatrix.cpp | 56 +++++++++++++++++++++++++++
 Cantera/src/numerics/SquareMatrix.h   |  7 ++++
 Cantera/src/numerics/ctlapack.h       |  9 +++--
 3 files changed, 68 insertions(+), 4 deletions(-)

diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp
index d93f55cd7..8ecc08972 100644
--- a/Cantera/src/numerics/SquareMatrix.cpp
+++ b/Cantera/src/numerics/SquareMatrix.cpp
@@ -179,6 +179,62 @@ namespace Cantera {
      return info;
   }
   //=====================================================================================================================
+  /*
+   * Solve Ax = b. Vector b is overwritten on exit with x.
+   */
+  int SquareMatrix::solveQR(double* b)
+  {
+    int info=0;
+    /*
+     * Check to see whether the matrix has been factored.
+     */
+    if (!m_factored) {
+      int retn = factorQR();
+      if (retn) {
+        return retn;
+      }
+    }
+    
+    int lwork = work.size(); 
+    if (lwork < m_nrows) {
+      work.resize(8 * m_nrows, 0.0);
+      lwork = 8 * m_nrows;
+    }
+
+    /*
+     * Solve the factored system
+     */
+    ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), b, m_nrows, 
+                        DATA_PTR(work), lwork, info);
+    if (info != 0) {
+      if (m_printLevel) {
+        writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
+      }
+      if (! m_useReturnErrorCode) {
+        throw CELapackError("SquareMatrix::solveQR()", "DORMQR returned INFO = " + int2str(info));
+      }
+    }
+    int lworkOpt = work[0];
+    if (lworkOpt <= lwork) {
+      work.resize(lworkOpt);
+    }
+
+    char dd = 'N';
+
+    ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1,  &(*(begin())), m_nrows, b,
+              m_nrows, info);
+    if (info != 0) {
+      if (m_printLevel) {
+        writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
+      }
+      if (! m_useReturnErrorCode) {
+        throw CELapackError("SquareMatrix::solveQR()", "DTRTRS returned INFO = " + int2str(info));
+      }
+    }
+
+    return info;
+  }
+  //=====================================================================================================================
 
 
 
diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h
index 79f77ec08..c8dac8895 100644
--- a/Cantera/src/numerics/SquareMatrix.h
+++ b/Cantera/src/numerics/SquareMatrix.h
@@ -86,6 +86,13 @@ namespace Cantera {
      */
     int factorQR();
 
+    //! Solves the linear problem Ax=b using the QR algorithm
+    //! returning x in the b spot
+    /*!
+     *  @param b  RHS to be solved.
+     */
+    int solveQR(doublereal *b); 
+
     /**
      * clear the factored flag
      */
diff --git a/Cantera/src/numerics/ctlapack.h b/Cantera/src/numerics/ctlapack.h
index fc50dfcab..28afaa25d 100644
--- a/Cantera/src/numerics/ctlapack.h
+++ b/Cantera/src/numerics/ctlapack.h
@@ -316,23 +316,24 @@ namespace Cantera {
 
   inline void ct_dormqr(ctlapack::side_t rlside, ctlapack::transpose_t trans, int m,
                         int n, int k, doublereal* a, int lda, doublereal *tau, doublereal *c, int ldc, 
-                        doublereal *work, const integer *lwork, int &info) {
+                        doublereal *work, int lwork, int &info) {
     char side = left_right[rlside];
     char tr = no_yes[trans];
     integer f_m = m;
     integer f_n = n;
     integer f_k = k;
+    integer f_lwork = lwork;
     integer f_lda = lda;
     integer f_ldc = ldc;
     integer f_info = info;
 #ifdef NO_FTN_STRING_LEN_AT_END
-    _DORMQR_(&side, &tr, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info);
+    _DORMQR_(&side, &tr, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info);
 #else
     ftnlen trsize = 1;
 #ifdef LAPACK_FTN_STRING_LEN_AT_END
-    _DORMQR_(&side, &tr, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info, trsize, trsize);
+    _DORMQR_(&side, &tr, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info, trsize, trsize);
 #else
-    _DORMQR_(&side, trsize, &tr, trsize, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info);
+    _DORMQR_(&side, trsize, &tr, trsize, &f_m, &f_n,  &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info);
 #endif
 #endif
     info = f_info;