From 155d7604b5976f92d92fdc852eec3e66122a4ff3 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Thu, 7 Apr 2011 16:14:00 +0000
Subject: [PATCH] Verified that the QR factorization works on one problem.

---
 Cantera/src/numerics/SquareMatrix.cpp | 37 ++++++++++++++------
 Cantera/src/numerics/SquareMatrix.h   | 49 +++++++++++++++++----------
 2 files changed, 58 insertions(+), 28 deletions(-)

diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp
index 8ecc08972..c8371d889 100644
--- a/Cantera/src/numerics/SquareMatrix.cpp
+++ b/Cantera/src/numerics/SquareMatrix.cpp
@@ -32,10 +32,12 @@ namespace Cantera {
   //====================================================================================================================
   SquareMatrix::SquareMatrix() :
     DenseMatrix(),
-    m_factored(0)
+    m_factored(0),
+    useQR_(0)
   {
   }
 
+  //====================================================================================================================
   // Constructor.
   /*
    * Create an \c n by \c n matrix, and initialize
@@ -46,21 +48,24 @@ namespace Cantera {
    */
   SquareMatrix::SquareMatrix(int n, doublereal v)  : 
     DenseMatrix(n, n, v),
-    m_factored(0)
+    m_factored(0),
+    useQR_(0)
   {
   }
-
+  //====================================================================================================================
   /*
    *
    * copy constructor
    */
   SquareMatrix::SquareMatrix(const SquareMatrix& y) :
     DenseMatrix(y),
-    m_factored(y.m_factored)
+    m_factored(y.m_factored),
+    useQR_(y.useQR_)
   {
   }
     
-  /**
+  //====================================================================================================================
+  /*
    * Assignment operator
    */
   SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) {
@@ -69,12 +74,18 @@ namespace Cantera {
     m_factored = y.m_factored;
     return *this;
   }
-    
-  /**
+  //====================================================================================================================
+  SquareMatrix::~SquareMatrix() {
+  }
+  //====================================================================================================================
+  /*
    * Solve Ax = b. Vector b is overwritten on exit with x.
    */
   int SquareMatrix::solve(double* b) 
   {
+    if (useQR_) {
+      return solveQR(b);
+    }
     int info=0;
     /*
      * Check to see whether the matrix has been factored.
@@ -122,10 +133,13 @@ namespace Cantera {
     DenseMatrix::resize(n, m, v);
   } 
   //====================================================================================================================
-  /**
+  /*
    * Factor A. A is overwritten with the LU decomposition of A.
    */
   int SquareMatrix::factor() {
+    if (useQR_) {
+      return factorQR();
+    }
     integer n = static_cast<int>(nRows());
     int info=0;
     m_factored = 1;
@@ -159,9 +173,10 @@ namespace Cantera {
   int SquareMatrix::factorQR() {
      if ((int) tau.size() < m_nrows)  {
        tau.resize(m_nrows, 0.0);
-       work.resize(9 * m_nrows, 0.0);
+       work.resize(8 * m_nrows, 0.0);
      }
      int info;
+     m_factored = 2;
      int lwork = work.size(); 
      ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info); 
      if (info != 0) {
@@ -173,7 +188,7 @@ namespace Cantera {
       }
      }
      int lworkOpt = work[0];
-     if (lworkOpt != lwork) {
+     if (lworkOpt > lwork) {
        work.resize(lworkOpt);
      }
      return info;
@@ -215,7 +230,7 @@ namespace Cantera {
       }
     }
     int lworkOpt = work[0];
-    if (lworkOpt <= lwork) {
+    if (lworkOpt > lwork) {
       work.resize(lworkOpt);
     }
 
diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h
index c8dac8895..1be99b62e 100644
--- a/Cantera/src/numerics/SquareMatrix.h
+++ b/Cantera/src/numerics/SquareMatrix.h
@@ -45,25 +45,35 @@ namespace Cantera {
      */
     SquareMatrix(int n, doublereal v = 0.0);
 
-    /**
-     * Copy Constructor
+    //! Copy Constructor
+    /*!
+     *  @param right Object to be copied 
      */
-    SquareMatrix(const SquareMatrix&);
-
-    /**
-     * Assignment operator
+    SquareMatrix(const SquareMatrix& right);
+   
+    //! Assignment operator
+    /*!
+     *  @param right  Object to be copied
      */
-    SquareMatrix& operator=(const SquareMatrix&);
+    SquareMatrix& operator=(const SquareMatrix& right);
 
 
-    /// Destructor. Does nothing.
-    virtual ~SquareMatrix(){}
+    //! Destructor. Does nothing.
+    virtual ~SquareMatrix();
 
-    /**
-     * Solves the Ax = b system returning x in the b spot.
+    
+     //! Solves the Ax = b system returning x in the b spot.
+    /*!
+     *  @param b  Vector for the rhs of the equation system
      */
-    int solve(double *b);
+    int solve(doublereal *b);
 
+    //! Resize the matrix
+    /*!
+     *  @param n Number of rows
+     *  @param m Number of columns
+     *  @param v double to fill the new space (defaults to zero)
+     */
     void resize(int n, int m, doublereal v = 0.0);
 
 
@@ -86,16 +96,14 @@ namespace Cantera {
      */
     int factorQR();
 
-    //! Solves the linear problem Ax=b using the QR algorithm
-    //! returning x in the b spot
+    //! 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
-     */
+    
+    //! clear the factored flag
     void clearFactorFlag();
     /**
      * set the factored flag
@@ -107,8 +115,15 @@ namespace Cantera {
      */
     int m_factored;
 
+    //! Work vector for QR algorithm
     vector_fp tau;
+
+    //! Work vector for QR algorithm
     vector_fp work;
+
+  public:
+    //!  Use the QR algorithm to factor and invert the matrix
+    int useQR_;
   };
 }