Cantera  3.1.0a1
Loading...
Searching...
No Matches
BandMatrix.cpp
Go to the documentation of this file.
1//! @file BandMatrix.cpp Banded matrices.
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
9
10#if CT_USE_LAPACK
12#else
13 #include "sunlinsol/sunlinsol_band.h"
14#endif
15
16#include <cstring>
17#include <fstream>
18
19using namespace std;
20
21namespace Cantera
22{
23
24// pImpl wrapper class for vector of Sundials index types to avoid needing to
25// include Sundials headers in BandMatrix.h
27#if CT_USE_LAPACK
28 vector<int> data;
29#else
30 vector<sunindextype> data;
31#endif
32};
33
35 m_ipiv{new PivData()}
36{
37}
38
39BandMatrix::~BandMatrix()
40{
41 // Needs to be defined here so m_ipiv can be deleted
42}
43
44BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, double v) :
45 m_n(n),
46 m_kl(kl),
47 m_ku(ku),
48 m_ipiv{new PivData()}
49{
50 data.resize(n*(2*kl + ku + 1));
51 ludata.resize(n*(2*kl + ku + 1));
52 fill(data.begin(), data.end(), v);
53 fill(ludata.begin(), ludata.end(), 0.0);
54 m_ipiv->data.resize(m_n);
55 m_colPtrs.resize(n);
56 m_lu_col_ptrs.resize(n);
57 size_t ldab = (2*kl + ku + 1);
58 for (size_t j = 0; j < n; j++) {
59 m_colPtrs[j] = &data[ldab * j];
60 m_lu_col_ptrs[j] = &ludata[ldab * j];
61 }
62}
63
66 data(y.data),
67 ludata(y.ludata),
68 m_n(y.m_n),
69 m_kl(y.m_kl),
70 m_ku(y.m_ku),
71 m_ipiv{new PivData()},
72 m_info(y.m_info)
73{
74 m_ipiv->data = y.m_ipiv->data;
75 m_colPtrs.resize(m_n);
76 m_lu_col_ptrs.resize(m_n);
77 size_t ldab = (2 *m_kl + m_ku + 1);
78 for (size_t j = 0; j < m_n; j++) {
79 m_colPtrs[j] = &data[ldab * j];
80 m_lu_col_ptrs[j] = &ludata[ldab * j];
81 }
82}
83
84BandMatrix& BandMatrix::operator=(const BandMatrix& y)
85{
86 if (&y == this) {
87 return *this;
88 }
89 GeneralMatrix::operator=(y);
90 m_n = y.m_n;
91 m_kl = y.m_kl;
92 m_ku = y.m_ku;
93 m_ipiv->data = y.m_ipiv->data;
94 data = y.data;
95 ludata = y.ludata;
96 m_colPtrs.resize(m_n);
97 m_lu_col_ptrs.resize(m_n);
98 size_t ldab = (2 * m_kl + m_ku + 1);
99 for (size_t j = 0; j < m_n; j++) {
100 m_colPtrs[j] = &data[ldab * j];
101 m_lu_col_ptrs[j] = &ludata[ldab * j];
102 }
103 m_info = y.m_info;
104 return *this;
105}
106
107void BandMatrix::resize(size_t n, size_t kl, size_t ku, double v)
108{
109 m_n = n;
110 m_kl = kl;
111 m_ku = ku;
112 data.resize(n*(2*kl + ku + 1));
113 ludata.resize(n*(2*kl + ku + 1));
114 m_ipiv->data.resize(m_n);
115 fill(data.begin(), data.end(), v);
116 m_colPtrs.resize(m_n);
117 m_lu_col_ptrs.resize(m_n);
118 size_t ldab = (2 * m_kl + m_ku + 1);
119 for (size_t j = 0; j < n; j++) {
120 m_colPtrs[j] = &data[ldab * j];
121 m_lu_col_ptrs[j] = &ludata[ldab * j];
122 }
123 m_factored = false;
124}
125
126void BandMatrix::bfill(double v)
127{
128 std::fill(data.begin(), data.end(), v);
129 m_factored = false;
130}
131
133{
134 std::fill(data.begin(), data.end(), 0.0);
135 m_factored = false;
136}
137
138double& BandMatrix::operator()(size_t i, size_t j)
139{
140 return value(i,j);
141}
142
143double BandMatrix::operator()(size_t i, size_t j) const
144{
145 return value(i,j);
146}
147
148double& BandMatrix::value(size_t i, size_t j)
149{
150 m_factored = false;
151 if (i + m_ku < j || i > j + m_kl) {
152 return m_zero;
153 }
154 return data[index(i,j)];
155}
156
157double BandMatrix::value(size_t i, size_t j) const
158{
159 if (i + m_ku < j || i > j + m_kl) {
160 return 0.0;
161 }
162 return data[index(i,j)];
163}
164
165size_t BandMatrix::index(size_t i, size_t j) const
166{
167 return (2*m_kl + m_ku)*j + m_kl + m_ku + i;
168}
169
170double BandMatrix::_value(size_t i, size_t j) const
171{
172 return data[index(i,j)];
173}
174
175size_t BandMatrix::nRows() const
176{
177 return m_n;
178}
179
181{
182 return m_n;
183}
184
186{
187 return m_kl;
188}
189
191{
192 return m_ku;
193}
194
195size_t BandMatrix::ldim() const
196{
197 return 2*m_kl + m_ku + 1;
198}
199
200void BandMatrix::mult(const double* b, double* prod) const
201{
202 for (size_t m = 0; m < m_n; m++) {
203 double sum = 0.0;
204 size_t start = (m >= m_kl) ? m - m_kl : 0;
205 size_t stop = std::min(m + m_ku + 1, m_n);
206 for (size_t j = start; j < stop; j++) {
207 sum += _value(m,j) * b[j];
208 }
209 prod[m] = sum;
210 }
211}
212
213void BandMatrix::leftMult(const double* const b, double* const prod) const
214{
215 for (size_t n = 0; n < m_n; n++) {
216 double sum = 0.0;
217 size_t start = (n >= m_ku) ? n - m_ku : 0;
218 size_t stop = std::min(n + m_kl + 1, m_n);
219 for (size_t i = start; i < stop; i++) {
220 sum += _value(i,n) * b[i];
221 }
222 prod[n] = sum;
223 }
224}
225
227{
228 ludata = data;
229#if CT_USE_LAPACK
230 ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(),
231 ludata.data(), ldim(), m_ipiv->data.data(), m_info);
232#else
233 long int nu = static_cast<long int>(nSuperDiagonals());
234 long int nl = static_cast<long int>(nSubDiagonals());
235 long int smu = nu + nl;
236 #if CT_SUNDIALS_VERSION >= 60
237 m_info = SUNDlsMat_bandGBTRF(m_lu_col_ptrs.data(),
238 static_cast<long int>(nColumns()),
239 nu, nl, smu, m_ipiv->data.data());
240 #else
241 m_info = bandGBTRF(m_lu_col_ptrs.data(), static_cast<long int>(nColumns()),
242 nu, nl, smu, m_ipiv->data.data());
243 #endif
244#endif
245 if (m_info != 0) {
246 throw Cantera::CanteraError("BandMatrix::factor",
247 "Factorization failed with DGBTRF error code {}.", m_info);
248 }
249 m_factored = true;
250 return m_info;
251}
252
253int BandMatrix::solve(const double* const b, double* const x)
254{
255 copy(b, b + m_n, x);
256 return solve(x);
257}
258
259int BandMatrix::solve(double* b, size_t nrhs, size_t ldb)
260{
261 if (!m_factored) {
262 factor();
263 }
264 if (ldb == 0) {
265 ldb = nColumns();
266 }
267#if CT_USE_LAPACK
268 ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
269 nSuperDiagonals(), nrhs, ludata.data(), ldim(),
270 m_ipiv->data.data(), b, ldb, m_info);
271#else
272 long int nu = static_cast<long int>(nSuperDiagonals());
273 long int nl = static_cast<long int>(nSubDiagonals());
274 long int smu = nu + nl;
275 double** a = m_lu_col_ptrs.data();
276 #if CT_SUNDIALS_VERSION >= 60
277 SUNDlsMat_bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl,
278 m_ipiv->data.data(), b);
279 #else
280 bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl,
281 m_ipiv->data.data(), b);
282 #endif
283 m_info = 0;
284#endif
285
286 if (m_info != 0) {
287 throw Cantera::CanteraError("BandMatrix::solve",
288 "Linear solve failed with DGBTRS error code {}.", m_info);
289 }
290 return m_info;
291}
292
293ostream& operator<<(ostream& s, const BandMatrix& m)
294{
295 for (size_t i = 0; i < m.nRows(); i++) {
296 s << m(i, 0);
297 for (size_t j = 1; j < m.nColumns(); j++) {
298 s << ", " << m(i,j);
299 }
300 s << endl;
301 }
302 return s;
303}
304
305double BandMatrix::rcond(double a1norm)
306{
307 iwork_.resize(m_n);
308 work_.resize(3 * m_n);
309
310 if (m_factored != 1) {
311 throw CanteraError("BandMatrix::rcond", "matrix isn't factored correctly");
312 }
313
314#if CT_USE_LAPACK
315 size_t ldab = (2 *m_kl + m_ku + 1);
316 int rinfo = 0;
317 double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(),
318 ldab, m_ipiv->data.data(), a1norm, work_.data(), iwork_.data(), rinfo);
319 if (rinfo != 0) {
320 throw CanteraError("BandMatrix::rcond", "DGBCON returned INFO = {}", rinfo);
321 }
322 return rcond;
323#else
324 throw CanteraError("BandMatrix::rcond", "not implemented when LAPACK is missing");
325#endif
326}
327
329{
330 return 0;
331}
332
334{
335 double value = 0.0;
336 for (size_t j = 0; j < m_n; j++) {
337 double sum = 0.0;
338 size_t start = (j >= m_ku) ? j - m_ku : 0;
339 size_t stop = std::min(j + m_kl + 1, m_n);
340 for (size_t i = start; i < stop; i++) {
341 sum += std::abs(_value(i,j));
342 }
343 value = std::max(sum, value);
344 }
345 return value;
346}
347
348size_t BandMatrix::checkRows(double& valueSmall) const
349{
350 valueSmall = 1.0E300;
351 size_t iSmall = npos;
352 for (size_t i = 0; i < m_n; i++) {
353 double valueS = 0.0;
354 size_t start = (i > m_kl) ? i - m_kl : 0;
355 size_t stop = std::min(i + m_ku + 1, m_n);
356 for (size_t j = start; j < stop; j++) {
357 valueS = std::max(fabs(_value(i,j)), valueS);
358 }
359 if (valueS < valueSmall) {
360 iSmall = i;
361 valueSmall = valueS;
362 if (valueSmall == 0.0) {
363 return iSmall;
364 }
365 }
366 }
367 return iSmall;
368}
369
370size_t BandMatrix::checkColumns(double& valueSmall) const
371{
372 valueSmall = 1.0E300;
373 size_t jSmall = npos;
374 for (size_t j = 0; j < m_n; j++) {
375 double valueS = 0.0;
376 size_t start = (j > m_ku) ? j - m_ku : 0;
377 size_t stop = std::min(j + m_kl + 1, m_n);
378 for (size_t i = start; i < stop; i++) {
379 valueS = std::max(fabs(_value(i,j)), valueS);
380 }
381 if (valueS < valueSmall) {
382 jSmall = j;
383 valueSmall = valueS;
384 if (valueSmall == 0.0) {
385 return jSmall;
386 }
387 }
388 }
389 return jSmall;
390}
391
392double* BandMatrix::ptrColumn(size_t j)
393{
394 return m_colPtrs[j];
395}
396
397double* const* BandMatrix::colPts()
398{
399 return &m_colPtrs[0];
400}
401
402}
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.
Definition BandMatrix.h:37
void leftMult(const double *const b, double *const prod) const override
Multiply b*A and write result to prod.
double *const * colPts() override
Return a vector of const pointers to the columns.
double & operator()(size_t i, size_t j) override
Index into the (i,j) element.
vector< double > ludata
Factorized data.
Definition BandMatrix.h:246
size_t m_ku
Number of super diagonals of the matrix.
Definition BandMatrix.h:255
size_t m_kl
Number of subdiagonals of the matrix.
Definition BandMatrix.h:252
double m_zero
value of zero
Definition BandMatrix.h:258
double rcond(double a1norm) override
Returns an estimate of the inverse of the condition number for the matrix.
size_t nSubDiagonals() const
Number of subdiagonals.
void resize(size_t n, size_t kl, size_t ku, double v=0.0)
Resize the matrix problem.
size_t ldim() const
Return the number of rows of storage needed for the band storage.
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.
void bfill(double v=0.0)
Fill or zero the matrix.
BandMatrix()
Base Constructor.
size_t checkColumns(double &valueSmall) const override
Check to see if we have any zero columns in the Jacobian.
vector< double * > m_colPtrs
Vector of column pointers.
Definition BandMatrix.h:266
double * ptrColumn(size_t j) override
Return a pointer to the top of column j.
vector< double > work_
Extra dp work array needed - size = 3n.
Definition BandMatrix.h:273
int factor() override
Perform an LU decomposition, the LAPACK routine DGBTRF is used.
unique_ptr< PivData > m_ipiv
Pivot vector.
Definition BandMatrix.h:263
size_t nColumns() const
Number of columns.
void mult(const double *b, double *prod) const override
Multiply A*b and write result to prod.
size_t checkRows(double &valueSmall) const override
Check to see if we have any zero rows in the Jacobian.
size_t nSuperDiagonals() const
Number of superdiagonals.
size_t nRows() const override
Return the number of rows in the matrix.
vector< int > iwork_
Extra work array needed - size = n.
Definition BandMatrix.h:270
vector< double > data
Matrix data.
Definition BandMatrix.h:243
int solve(const double *const b, double *const x)
Solve the matrix problem Ax = b.
void zero() override
Zero the matrix elements.
size_t m_n
Number of rows and columns of the matrix.
Definition BandMatrix.h:249
double oneNorm() const override
Returns the one norm of the matrix.
double _value(size_t i, size_t j) const
Return the value of the (i,j) element for (i,j) within the bandwidth.
double & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
int factorAlgorithm() const override
Returns the factor algorithm used.
Base class for exceptions thrown by Cantera classes.
int m_factored
Indicates whether the matrix is factored.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
std::ostream & operator<<(std::ostream &s, const Array2D &m)
Output the current contents of the Array2D object.
Definition Array.cpp:100
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...