SquareMatrix.cpp

/**
 *  @file DenseMatrix.cpp
 */

/*
 * Copyright 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 "cantera/base/ct_defs.h"
#include "cantera/base/stringUtils.h"
#include "cantera/numerics/ctlapack.h"
#include "cantera/numerics/SquareMatrix.h"
#include "cantera/base/global.h"

#include <iostream>
#include <vector>
#include <cstring>

using namespace std;

namespace Cantera
{

//====================================================================================================================
SquareMatrix::SquareMatrix() :
    DenseMatrix(),
    GeneralMatrix(0),
    m_factored(0),
    a1norm_(0.0),
    useQR_(0)
{
}

//====================================================================================================================
// Constructor.
/*
 * Create an \c n by \c n matrix, and initialize
 * all elements to \c v.
 *
 * @param n   size of the square matrix
 * @param v   initial value of all matrix components.
 */
SquareMatrix::SquareMatrix(size_t n, doublereal v)  :
    DenseMatrix(n, n, v),
    GeneralMatrix(0),
    m_factored(0),
    a1norm_(0.0),
    useQR_(0)

{
}
//====================================================================================================================
/*
 *
 * copy constructor
 */
SquareMatrix::SquareMatrix(const SquareMatrix& y) :
    DenseMatrix(y),
    GeneralMatrix(0),
    m_factored(y.m_factored),
    a1norm_(y.a1norm_),
    useQR_(y.useQR_)
{
}

//====================================================================================================================
/*
 * Assignment operator
 */
SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y)
{
    if (&y == this) {
        return *this;
    }
    DenseMatrix::operator=(y);
    GeneralMatrix::operator=(y);
    m_factored = y.m_factored;
    a1norm_ = y.a1norm_;
    useQR_ = y.useQR_;
    return *this;
}
//====================================================================================================================
SquareMatrix::~SquareMatrix()
{
}
//====================================================================================================================
/*
 * Solve Ax = b. Vector b is overwritten on exit with x.
 */
int SquareMatrix::solve(doublereal* b)
{
    if (useQR_) {
        return solveQR(b);
    }
    int info=0;
    /*
     * Check to see whether the matrix has been factored.
     */
    if (!m_factored) {
        int retn = factor();
        if (retn) {
            return retn;
        }
    }
    /*
     * 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) {
        if (m_printLevel) {
            writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);
        }
        if (! m_useReturnErrorCode) {
            throw CELapackError("SquareMatrix::solve()", "DGETRS returned INFO = " + int2str(info));
        }
    }
    return info;
}
//====================================================================================================================
/*
 * Set all entries to zero
 */
void SquareMatrix::zero()
{
    size_t n = nRows();
    if (n > 0) {
        size_t 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));
    }
}
//====================================================================================================================
void SquareMatrix::resize(size_t n, size_t m, doublereal v)
{
    DenseMatrix::resize(n, m, v);
}

//====================================================================================================================
// Multiply A*b and write result to prod.
/*
 *  @param b    Vector to do the rh multiplication
 *  @param prod OUTPUT vector to receive the result
 */
void  SquareMatrix::mult(const doublereal* b, doublereal* prod) const
{
    DenseMatrix::mult(b, prod);
}
//====================================================================================================================
// Multiply b*A and write result to prod.
/*
 *  @param b    Vector to do the lh multiplication
 *  @param prod OUTPUT vector to receive the result
 */
void  SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
{
    DenseMatrix::leftMult(b, prod);
}
//====================================================================================================================
/*
 * Factor A. A is overwritten with the LU decomposition of A.
 */
int SquareMatrix::factor()
{
    if (useQR_) {
        return factorQR();
    }
    a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
    integer n = static_cast<int>(nRows());
    int info=0;
    m_factored = 1;
    ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()), DATA_PTR(ipiv()), info);
    if (info != 0) {
        if (m_printLevel) {
            writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info);
        }
        if (! m_useReturnErrorCode) {
            throw CELapackError("SquareMatrix::factor()", "DGETRS returned INFO = "+int2str(info));
        }
    }
    return info;
}
//=====================================================================================================================
/*
 * clear the factored flag
 */
void SquareMatrix::clearFactorFlag()
{
    m_factored = 0;
}
//=====================================================================================================================
/*
 * set the factored flag
 */
void SquareMatrix::setFactorFlag()
{
    m_factored = 1;
}
//=====================================================================================================================
int SquareMatrix::factorQR()
{
    if (tau.size() < m_nrows)  {
        tau.resize(m_nrows, 0.0);
        work.resize(8 * m_nrows, 0.0);
    }
    a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
    int info = 0;
    m_factored = 2;
    size_t lwork = work.size();
    ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
    if (info != 0) {
        if (m_printLevel) {
            writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
        }
        if (! m_useReturnErrorCode) {
            throw CELapackError("SquareMatrix::factorQR()", "DGEQRF returned INFO = " + int2str(info));
        }
    }
    size_t lworkOpt = static_cast<size_t>(work[0]);
    if (lworkOpt > lwork) {
        work.resize(lworkOpt);
    }


    return info;
}
//=====================================================================================================================
/*
 * Solve Ax = b. Vector b is overwritten on exit with x.
 */
