13 #include "cantera/numerics/eigen_dense.h" 20 m_useReturnErrorCode(0),
27 m_useReturnErrorCode(0),
30 m_ipiv.resize(std::max(n, m));
33 for (
size_t j = 0; j < m; j++) {
41 m_useReturnErrorCode(0),
47 for (
size_t j = 0; j <
m_ncols; j++) {
53 DenseMatrix& DenseMatrix::operator=(
const DenseMatrix& y)
58 Array2D::operator=(y);
61 for (
size_t j = 0; j <
m_ncols; j++) {
72 m_ipiv.resize(std::max(n,m));
75 for (
size_t j = 0; j <
m_ncols; j++) {
81 doublereal*
const* DenseMatrix::colPts()
91 void DenseMatrix::mult(
const double* b,
double* prod)
const 94 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
95 static_cast<int>(
nRows()),
97 static_cast<int>(
nRows()), b, 1, 0.0, prod, 1);
100 MappedVector bm(const_cast<double*>(b),
nColumns());
101 MappedVector pm(prod,
nRows());
110 "Inner matrix dimensions do not agree: {} != {}",
nColumns(), B.
nRows());
114 "Output matrix has wrong dimensions: {}x{} != {}x{}",
118 doublereal*
const* prodcols = prod.colPts();
119 for (
size_t col=0; col < B.
nColumns(); ++col) {
122 mult(bcols[col], prodcols[col]);
128 for (
size_t n = 0; n <
nColumns(); n++) {
130 for (
size_t i = 0; i <
nRows(); i++) {
131 sum +=
value(i,n)*b[i];
146 writelogf(
"solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
148 throw CanteraError(
"solve(DenseMatrix& A, double* b)",
"Can only solve a square matrix");
160 writelogf(
"solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has" 161 " been completed, but the factor U is exactly singular, and division by zero will occur if " 162 "it is used to solve a system of equations.\n", info);
166 "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has" 167 " been completed, but the factor U is exactly singular, and division by zero will occur if " 168 "it is used to solve a system of equations.", info);
171 }
else if (info < 0) {
173 writelogf(
"solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
177 "DGETRF returned INFO = {}. The argument i has an illegal value", info);
184 writelogf(
"solve(DenseMatrix& A, double* b): DGETRS returned INFO = %d\n", info);
187 throw CanteraError(
"solve(DenseMatrix& A, double* b)",
"DGETRS returned INFO = {}", info);
193 auto lu = Am.partialPivLu();
195 auto lu = Am.fullPivLu();
196 if (lu.nonzeroPivots() <
static_cast<long int>(A.
nColumns())) {
198 "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}",
202 for (
size_t i = 0; i < nrhs; i++) {
203 MappedVector bm(b + ldb*i, A.
nColumns());
223 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
224 static_cast<int>(A.
nRows()), static_cast<int>(A.
nColumns()), 1.0,
227 MappedMatrix Am(&const_cast<DenseMatrix&>(A)(0,0), A.
nRows(), A.
nColumns());
228 MappedVector bm(const_cast<double*>(b), A.
nColumns());
229 MappedVector pm(prod, A.
nRows());
238 integer n =
static_cast<int>(nn !=
npos ? nn : A.
nRows());
243 writelogf(
"invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info);
246 throw CanteraError(
"invert(DenseMatrix& A, int nn)",
"DGETRS returned INFO = {}", info);
252 integer lwork =
static_cast<int>(work.size());
254 &A.
ipiv()[0], &work[0], lwork, info);
257 writelogf(
"invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info);
260 throw CanteraError(
"invert(DenseMatrix& A, int nn)",
"DGETRI returned INFO={}", info);
268 Am.topLeftCorner(nn, nn) = Am.topLeftCorner(nn, nn).inverse();
size_t nRows() const
Number of rows.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
vector_fp m_data
Data stored in a single array.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
const size_t npos
index returned by functions to indicate "no position"
vector_int m_ipiv
Vector of pivots. Length is equal to the max of m and n.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void increment(const DenseMatrix &A, const double *b, double *prod)
Multiply A*b and add it to the result in prod. Uses BLAS routine DGEMV.
int m_printLevel
Print Level.
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
std::vector< int > vector_int
Vector of ints.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
virtual void leftMult(const double *const b, double *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
std::vector< doublereal * > m_colPts
Vector of column pointers.
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
size_t m_ncols
Number of columns.
size_t m_nrows
Number of rows.
Base class for exceptions thrown by Cantera classes.
doublereal & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
const doublereal *const * const_colPts() const
Return a const vector of const pointers to the columns.
Contains declarations for string manipulation functions within Cantera.
DenseMatrix()
Default Constructor.
int m_useReturnErrorCode
Error Handling Flag.
Namespace for the Cantera kernel.
size_t nColumns() const
Number of columns.
vector_int & ipiv()
Return a changeable value of the pivot vector.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...