Cantera  2.0
DenseMatrix.cpp
Go to the documentation of this file.
1 /**
2  * @file DenseMatrix.cpp
3  */
4 
5 // Copyright 2001 California Institute of Technology
6 
7 #include "cantera/base/ct_defs.h"
8 #include "cantera/numerics/ctlapack.h"
9 #include "cantera/base/utilities.h"
10 #include "cantera/numerics/DenseMatrix.h"
11 #include "cantera/base/stringUtils.h"
12 #include "cantera/base/global.h"
13 
14 namespace Cantera
15 {
16 //====================================================================================================================
17 // Default Constructor
18 DenseMatrix::DenseMatrix() :
19  Array2D(0,0,0.0),
20  m_ipiv(0),
21  m_useReturnErrorCode(0),
22  m_printLevel(0)
23 {
24 }
25 //====================================================================================================================
26 /*
27  * Constructor. Create an \c n by \c m matrix, and initialize
28  * all elements to \c v.
29  */
30 DenseMatrix::DenseMatrix(size_t n, size_t m, doublereal v) :
31  Array2D(n, m, v),
32  m_ipiv(0),
33  m_useReturnErrorCode(0),
34  m_printLevel(0)
35 {
36  m_ipiv.resize(std::max(n, m));
37  m_colPts.resize(m);
38  if (!m_data.empty()) {
39  for (size_t j = 0; j < m; j++) {
40  m_colPts[j] = &(m_data[m_nrows*j]);
41  }
42  }
43 }
44 //====================================================================================================================
45 // Copy constructor
46 /*
47  * @param y Object to be copied
48  */
49 DenseMatrix::DenseMatrix(const DenseMatrix& y) :
50  Array2D(y),
51  m_ipiv(0),
52  m_useReturnErrorCode(0),
53  m_printLevel(0)
54 {
55  m_ipiv = y.ipiv();
56  m_colPts.resize(m_ncols);
57  if (!m_data.empty()) {
58  for (size_t j = 0; j < m_ncols; j++) {
59  m_colPts[j] = &(m_data[m_nrows*j]);
60  }
61  }
62 }
63 //====================================================================================================================
64 // assignment
65 DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y)
66 {
67  if (&y == this) {
68  return *this;
69  }
70  Array2D::operator=(y);
71  m_ipiv = y.ipiv();
72  m_colPts.resize(m_ncols);
73  for (size_t j = 0; j < m_ncols; j++) {
74  m_colPts[j] = &(m_data[m_nrows*j]);
75  }
76  m_useReturnErrorCode = y.m_useReturnErrorCode;
77  m_printLevel = y.m_printLevel;
78  return *this;
79 }
80 //====================================================================================================================
81 // Destructor. Does nothing.
82 DenseMatrix::~DenseMatrix()
83 {
84 }
85 //====================================================================================================================
86 void DenseMatrix::resize(size_t n, size_t m, doublereal v)
87 {
88  Array2D::resize(n,m,v);
89  m_ipiv.resize(std::max(n,m));
90  m_colPts.resize(m_ncols);
91  if (!m_data.empty()) {
92  for (size_t j = 0; j < m_ncols; j++) {
93  m_colPts[j] = &(m_data[m_nrows*j]);
94  }
95  }
96 }
97 //====================================================================================================================
98 doublereal* const* DenseMatrix::colPts()
99 {
100  return &(m_colPts[0]);
101 }
102 //====================================================================================================================
103 const doublereal* const* DenseMatrix::const_colPts() const
104 {
105  return &(m_colPts[0]);
106 }
107 //====================================================================================================================
108 void DenseMatrix::mult(const double* b, double* prod) const
109 {
110  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
111  static_cast<int>(nRows()),
112  static_cast<int>(nRows()), 1.0, ptrColumn(0), //begin(),
113  static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
114 }
115 //====================================================================================================================
116 void DenseMatrix::leftMult(const double* const b, double* const prod) const
117 {
118  size_t nc = nColumns();
119  size_t nr = nRows();
120  double sum = 0.0;
121  for (size_t n = 0; n < nc; n++) {
122  sum = 0.0;
123  for (size_t i = 0; i < nr; i++) {
124  sum += value(i,n)*b[i];
125  }
126  prod[n] = sum;
127  }
128 }
129 //====================================================================================================================
130 vector_int& DenseMatrix::ipiv()
131 {
132  return m_ipiv;
133 }
134 //====================================================================================================================
135 int solve(DenseMatrix& A, double* b)
136 {
137  int info = 0;
138  if (A.nColumns() != A.nRows()) {
139  if (A.m_printLevel) {
140  writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
141  }
142  throw CELapackError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix");
143  }
144  ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
145  A.nRows(), &A.ipiv()[0], info);
146  if (info > 0) {
147  if (A.m_printLevel) {
148  writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has"
149  " been completed, but the factor U is exactly singular, and division by zero will occur if "
150  "it is used to solve a system of equations.\n", info);
151  }
152  if (!A.m_useReturnErrorCode) {
153  throw CELapackError("solve(DenseMatrix& A, double* b)",
154  "DGETRF returned INFO = "+int2str(info) + ". U(i,i) is exactly zero. The factorization has"
155  " been completed, but the factor U is exactly singular, and division by zero will occur if "
156  "it is used to solve a system of equations.");
157  }
158  return info;
159  } else if (info < 0) {
160  if (A.m_printLevel) {
161  writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
162  }
163 
164  throw CELapackError("solve(DenseMatrix& A, double* b)",
165  "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value");
166  }
167 
168  ct_dgetrs(ctlapack::NoTranspose, A.nRows(), 1, A.ptrColumn(0),
169  A.nRows(), &A.ipiv()[0], b, A.nColumns(), info);
170  if (info != 0) {
171  if (A.m_printLevel) {
172  writelogf("solve(DenseMatrix& A, double* b): DGETRS returned INFO = %d\n", info);
173  }
174  if (info < 0 || !A.m_useReturnErrorCode) {
175  throw CELapackError("solve(DenseMatrix& A, double* b)", "DGETRS returned INFO = "+int2str(info));
176  }
177  }
178  return info;
179 }
180 //====================================================================================================================
181 int solve(DenseMatrix& A, DenseMatrix& b)
182 {
183  int info = 0;
184  if (A.nColumns() != A.nRows()) {
185  if (A.m_printLevel) {
186  writelogf("solve(DenseMatrix& A, DenseMatrix& b): Can only solve a square matrix\n");
187  }
188  if (! A.m_useReturnErrorCode) {
189  throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "Can only solve a square matrix");
190  }
191  return -1;
192  }
193  ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
194  A.nRows(), &A.ipiv()[0], info);
195  if (info != 0) {
196  if (info > 0) {
197  if (A.m_printLevel) {
198  writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has"
199  " been completed, but the factor U is exactly singular, and division by zero will occur if "
200  "it is used to solve a system of equations.\n", info);
201  }
202  if (! A.m_useReturnErrorCode) {
203  throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)",
204  "DGETRF returned INFO = "+int2str(info) + ". U(i,i) is exactly zero. The factorization has"
205  " been completed, but the factor U is exactly singular, and division by zero will occur if "
206  "it is used to solve a system of equations.");
207  }
208  } else {
209  if (A.m_printLevel) {
210  writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
211  }
212  if (! A.m_useReturnErrorCode) {
213  throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)",
214  "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value");
215  }
216  }
217  return info;
218  }
219 
220  ct_dgetrs(ctlapack::NoTranspose, A.nRows(), b.nColumns(), A.ptrColumn(0),
221  A.nRows(), &A.ipiv()[0], b.ptrColumn(0), b.nRows(), info);
222  if (info != 0) {
223  if (A.m_printLevel) {
224  writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRS returned INFO = %d\n", info);
225  }
226  if (! A.m_useReturnErrorCode) {
227  throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "DGETRS returned INFO = "+int2str(info));
228  }
229  }
230 
231  return info;
232 }
233 //====================================================================================================================
234 void multiply(const DenseMatrix& A, const double* const b, double* const prod)
235 {
236  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
237  static_cast<int>(A.nRows()), static_cast<int>(A.nColumns()), 1.0,
238  A.ptrColumn(0), static_cast<int>(A.nRows()), b, 1, 0.0, prod, 1);
239 }
240 //====================================================================================================================
241 void increment(const DenseMatrix& A, const double* b, double* prod)
242 {
243  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
244  static_cast<int>(A.nRows()), static_cast<int>(A.nRows()), 1.0,
245  A.ptrColumn(0), static_cast<int>(A.nRows()), b, 1, 1.0, prod, 1);
246 }
247 //====================================================================================================================
248 int invert(DenseMatrix& A, size_t nn)
249 {
250  integer n = static_cast<int>(nn != npos ? nn : A.nRows());
251  int info=0;
252  ct_dgetrf(n, n, A.ptrColumn(0), static_cast<int>(A.nRows()),
253  &A.ipiv()[0], info);
254  if (info != 0) {
255  if (A.m_printLevel) {
256  writelogf("invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info);
257  }
258  if (! A.m_useReturnErrorCode) {
259  throw CELapackError("invert(DenseMatrix& A, int nn)", "DGETRS returned INFO = "+int2str(info));
260  }
261  return info;
262  }
263 
264  vector_fp work(n);
265  integer lwork = static_cast<int>(work.size());
266  ct_dgetri(n, A.ptrColumn(0), static_cast<int>(A.nRows()),
267  &A.ipiv()[0], &work[0], lwork, info);
268  if (info != 0) {
269  if (A.m_printLevel) {
270  writelogf("invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info);
271  }
272  if (! A.m_useReturnErrorCode) {
273  throw CELapackError("invert(DenseMatrix& A, int nn)", "DGETRI returned INFO="+int2str(info));
274  }
275  }
276  return info;
277 }
278 //====================================================================================================================
279 }
280 
Generated by   doxygen 1.8.2