int SquareMatrix::solveQR(doublereal* b)
{
    int info=0;
    /*
     * Check to see whether the matrix has been factored.
     */
    if (!m_factored) {
        int retn = factorQR();
        if (retn) {
            return retn;
        }
    }

    size_t 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));
        }
    }
    size_t lworkOpt = static_cast<size_t>(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;
}
//=====================================================================================================================
doublereal SquareMatrix::rcond(doublereal anorm)
{

    if (iwork_.size() < m_nrows) {
        iwork_.resize(m_nrows);
    }
    if (work.size() <4 * m_nrows) {
        work.resize(4 * m_nrows);
    }
    doublereal rcond = 0.0;
    if (m_factored != 1) {
        throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly");
    }

    //  doublereal anorm = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));


    int rinfo = 0;
    rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work),
                      DATA_PTR(iwork_), rinfo);
    if (rinfo != 0) {
        if (m_printLevel) {
            writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo);
        }
        if (! m_useReturnErrorCode) {
            throw CELapackError("SquareMatrix::rcond()", "DGECON returned INFO = " + int2str(rinfo));
        }
    }
    return rcond;
}
//=====================================================================================================================
doublereal SquareMatrix::oneNorm() const
{
    return a1norm_;
}
//=====================================================================================================================
doublereal SquareMatrix::rcondQR()
{

    if (iwork_.size() < m_nrows) {
        iwork_.resize(m_nrows);
    }
    if (work.size() <3 * m_nrows) {
        work.resize(3 * m_nrows);
    }
    doublereal rcond = 0.0;
    if (m_factored != 2) {
        throw CELapackError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
    }

    int rinfo = 0;
    rcond =  ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work),
                       DATA_PTR(iwork_), rinfo);
    if (rinfo != 0) {
        if (m_printLevel) {
            writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
        }
        if (! m_useReturnErrorCode) {
            throw CELapackError("SquareMatrix::rcondQR()", "DTRCON returned INFO = " + int2str(rinfo));
        }
    }
    return rcond;
}
//=====================================================================================================================
void SquareMatrix::useFactorAlgorithm(int fAlgorithm)
{
    useQR_ = fAlgorithm;
}
//=====================================================================================================================
int SquareMatrix::factorAlgorithm() const
{
    return (int) useQR_;
}
//=====================================================================================================================
bool SquareMatrix::factored() const
{
    return (m_factored != 0);
}
//=====================================================================================================================
// Return a pointer to the top of column j, columns are contiguous in memory
/*
 *  @param j   Value of the column
 *
 *  @return  Returns a pointer to the top of the column
 */
doublereal* SquareMatrix::ptrColumn(size_t j)
{
    return Array2D::ptrColumn(j);
}
//=====================================================================================================================
// Copy the data from one array into another without doing any checking
/*
 *  This differs from the assignment operator as no resizing is done and memcpy() is used.
 *  @param y Array to be copied
 */
void  SquareMatrix::copyData(const GeneralMatrix& y)
{
    const SquareMatrix* yy_ptr = dynamic_cast<const SquareMatrix*>(& y);
    Array2D::copyData(*yy_ptr);
}
//=====================================================================================================================
size_t  SquareMatrix::nRows() const
{
    return m_nrows;
}
//=====================================================================================================================
size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const
{
    return m_nrows;
}
//=====================================================================================================================
GeneralMatrix* SquareMatrix::duplMyselfAsGeneralMatrix() const
{
    SquareMatrix* dd = new SquareMatrix(*this);
    return static_cast<GeneralMatrix*>(dd);
}
//=====================================================================================================================
// Return an iterator pointing to the first element
vector_fp::iterator SquareMatrix::begin()
{
    return m_data.begin();
}
//=====================================================================================================================
// Return a const iterator pointing to the first element
vector_fp::const_iterator SquareMatrix::begin() const
{
    return m_data.begin();
}
//=====================================================================================================================
// Return a vector of const pointers to the columns
/*
 *  Note the value of the pointers are protected by their being const.
 *  However, the value of the matrix is open to being changed.
 *
 *   @return returns a vector of pointers to the top of the columns
 *           of the matrices.
 */
doublereal*   const* SquareMatrix::colPts()
{
    return DenseMatrix::colPts();
}
//=====================================================================================================================

size_t SquareMatrix::checkRows(doublereal& valueSmall) const
{
    valueSmall = 1.0E300;
    size_t iSmall = npos;
    for (size_t i = 0; i < m_nrows; i++) {
        double valueS = 0.0;
        for (size_t j = 0; j < m_nrows; j++) {
            if (fabs(value(i,j)) > valueS) {
                valueS = fabs(value(i,j));
            }
        }
        if (valueS < valueSmall) {
            iSmall = i;
            valueSmall = valueS;
        }
    }
    return iSmall;
}
//=====================================================================================================================
size_t SquareMatrix::checkColumns(doublereal& valueSmall) const
{
    valueSmall = 1.0E300;
    size_t jSmall = npos;
    for (size_t j = 0; j < m_nrows; j++) {
        double valueS = 0.0;
        for (size_t i = 0; i < m_nrows; i++) {
            if (fabs(value(i,j)) > valueS) {
                valueS = fabs(value(i,j));
            }
        }
        if (valueS < valueSmall) {
            jSmall = j;
            valueSmall = valueS;
        }
    }
    return jSmall;
}
//=====================================================================================================================


}