14 #include "cantera/numerics/eigen_dense.h"
23 m_ipiv.resize(std::max(n, m));
26 for (
size_t j = 0; j < m; j++) {
38 for (
size_t j = 0; j <
m_ncols; j++) {
44DenseMatrix& DenseMatrix::operator=(
const DenseMatrix& y)
49 Array2D::operator=(y);
52 for (
size_t j = 0; j <
m_ncols; j++) {
63 m_ipiv.resize(std::max(n,m));
66 for (
size_t j = 0; j <
m_ncols; j++) {
72double*
const* DenseMatrix::colPts()
82void DenseMatrix::mult(
const double* b,
double* prod)
const
85 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
86 static_cast<int>(
nRows()),
88 static_cast<int>(
nRows()), b, 1, 0.0, prod, 1);
91 MappedVector bm(
const_cast<double*
>(b),
nColumns());
92 MappedVector pm(prod,
nRows());
101 "Inner matrix dimensions do not agree: {} != {}",
nColumns(), B.
nRows());
105 "Output matrix has wrong dimensions: {}x{} != {}x{}",
109 double*
const* prodcols = prod.colPts();
110 for (
size_t col=0; col < B.
nColumns(); ++col) {
113 mult(bcols[col], prodcols[col]);
119 for (
size_t n = 0; n <
nColumns(); n++) {
121 for (
size_t i = 0; i <
nRows(); i++) {
122 sum +=
value(i,n)*b[i];
137 writelogf(
"solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
139 throw CanteraError(
"solve(DenseMatrix& A, double* b)",
"Can only solve a square matrix");
151 writelogf(
"solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has"
152 " been completed, but the factor U is exactly singular, and division by zero will occur if "
153 "it is used to solve a system of equations.\n", info);
157 "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has"
158 " been completed, but the factor U is exactly singular, and division by zero will occur if "
159 "it is used to solve a system of equations.", info);
162 }
else if (info < 0) {
164 writelogf(
"solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
168 "DGETRF returned INFO = {}. The argument i has an illegal value", info);
175 writelog(
"solve(DenseMatrix& A, double* b): DGETRS returned INFO = {}\n", info);
179 "DGETRS returned INFO = {}", info);
185 auto lu = Am.partialPivLu();
187 auto lu = Am.fullPivLu();
188 if (lu.nonzeroPivots() <
static_cast<long int>(A.
nColumns())) {
190 "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}",
194 for (
size_t i = 0; i < nrhs; i++) {
195 MappedVector bm(b + ldb*i, A.
nColumns());
215 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
216 static_cast<int>(A.
nRows()),
static_cast<int>(A.
nColumns()), 1.0,
220 MappedVector bm(
const_cast<double*
>(b), A.
nColumns());
221 MappedVector pm(prod, A.
nRows());
230 integer n =
static_cast<int>(nn !=
npos ? nn : A.
nRows());
235 writelogf(
"invert(DenseMatrix& A, size_t nn): DGETRS returned INFO = %d\n", info);
238 throw CanteraError(
"invert(DenseMatrix& A, size_t nn)",
"DGETRS returned INFO = {}", info);
243 vector<double> work(n);
244 integer lwork =
static_cast<int>(work.size());
246 &A.
ipiv()[0], &work[0], lwork, info);
249 writelogf(
"invert(DenseMatrix& A, size_t nn): DGETRS returned INFO = %d\n", info);
252 throw CanteraError(
"invert(DenseMatrix& A, size_t nn)",
"DGETRI returned INFO={}", info);
260 Am.topLeftCorner(nn, nn) = Am.topLeftCorner(nn, nn).inverse();
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
vector< double > m_data
Data stored in a single array.
size_t m_nrows
Number of rows.
size_t nRows() const
Number of rows.
size_t m_ncols
Number of columns.
size_t nColumns() const
Number of columns.
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
virtual void leftMult(const double *const b, double *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
int m_printLevel
Print Level.
int m_useReturnErrorCode
Error Handling Flag.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
vector< int > & ipiv()
Return a changeable value of the pivot vector.
DenseMatrix()=default
Default Constructor.
const double *const * const_colPts() const
Return a const vector of const pointers to the columns.
vector< double * > m_colPts
Vector of column pointers.
vector< int > m_ipiv
Vector of pivots. Length is equal to the max of m and n.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Namespace for the Cantera kernel.
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.
const size_t npos
index returned by functions to indicate "no position"
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.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Contains declarations for string manipulation functions within Cantera.