From 0a98f73828788df77ed42762dd529ec40c8dc143 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Tue, 27 Jan 2009 15:59:13 +0000
Subject: [PATCH] Added a new object SquareMatrix, that can solve Ax=b.

---
 Cantera/src/numerics/SquareMatrix.cpp | 118 ++++++++++++++++++++++++++
 Cantera/src/numerics/SquareMatrix.h   |  96 +++++++++++++++++++++
 2 files changed, 214 insertions(+)
 create mode 100644 Cantera/src/numerics/SquareMatrix.cpp
 create mode 100644 Cantera/src/numerics/SquareMatrix.h

diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp
new file mode 100644
index 000000000..8b34c083f
--- /dev/null
+++ b/Cantera/src/numerics/SquareMatrix.cpp
@@ -0,0 +1,118 @@
+/**
+ *  @file DenseMatrix.cpp
+ *
+ */
+/*
+ * $Revision$
+ * $Date$
+ */
+/*
+ * Copywrite 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
+ * retains certain rights in this software.
+ * See file License.txt for licensing information.
+ */
+
+#include "ct_defs.h"
+#include "stringUtils.h"
+#include "ctlapack.h"
+#include "SquareMatrix.h"
+
+#include <iostream>
+#include <vector>
+
+#include <cstring>
+
+using namespace std;
+
+namespace Cantera {
+  /**
+   *
+   * copy constructor
+   */
+  SquareMatrix::SquareMatrix(const SquareMatrix& y) :
+    DenseMatrix(y),
+    m_factored(y.m_factored)
+  {
+  }
+    
+  /**
+   * Assignment operator
+   */
+  SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) {
+    if (&y == this) return *this;
+    DenseMatrix::operator=(y);
+    m_factored = y.m_factored;
+    return *this;
+  }
+    
+  /**
+   * Solve Ax = b. Vector b is overwritten on exit with x.
+   */
+  int SquareMatrix::solve(double* b) 
+  {
+    int info=0;
+    /*
+     * Check to see whether the matrix has been factored.
+     */
+    if (!m_factored) {
+      factor();
+    }
+    /*
+     * Solve the factored system
+     */
+    ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
+	      1, &(*(begin())), static_cast<int>(nRows()), 
+	      DATA_PTR(ipiv()), b, static_cast<int>(nColumns()), info);
+    if (info != 0)
+      throw CanteraError("SquareMatrix::solve",
+			 "DGETRS returned INFO = "+int2str(info));
+    return 0;
+  }
+
+  /**
+   * Set all entries to zero
+   */
+  void SquareMatrix::zero() {
+    int n = static_cast<int>(nRows());
+    if (n > 0) {
+      int nn = n * n;
+      double *sm = &m_data[0];
+      /*
+       * Using memset is the fastest way to zero a contiguous
+       * section of memory.
+       */
+      (void) memset((void *) sm, 0, nn * sizeof(double));
+    }
+  }
+   
+  /**
+   * Factor A. A is overwritten with the LU decomposition of A.
+   */
+  int SquareMatrix::factor() {
+    integer n = static_cast<int>(nRows());
+    int info=0;
+    m_factored = true;
+    ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()),
+	      DATA_PTR(ipiv()), info);
+    if (info != 0) {
+      cout << "Singular matrix, info = " << info << endl;
+      throw CanteraError("invert",
+			 "DGETRF returned INFO="+int2str(info));
+    }
+    return 0;
+  }
+  /*
+   * clear the factored flag
+   */
+  void SquareMatrix::clearFactorFlag() {
+    m_factored = false;
+  }
+  /**
+   * set the factored flag
+   */
+  void SquareMatrix::setFactorFlag() {
+    m_factored = true;
+  }
+}
+
diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h
new file mode 100644
index 000000000..76bf19c00
--- /dev/null
+++ b/Cantera/src/numerics/SquareMatrix.h
@@ -0,0 +1,96 @@
+/**
+ *  @file SquareMatrix.h
+ * Dense, Square (not sparse) matrices.
+ */
+
+/*
+ *  $Date$
+ *  $Revision$
+ */
+/*
+ * Copywrite 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
+ * retains certain rights in this software.
+ * See file License.txt for licensing information.
+ */
+
+#ifndef CT_SQUAREMATRIX_H
+#define CT_SQUAREMATRIX_H
+
+#include "DenseMatrix.h"
+
+namespace Cantera { 
+
+  /**
+   *  A class for full (non-sparse) matrices with Fortran-compatible
+   *  data storage. Adds matrix operations to class Array2D.
+   */
+  class SquareMatrix: public DenseMatrix {
+
+  public:
+
+    SquareMatrix():
+      DenseMatrix(),
+      m_factored(false)
+    {
+    }
+
+    /** 
+     * Constructor. Create an \c n by \c n matrix, and initialize
+     * all elements to \c v.
+     */
+    SquareMatrix(int n, doublereal v = 0.0)  : 
+      DenseMatrix(n, n, v),
+      m_factored(false)
+    {
+    }
+
+    /**
+     * Copy Constructor
+     */
+    SquareMatrix(const SquareMatrix&);
+
+    /**
+     * Assignment operator
+     */
+    SquareMatrix& operator=(const SquareMatrix&);
+
+
+    /// Destructor. Does nothing.
+    virtual ~SquareMatrix(){}
+
+    /**
+     * Solves the Ax = b system returning x in the b spot.
+     */
+    int solve(double *b);
+
+    /**
+     * Zero the matrix
+     */
+    void zero();
+
+    /**
+     * Factors the A matrix, overwriting A. We flip m_factored
+     * boolean to indicate that the matrix is now A-1.
+     */
+    int factor();
+    /**
+     * clear the factored flag
+     */
+    void clearFactorFlag();
+    /**
+     * set the factored flag
+     */
+    void setFactorFlag();
+
+    /*
+     * the factor flag
+     */
+    bool m_factored;
+  };
+}
+
+#endif
+
+
+