Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
5 
9 
10 #if CT_USE_LAPACK
12 #else
13  #if CT_SUNDIALS_USE_LAPACK
14  #if CT_SUNDIALS_VERSION >= 30
15  #include "sunlinsol/sunlinsol_lapackband.h"
16  #else
17  #include "cvodes/cvodes_lapack.h"
18  #endif
19  #else
20  #if CT_SUNDIALS_VERSION >= 30
21  #include "sunlinsol/sunlinsol_band.h"
22  #else
23  #include "cvodes/cvodes_dense.h"
24  #include "cvodes/cvodes_band.h"
25  #endif
26  #endif
27 #endif
28 
29 #include <cstring>
30 #include <fstream>
31 
32 using namespace std;
33 
34 namespace Cantera
35 {
36 
37 // pImpl wrapper class for vector of Sundials index types to avoid needing to
38 // include Sundials headers in BandMatrix.h
39 struct BandMatrix::PivData {
40 #if CT_USE_LAPACK
41  vector_int data;
42 #elif CT_SUNDIALS_VERSION >= 30
43  std::vector<sunindextype> data;
44 #else
45  std::vector<long int> data;
46 #endif
47 };
48 
49 BandMatrix::BandMatrix() :
50  m_n(0),
51  m_kl(0),
52  m_ku(0),
53  m_zero(0.0),
54  m_ipiv{new PivData()},
55  m_info(0)
56 {
57 }
58 
59 BandMatrix::~BandMatrix()
60 {
61  // Needs to be defined here so m_ipiv can be deleted
62 }
63 
64 BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) :
65  m_n(n),
66  m_kl(kl),
67  m_ku(ku),
68  m_zero(0.0),
69  m_ipiv{new PivData()},
70  m_info(0)
71 {
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);
77  m_colPtrs.resize(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];
83  }
84 }
85 
86 BandMatrix::BandMatrix(const BandMatrix& y) :
87  GeneralMatrix(y),
88  m_n(0),
89  m_kl(0),
90  m_ku(0),
91  m_zero(0.0),
92  m_ipiv{new PivData()},
93  m_info(y.m_info)
94 {
95  m_n = y.m_n;
96  m_kl = y.m_kl;
97  m_ku = y.m_ku;
98  data = y.data;
99  ludata = y.ludata;
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];
107  }
108 }
109 
110 BandMatrix& BandMatrix::operator=(const BandMatrix& y)
111 {
112  if (&y == this) {
113  return *this;
114  }
115  GeneralMatrix::operator=(y);
116  m_n = y.m_n;
117  m_kl = y.m_kl;
118  m_ku = y.m_ku;
119  m_ipiv->data = y.m_ipiv->data;
120  data = y.data;
121  ludata = y.ludata;
122  m_colPtrs.resize(m_n);
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++) {
126  m_colPtrs[j] = &data[ldab * j];
127  m_lu_col_ptrs[j] = &ludata[ldab * j];
128  }
129  m_info = y.m_info;
130  return *this;
131 }
132 
133 void BandMatrix::resize(size_t n, size_t kl, size_t ku, doublereal v)
134 {
135  m_n = n;
136  m_kl = kl;
137  m_ku = ku;
138  data.resize(n*(2*kl + ku + 1));
139  ludata.resize(n*(2*kl + ku + 1));
140  m_ipiv->data.resize(m_n);
141  fill(data.begin(), data.end(), v);
142  m_colPtrs.resize(m_n);
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++) {
146  m_colPtrs[j] = &data[ldab * j];
147  m_lu_col_ptrs[j] = &ludata[ldab * j];
148  }
149  m_factored = false;
150 }
151 
152 void BandMatrix::bfill(doublereal v)
153 {
154  std::fill(data.begin(), data.end(), v);
155  m_factored = false;
156 }
157 
159 {
160  std::fill(data.begin(), data.end(), 0.0);
161  m_factored = false;
162 }
163 
164 doublereal& BandMatrix::operator()(size_t i, size_t j)
165 {
166  return value(i,j);
167 }
168 
169 doublereal BandMatrix::operator()(size_t i, size_t j) const
170 {
171  return value(i,j);
172 }
173 
174 doublereal& BandMatrix::value(size_t i, size_t j)
175 {
176  m_factored = false;
177  if (i + m_ku < j || i > j + m_kl) {
178  return m_zero;
179  }
180  return data[index(i,j)];
181 }
182 
183 doublereal BandMatrix::value(size_t i, size_t j) const
184 {
185  if (i + m_ku < j || i > j + m_kl) {
186  return 0.0;
187  }
188  return data[index(i,j)];
189 }
190 
191 size_t BandMatrix::index(size_t i, size_t j) const
192 {
193  return (2*m_kl + m_ku)*j + m_kl + m_ku + i;
194 }
195 
196 doublereal BandMatrix::_value(size_t i, size_t j) const
197 {
198  return data[index(i,j)];
199 }
200 
201 size_t BandMatrix::nRows() const
202 {
203  return m_n;
204 }
205 
206 size_t BandMatrix::nColumns() const
207 {
208  return m_n;
209 }
210 
212 {
213  return m_kl;
214 }
215 
217 {
218  return m_ku;
219 }
220 
221 size_t BandMatrix::ldim() const
222 {
223  return 2*m_kl + m_ku + 1;
224 }
225 
226 void BandMatrix::mult(const doublereal* b, doublereal* prod) const
227 {
228  for (size_t m = 0; m < m_n; m++) {
229  double sum = 0.0;
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];
234  }
235  prod[m] = sum;
236  }
237 }
238 
239 void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
240 {
241  for (size_t n = 0; n < m_n; n++) {
242  double sum = 0.0;
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];
247  }
248  prod[n] = sum;
249  }
250 }
251 
253 {
254  ludata = data;
255 #if CT_USE_LAPACK
256  ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(),
257  ludata.data(), ldim(), m_ipiv->data.data(), m_info);
258 #else
259  long int nu = static_cast<long int>(nSuperDiagonals());
260  long int nl = static_cast<long int>(nSubDiagonals());
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());
264 #endif
265  if (m_info != 0) {
266  throw Cantera::CanteraError("BandMatrix::factor",
267  "Factorization failed with DGBTRF error code {}.", m_info);
268  }
269  m_factored = true;
270  return m_info;
271 }
272 
273 int BandMatrix::solve(const doublereal* const b, doublereal* const x)
274 {
275  copy(b, b + m_n, x);
276  return solve(x);
277 }
278 
279 int BandMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
280 {
281  if (!m_factored) {
282  factor();
283  }
284  if (ldb == 0) {
285  ldb = nColumns();
286  }
287 #if CT_USE_LAPACK
288  ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
289  nSuperDiagonals(), nrhs, ludata.data(), ldim(),
290  m_ipiv->data.data(), b, ldb, m_info);
291 #else
292  long int nu = static_cast<long int>(nSuperDiagonals());
293  long int nl = static_cast<long int>(nSubDiagonals());
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);
297  m_info = 0;
298 #endif
299 
300  if (m_info != 0) {
301  throw Cantera::CanteraError("BandMatrix::solve",
302  "Linear solve failed with DGBTRS error code {}.", m_info);
303  }
304  return m_info;
305 }
306 
307 vector_fp::iterator BandMatrix::begin()
308 {
309  m_factored = false;
310  return data.begin();
311 }
312 
313 vector_fp::iterator BandMatrix::end()
314 {
315  m_factored = false;
316  return data.end();
317 }
318 
319 vector_fp::const_iterator BandMatrix::begin() const
320 {
321  return data.begin();
322 }
323 
324 vector_fp::const_iterator BandMatrix::end() const
325 {
326  return data.end();
327 }
328 
329 ostream& operator<<(ostream& s, const BandMatrix& m)
330 {
331  for (size_t i = 0; i < m.nRows(); i++) {
332  s << m(i, 0);
333  for (size_t j = 1; j < m.nColumns(); j++) {
334  s << ", " << m(i,j);
335  }
336  s << endl;
337  }
338  return s;
339 }
340 
341 doublereal BandMatrix::rcond(doublereal a1norm)
342 {
343  iwork_.resize(m_n);
344  work_.resize(3 * m_n);
345 
346  if (m_factored != 1) {
347  throw CanteraError("BandMatrix::rcond()", "matrix isn't factored correctly");
348  }
349 
350 #if CT_USE_LAPACK
351  size_t ldab = (2 *m_kl + m_ku + 1);
352  int rinfo = 0;
353  double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(),
354  ldab, m_ipiv->data.data(), a1norm, work_.data(), iwork_.data(), rinfo);
355  if (rinfo != 0) {
356  throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = {}", rinfo);
357  }
358  return rcond;
359 #else
360  throw CanteraError("BandMatrix::rcond", "not implemented when LAPACK is missing");
361 #endif
362 }
363 
365 {
366  return 0;
367 }
368 
369 doublereal BandMatrix::oneNorm() const
370 {
371  double value = 0.0;
372  for (size_t j = 0; j < m_n; j++) {
373  double sum = 0.0;
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));
378  }
379  value = std::max(sum, value);
380  }
381  return value;
382 }
383 
384 size_t BandMatrix::checkRows(doublereal& valueSmall) const
385 {
386  valueSmall = 1.0E300;
387  size_t iSmall = npos;
388  for (size_t i = 0; i < m_n; i++) {
389  double valueS = 0.0;
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);
394  }
395  if (valueS < valueSmall) {
396  iSmall = i;
397  valueSmall = valueS;
398  if (valueSmall == 0.0) {
399  return iSmall;
400  }
401  }
402  }
403  return iSmall;
404 }
405 
406 size_t BandMatrix::checkColumns(doublereal& valueSmall) const
407 {
408  valueSmall = 1.0E300;
409  size_t jSmall = npos;
410  for (size_t j = 0; j < m_n; j++) {
411  double valueS = 0.0;
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);
416  }
417  if (valueS < valueSmall) {
418  jSmall = j;
419  valueSmall = valueS;
420  if (valueSmall == 0.0) {
421  return jSmall;
422  }
423  }
424  }
425  return jSmall;
426 }
427 
428 doublereal* BandMatrix::ptrColumn(size_t j)
429 {
430  return m_colPtrs[j];
431 }
432 
433 doublereal* const* BandMatrix::colPts()
434 {
435  return &m_colPtrs[0];
436 }
437 
438 }
size_t m_kl
Number of subdiagonals of the matrix.
Definition: BandMatrix.h:274
virtual vector_fp::iterator begin()
Returns an iterator for the start of the band storage data.
Definition: BandMatrix.cpp:307
vector_int iwork_
Extra work array needed - size = n.
Definition: BandMatrix.h:292
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j.
Definition: BandMatrix.cpp:428
vector_fp work_
Extra dp work array needed - size = 3n.
Definition: BandMatrix.h:295
BandMatrix()
Base Constructor.
Definition: BandMatrix.cpp:49
doublereal _value(size_t i, size_t j) const
Return the value of the (i,j) element for (i,j) within the bandwidth.
Definition: BandMatrix.cpp:196
virtual size_t checkColumns(doublereal &valueSmall) const
Check to see if we have any zero columns in the Jacobian.
Definition: BandMatrix.cpp:406
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"
Definition: ct_defs.h:165
int m_factored
Indicates whether the matrix is factored.
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
Definition: BandMatrix.cpp:273
doublereal & operator()(size_t i, size_t j)
Index into the (i,j) element.
Definition: BandMatrix.cpp:164
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...
Definition: BandMatrix.cpp:191
virtual void leftMult(const doublereal *const b, doublereal *const prod) const
Multiply b*A and write result to prod.
Definition: BandMatrix.cpp:239
STL namespace.
vector_fp data
Matrix data.
Definition: BandMatrix.h:261
size_t m_n
Number of rows and columns of the matrix.
Definition: BandMatrix.h:271
size_t nSubDiagonals() const
Number of subdiagonals.
Definition: BandMatrix.cpp:211
void resize(size_t n, size_t kl, size_t ku, doublereal v=0.0)
Resize the matrix problem.
Definition: BandMatrix.cpp:133
size_t nColumns() const
Number of columns.
Definition: BandMatrix.cpp:206
virtual doublereal *const * colPts()
Return a vector of const pointers to the columns.
Definition: BandMatrix.cpp:433
virtual doublereal oneNorm() const
Returns the one norm of the matrix.
Definition: BandMatrix.cpp:369
virtual int factorAlgorithm() const
Returns the factor algorithm used.
Definition: BandMatrix.cpp:364
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
std::unique_ptr< PivData > m_ipiv
Pivot vector.
Definition: BandMatrix.h:282
size_t nSuperDiagonals() const
Number of superdiagonals.
Definition: BandMatrix.cpp:216
virtual doublereal rcond(doublereal a1norm)
Returns an estimate of the inverse of the condition number for the matrix.
Definition: BandMatrix.cpp:341
virtual size_t checkRows(doublereal &valueSmall) const
Check to see if we have any zero rows in the Jacobian.
Definition: BandMatrix.cpp:384
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Definition: BandMatrix.cpp:174
vector_fp ludata
Factorized data.
Definition: BandMatrix.h:268
size_t m_ku
Number of super diagonals of the matrix.
Definition: BandMatrix.h:277
virtual size_t nRows() const
Return the number of rows in the matrix.
Definition: BandMatrix.cpp:201
size_t ldim() const
Return the number of rows of storage needed for the band storage.
Definition: BandMatrix.cpp:221
void bfill(doublereal v=0.0)
Fill or zero the matrix.
Definition: BandMatrix.cpp:152
doublereal m_zero
value of zero
Definition: BandMatrix.h:280
Contains declarations for string manipulation functions within Cantera.
int factor()
Perform an LU decomposition, the LAPACK routine DGBTRF is used.
Definition: BandMatrix.cpp:252
virtual void zero()
Zero the matrix elements.
Definition: BandMatrix.cpp:158
std::vector< doublereal * > m_colPtrs
Vector of column pointers.
Definition: BandMatrix.h:288
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
virtual void mult(const doublereal *b, doublereal *prod) const
Multiply A*b and write result to prod.
Definition: BandMatrix.cpp:226
vector_fp::iterator end()
Returns an iterator for the end of the band storage data.
Definition: BandMatrix.cpp:313
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:34