13 #if CT_SUNDIALS_USE_LAPACK 14 #if CT_SUNDIALS_VERSION >= 30 15 #include "sunlinsol/sunlinsol_lapackband.h" 17 #include "cvodes/cvodes_lapack.h" 20 #if CT_SUNDIALS_VERSION >= 30 21 #include "sunlinsol/sunlinsol_band.h" 23 #include "cvodes/cvodes_dense.h" 24 #include "cvodes/cvodes_band.h" 39 struct BandMatrix::PivData {
42 #elif CT_SUNDIALS_VERSION >= 30 43 std::vector<sunindextype> data;
45 std::vector<long int> data;
49 BandMatrix::BandMatrix() :
54 m_ipiv{
new PivData()},
59 BandMatrix::~BandMatrix()
69 m_ipiv{
new PivData()},
72 data.resize(n*(2*kl + ku + 1));
73 ludata.resize(n*(2*kl + ku + 1));
74 fill(data.begin(), data.end(), v);
75 fill(ludata.begin(), ludata.end(), 0.0);
76 m_ipiv->data.resize(m_n);
78 m_lu_col_ptrs.resize(n);
79 size_t ldab = (2*kl + ku + 1);
80 for (
size_t j = 0; j < n; j++) {
81 m_colPtrs[j] = &data[ldab * j];
82 m_lu_col_ptrs[j] = &ludata[ldab * j];
92 m_ipiv{
new PivData()},
100 m_ipiv->data = y.m_ipiv->data;
101 m_colPtrs.resize(m_n);
102 m_lu_col_ptrs.resize(m_n);
103 size_t ldab = (2 *m_kl + m_ku + 1);
104 for (
size_t j = 0; j < m_n; j++) {
105 m_colPtrs[j] = &data[ldab * j];
106 m_lu_col_ptrs[j] = &ludata[ldab * j];
110 BandMatrix& BandMatrix::operator=(
const BandMatrix& y)
115 GeneralMatrix::operator=(y);
119 m_ipiv->data = y.m_ipiv->data;
123 m_lu_col_ptrs.resize(
m_n);
124 size_t ldab = (2 *
m_kl +
m_ku + 1);
125 for (
size_t j = 0; j <
m_n; j++) {
127 m_lu_col_ptrs[j] = &
ludata[ldab * j];
138 data.resize(n*(2*kl + ku + 1));
139 ludata.resize(n*(2*kl + ku + 1));
143 m_lu_col_ptrs.resize(
m_n);
144 size_t ldab = (2 *
m_kl +
m_ku + 1);
145 for (
size_t j = 0; j < n; j++) {
147 m_lu_col_ptrs[j] = &
ludata[ldab * j];
154 std::fill(
data.begin(),
data.end(), v);
160 std::fill(
data.begin(),
data.end(), 0.0);
177 if (i + m_ku < j || i > j +
m_kl) {
185 if (i + m_ku < j || i > j +
m_kl) {
228 for (
size_t m = 0; m <
m_n; m++) {
230 size_t start = (m >=
m_kl) ? m -
m_kl : 0;
231 size_t stop = std::min(m +
m_ku + 1,
m_n);
232 for (
size_t j = start; j < stop; j++) {
233 sum +=
_value(m,j) * b[j];
241 for (
size_t n = 0; n <
m_n; n++) {
243 size_t start = (n >=
m_ku) ? n -
m_ku : 0;
244 size_t stop = std::min(n +
m_kl + 1,
m_n);
245 for (
size_t i = start; i < stop; i++) {
246 sum +=
_value(i,n) * b[i];
261 long int smu = nu + nl;
262 m_info = bandGBTRF(m_lu_col_ptrs.data(),
static_cast<long int>(
nColumns()),
263 nu, nl, smu,
m_ipiv->data.data());
267 "Factorization failed with DGBTRF error code {}.", m_info);
290 m_ipiv->data.data(), b, ldb, m_info);
294 long int smu = nu + nl;
295 double** a = m_lu_col_ptrs.data();
296 bandGBTRS(a, static_cast<long int>(
nColumns()), smu, nl,
m_ipiv->data.data(), b);
302 "Linear solve failed with DGBTRS error code {}.", m_info);
331 for (
size_t i = 0; i < m.
nRows(); i++) {
333 for (
size_t j = 1; j < m.nColumns(); j++) {
347 throw CanteraError(
"BandMatrix::rcond()",
"matrix isn't factored correctly");
356 throw CanteraError(
"BandMatrix::rcond()",
"DGBCON returned INFO = {}", rinfo);
360 throw CanteraError(
"BandMatrix::rcond",
"not implemented when LAPACK is missing");
372 for (
size_t j = 0; j <
m_n; j++) {
374 size_t start = (j >=
m_ku) ? j -
m_ku : 0;
375 size_t stop = std::min(j +
m_kl + 1,
m_n);
376 for (
size_t i = start; i < stop; i++) {
377 sum += std::abs(
_value(i,j));
386 valueSmall = 1.0E300;
387 size_t iSmall =
npos;
388 for (
size_t i = 0; i <
m_n; i++) {
390 size_t start = (i >
m_kl) ? i -
m_kl : 0;
391 size_t stop = std::min(i +
m_ku + 1,
m_n);
392 for (
size_t j = start; j < stop; j++) {
393 valueS = std::max(fabs(
_value(i,j)), valueS);
395 if (valueS < valueSmall) {
398 if (valueSmall == 0.0) {
408 valueSmall = 1.0E300;
409 size_t jSmall =
npos;
410 for (
size_t j = 0; j <
m_n; j++) {
412 size_t start = (j >
m_ku) ? j -
m_ku : 0;
413 size_t stop = std::min(j +
m_kl + 1,
m_n);
414 for (
size_t i = start; i < stop; i++) {
415 valueS = std::max(fabs(
_value(i,j)), valueS);
417 if (valueS < valueSmall) {
420 if (valueSmall == 0.0) {
size_t m_kl
Number of subdiagonals of the matrix.
virtual vector_fp::iterator begin()
Returns an iterator for the start of the band storage data.
vector_int iwork_
Extra work array needed - size = n.
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j.
vector_fp work_
Extra dp work array needed - size = 3n.
BandMatrix()
Base Constructor.
doublereal _value(size_t i, size_t j) const
Return the value of the (i,j) element for (i,j) within the bandwidth.
virtual size_t checkColumns(doublereal &valueSmall) const
Check to see if we have any zero columns in the Jacobian.
Various templated functions that carry out common vector operations (see Templated Utility Functions)...
const size_t npos
index returned by functions to indicate "no position"
int m_factored
Indicates whether the matrix is factored.
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
doublereal & operator()(size_t i, size_t j)
Index into the (i,j) element.
size_t index(size_t i, size_t j) const
Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array...
virtual void leftMult(const doublereal *const b, doublereal *const prod) const
Multiply b*A and write result to prod.
vector_fp data
Matrix data.
size_t m_n
Number of rows and columns of the matrix.
size_t nSubDiagonals() const
Number of subdiagonals.
void resize(size_t n, size_t kl, size_t ku, doublereal v=0.0)
Resize the matrix problem.
size_t nColumns() const
Number of columns.
virtual doublereal *const * colPts()
Return a vector of const pointers to the columns.
virtual doublereal oneNorm() const
Returns the one norm of the matrix.
virtual int factorAlgorithm() const
Returns the factor algorithm used.
std::vector< int > vector_int
Vector of ints.
std::unique_ptr< PivData > m_ipiv
Pivot vector.
size_t nSuperDiagonals() const
Number of superdiagonals.
virtual doublereal rcond(doublereal a1norm)
Returns an estimate of the inverse of the condition number for the matrix.
virtual size_t checkRows(doublereal &valueSmall) const
Check to see if we have any zero rows in the Jacobian.
Base class for exceptions thrown by Cantera classes.
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
vector_fp ludata
Factorized data.
size_t m_ku
Number of super diagonals of the matrix.
virtual size_t nRows() const
Return the number of rows in the matrix.
size_t ldim() const
Return the number of rows of storage needed for the band storage.
void bfill(doublereal v=0.0)
Fill or zero the matrix.
doublereal m_zero
value of zero
Contains declarations for string manipulation functions within Cantera.
int factor()
Perform an LU decomposition, the LAPACK routine DGBTRF is used.
virtual void zero()
Zero the matrix elements.
std::vector< doublereal * > m_colPtrs
Vector of column pointers.
Namespace for the Cantera kernel.
virtual void mult(const doublereal *b, doublereal *prod) const
Multiply A*b and write result to prod.
vector_fp::iterator end()
Returns an iterator for the end of the band storage data.
Declarations for the class BandMatrix which is a child class of GeneralMatrix for banded matrices han...
A class for banded matrices, involving matrix inversion processes.