www.mooseframework.org
ColumnMajorMatrix.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 // MOOSE includes
16 #include "ColumnMajorMatrix.h"
17 
18 #include "libmesh/petsc_macro.h"
19 
20 // PETSc includes
21 #include <petscsys.h>
22 #include <petscblaslapack.h>
23 
24 #if !defined(LIBMESH_HAVE_PETSC)
25 extern "C" void FORTRAN_CALL(dsyev)(...);
26 extern "C" void FORTRAN_CALL(dgeev)(...);
27 extern "C" void FORTRAN_CALL(dgetrf)(...);
28 #endif
29 #if !defined(LIBMESH_HAVE_PETSC) || PETSC_VERSION_LESS_THAN(3, 5, 0)
30 extern "C" void FORTRAN_CALL(dgetri)(...); // matrix inversion routine from LAPACK
31 #endif
32 
33 ColumnMajorMatrix::ColumnMajorMatrix(unsigned int rows, unsigned int cols)
34  : _n_rows(rows), _n_cols(cols), _n_entries(rows * cols), _values(rows * cols, 0.0)
35 {
36  _values.resize(rows * cols);
37 }
38 
40  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
41 {
42  *this = rhs;
43 }
44 
45 ColumnMajorMatrix::ColumnMajorMatrix(const TypeTensor<Real> & rhs)
46  : _n_rows(LIBMESH_DIM),
47  _n_cols(LIBMESH_DIM),
48  _n_entries(LIBMESH_DIM * LIBMESH_DIM),
49  _values(LIBMESH_DIM * LIBMESH_DIM)
50 {
51  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
52  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
53  (*this)(i, j) = rhs(i, j);
54 }
55 
56 ColumnMajorMatrix::ColumnMajorMatrix(const DenseMatrix<Real> & rhs)
57  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
58 {
59  *this = rhs;
60 }
61 
62 ColumnMajorMatrix::ColumnMajorMatrix(const DenseVector<Real> & rhs)
63  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
64 {
65  *this = rhs;
66 }
67 
68 ColumnMajorMatrix::ColumnMajorMatrix(const TypeVector<Real> & col1,
69  const TypeVector<Real> & col2,
70  const TypeVector<Real> & col3)
71  : _n_rows(LIBMESH_DIM),
72  _n_cols(LIBMESH_DIM),
73  _n_entries(LIBMESH_DIM * LIBMESH_DIM),
74  _values(LIBMESH_DIM * LIBMESH_DIM)
75 {
76  unsigned int entry = 0;
77  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
78  _values[entry++] = col1(i);
79 
80  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
81  _values[entry++] = col2(i);
82 
83  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
84  _values[entry++] = col3(i);
85 }
86 
89 {
90  mooseAssert(_n_rows == rhs._n_cols && _n_cols == rhs._n_rows,
91  "Matrices must be the same shape for a kronecker product!");
92 
93  ColumnMajorMatrix ret_matrix(_n_rows * _n_rows, rhs._n_cols * rhs._n_cols);
94 
95  for (unsigned int i = 0; i < _n_rows; i++)
96  for (unsigned int j = 0; j < _n_cols; j++)
97  for (unsigned int k = 0; k < rhs._n_rows; k++)
98  for (unsigned int l = 0; l < rhs._n_cols; l++)
99  ret_matrix(((i * _n_rows) + k), ((j * _n_cols) + l)) = (*this)(i, j) * rhs(k, l);
100 
101  return ret_matrix;
102 }
103 
105 ColumnMajorMatrix::operator=(const DenseMatrix<Real> & rhs)
106 {
107  mooseAssert(_n_rows == rhs.m(), "different number of rows");
108  mooseAssert(_n_cols == rhs.n(), "different number of cols");
109 
110  _n_rows = rhs.m();
111  _n_cols = rhs.n();
112  _n_entries = rhs.m() * rhs.n();
113  _values.resize(rhs.m() * rhs.n());
114 
115  for (unsigned int j = 0; j < _n_cols; ++j)
116  for (unsigned int i = 0; i < _n_rows; ++i)
117  (*this)(i, j) = rhs(i, j);
118 
119  return *this;
120 }
121 
123 ColumnMajorMatrix::operator=(const DenseVector<Real> & rhs)
124 {
125  mooseAssert(_n_rows == rhs.size(), "different number of rows");
126  mooseAssert(_n_cols == 1, "different number of cols");
127 
128  _n_rows = rhs.size();
129  _n_cols = 1;
130  _n_entries = rhs.size();
131  _values.resize(rhs.size());
132 
133  for (unsigned int i = 0; i < _n_rows; ++i)
134  (*this)(i) = rhs(i);
135 
136  return *this;
137 }
138 
139 void
141 {
142  mooseAssert(_n_rows == _n_cols, "Cannot solve eigen system of a non-square matrix!");
143 
144  char jobz = 'V';
145  char uplo = 'U';
146  int n = _n_rows;
147  int return_value = 0;
148 
149  eval._n_rows = _n_rows;
150  eval._n_cols = 1;
151  eval._n_entries = _n_rows;
152  eval._values.resize(_n_rows);
153 
154  evec = *this;
155 
156  Real * eval_data = eval.rawData();
157  Real * evec_data = evec.rawData();
158 
159  int buffer_size = n * 64;
160  std::vector<Real> buffer(buffer_size);
161 
162 #if !defined(LIBMESH_HAVE_PETSC)
163  FORTRAN_CALL(dsyev)
164  (&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
165 #else
166  LAPACKsyev_(&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
167 #endif
168 
169  if (return_value)
170  mooseError("error in lapack eigen solve");
171 }
172 
173 void
175  ColumnMajorMatrix & eval_img,
176  ColumnMajorMatrix & evec_right,
177  ColumnMajorMatrix & evec_left) const
178 {
179  mooseAssert(_n_rows == _n_cols, "Cannot solve eigen system of a non-square matrix!");
180 
181  ColumnMajorMatrix a(*this);
182 
183  char jobvl = 'V';
184  char jobvr = 'V';
185  int n = _n_rows;
186  int return_value = 0;
187 
188  eval_real._n_rows = _n_rows;
189  eval_real._n_cols = 1;
190  eval_real._n_entries = _n_rows;
191  eval_real._values.resize(_n_rows);
192 
193  eval_img._n_rows = _n_rows;
194  eval_img._n_cols = 1;
195  eval_img._n_entries = _n_rows;
196  eval_img._values.resize(_n_rows);
197 
198  Real * a_data = a.rawData();
199  Real * eval_r = eval_real.rawData();
200  Real * eval_i = eval_img.rawData();
201  Real * evec_ri = evec_right.rawData();
202  Real * evec_le = evec_left.rawData();
203 
204  int buffer_size = n * 64;
205  std::vector<Real> buffer(buffer_size);
206 
207 #if !defined(LIBMESH_HAVE_PETSC)
208  FORTRAN_CALL(dgeev)
209  (&jobvl,
210  &jobvr,
211  &n,
212  a_data,
213  &n,
214  eval_r,
215  eval_i,
216  evec_le,
217  &n,
218  evec_ri,
219  &n,
220  &buffer[0],
221  &buffer_size,
222  &return_value);
223 #else
224  LAPACKgeev_(&jobvl,
225  &jobvr,
226  &n,
227  a_data,
228  &n,
229  eval_r,
230  eval_i,
231  evec_le,
232  &n,
233  evec_ri,
234  &n,
235  &buffer[0],
236  &buffer_size,
237  &return_value);
238 #endif
239 
240  if (return_value)
241  mooseError("error in lapack eigen solve");
242 }
243 
244 void
246 {
247  mooseAssert(_n_rows == _n_cols, "The Matrix being exponentiated is not square");
248  ColumnMajorMatrix a(*this);
249  ColumnMajorMatrix evals_real(_n_rows, 1), evals_img(_n_rows, 1), evals_real2(_n_rows, _n_cols);
250  ColumnMajorMatrix evec_right(_n_rows, _n_cols), evec_left(_n_rows, _n_cols);
251  ColumnMajorMatrix evec_right_inverse(_n_rows, _n_cols);
252 
253  a.eigenNonsym(evals_real, evals_img, evec_right, evec_left);
254 
255  for (unsigned int i = 0; i < _n_rows; i++)
256  evals_real2(i, i) = std::exp(evals_real(i, 0));
257 
258  evec_right.inverse(evec_right_inverse);
259 
260  z = evec_right * evals_real2 * evec_right_inverse;
261 }
262 
263 void
265 {
266  mooseAssert(_n_rows == _n_cols, "Cannot solve for inverse of a non-square matrix!");
267  mooseAssert(_n_rows == invA._n_cols && _n_cols == invA._n_rows,
268  "Matrices must be the same size for matrix inverse!");
269 
270  int n = _n_rows;
271  int return_value = 0;
272 
273  invA = *this;
274 
275  std::vector<PetscBLASInt> ipiv(n);
276  Real * invA_data = invA.rawData();
277 
278  int buffer_size = n * 64;
279  std::vector<Real> buffer(buffer_size);
280 
281 #if !defined(LIBMESH_HAVE_PETSC)
282  FORTRAN_CALL(dgetrf)(&n, &n, invA_data, &n, &ipiv[0], &return_value);
283 #else
284  LAPACKgetrf_(&n, &n, invA_data, &n, &ipiv[0], &return_value);
285 #endif
286 
287 #if !defined(LIBMESH_HAVE_PETSC) || PETSC_VERSION_LESS_THAN(3, 5, 0)
288  FORTRAN_CALL(dgetri)(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
289 #else
290  LAPACKgetri_(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
291 #endif
292 
293  if (return_value)
294  mooseException("Error in LAPACK matrix-inverse calculation");
295 }
ColumnMajorMatrix kronecker(const ColumnMajorMatrix &rhs) const
Kronecker Product.
unsigned int n() const
Returns the number of rows.
Real * rawData()
Returns a reference to the raw data pointer.
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
unsigned int _n_entries
This class defines a Tensor that can change its shape.
ColumnMajorMatrix(const unsigned int rows=LIBMESH_DIM, const unsigned int cols=LIBMESH_DIM)
Constructor that sets an initial number of entries and shape.
std::vector< Real > _values
ColumnMajorMatrix & operator=(const TypeTensor< Real > &rhs)
Sets the values in this tensor to the values on the RHS.
void FORTRAN_CALL() dgetri(...)
void FORTRAN_CALL() dgetrf(...)
void eigenNonsym(ColumnMajorMatrix &eval_real, ColumnMajorMatrix &eval_img, ColumnMajorMatrix &evec_right, ColumnMajorMatrix &eve_left) const
Returns eigen system solve for a non-symmetric real matrix.
void FORTRAN_CALL() dsyev(...)
void exp(ColumnMajorMatrix &z) const
Returns matrix that is the exponential of the matrix this was called on.
void FORTRAN_CALL() dgeev(...)
void inverse(ColumnMajorMatrix &invA) const
Returns inverse of a general matrix.
void eigen(ColumnMajorMatrix &eval, ColumnMajorMatrix &evec) const
Returns eigen system solve for a symmetric real matrix.