118 lines
2.5 KiB
C++
118 lines
2.5 KiB
C++
/**
|
|
* @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;
|
|
}
|
|
}
|
|
|