www.mooseframework.org
ColumnMajorMatrix.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 // MOOSE includes
11 #include "ColumnMajorMatrix.h"
12 
13 #include "DualRealOps.h"
14 
15 #include "libmesh/petsc_macro.h"
16 
17 // PETSc includes
18 #include <petscsys.h>
19 #include <petscblaslapack.h>
20 
21 template <typename T>
22 ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(unsigned int rows, unsigned int cols)
23  : _n_rows(rows), _n_cols(cols), _n_entries(rows * cols), _values(rows * cols, 0.0)
24 {
25  _values.resize(rows * cols);
26 }
27 
28 template <typename T>
30  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
31 {
32  *this = rhs;
33 }
34 
35 template <typename T>
37  : _n_rows(LIBMESH_DIM),
38  _n_cols(LIBMESH_DIM),
39  _n_entries(LIBMESH_DIM * LIBMESH_DIM),
40  _values(LIBMESH_DIM * LIBMESH_DIM)
41 {
42  for (const auto j : make_range(Moose::dim))
43  for (const auto i : make_range(Moose::dim))
44  (*this)(i, j) = rhs(i, j);
45 }
46 
47 template <typename T>
49  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
50 {
51  *this = rhs;
52 }
53 
54 template <typename T>
56  : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
57 {
58  *this = rhs;
59 }
60 
61 template <typename T>
63  const TypeVector<T> & col2,
64  const TypeVector<T> & col3)
65  : _n_rows(LIBMESH_DIM),
66  _n_cols(LIBMESH_DIM),
67  _n_entries(LIBMESH_DIM * LIBMESH_DIM),
68  _values(LIBMESH_DIM * LIBMESH_DIM)
69 {
70  unsigned int entry = 0;
71  for (const auto i : make_range(Moose::dim))
72  _values[entry++] = col1(i);
73 
74  for (const auto i : make_range(Moose::dim))
75  _values[entry++] = col2(i);
76 
77  for (const auto i : make_range(Moose::dim))
78  _values[entry++] = col3(i);
79 }
80 
81 template <typename T>
84 {
85  rhs.checkSquareness();
86 
87  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows * rhs._n_rows, _n_cols * rhs._n_cols);
88 
89  for (unsigned int i = 0; i < _n_rows; i++)
90  for (unsigned int j = 0; j < _n_cols; j++)
91  for (unsigned int k = 0; k < rhs._n_rows; k++)
92  for (unsigned int l = 0; l < rhs._n_cols; l++)
93  ret_matrix(((i * _n_rows) + k), ((j * _n_cols) + l)) = (*this)(i, j) * rhs(k, l);
94 
95  return ret_matrix;
96 }
97 
98 template <typename T>
100 ColumnMajorMatrixTempl<T>::operator=(const DenseMatrix<T> & rhs)
101 {
102  if (_n_rows != rhs.m() || _n_cols != rhs.n())
103  mooseError("ColumnMajorMatrix and DenseMatrix should be of the same shape.");
104 
105  _n_rows = rhs.m();
106  _n_cols = rhs.n();
107  _n_entries = rhs.m() * rhs.n();
108  _values.resize(rhs.m() * rhs.n());
109 
110  for (unsigned int j = 0; j < _n_cols; ++j)
111  for (unsigned int i = 0; i < _n_rows; ++i)
112  (*this)(i, j) = rhs(i, j);
113 
114  return *this;
115 }
116 
117 template <typename T>
119 ColumnMajorMatrixTempl<T>::operator=(const DenseVector<T> & rhs)
120 {
121  if (_n_rows != rhs.size() || _n_cols != 1)
122  mooseError("ColumnMajorMatrix and DenseVector should be of the same shape.");
123 
124  _n_rows = rhs.size();
125  _n_cols = 1;
126  _n_entries = rhs.size();
127  _values.resize(rhs.size());
128 
129  for (unsigned int i = 0; i < _n_rows; ++i)
130  (*this)(i) = rhs(i);
131 
132  return *this;
133 }
134 
135 template <typename T>
136 void
138  ColumnMajorMatrixTempl<T> & evec) const
139 {
140  this->checkSquareness();
141 
142  char jobz = 'V';
143  char uplo = 'U';
144  PetscBLASInt n = _n_rows;
145  PetscBLASInt return_value = 0;
146 
147  eval._n_rows = _n_rows;
148  eval._n_cols = 1;
149  eval._n_entries = _n_rows;
150  eval._values.resize(_n_rows);
151 
152  evec = *this;
153 
154  T * eval_data = eval.rawData();
155  T * evec_data = evec.rawData();
156 
157  PetscBLASInt buffer_size = n * 64;
158  std::vector<T> buffer(buffer_size);
159 
160  LAPACKsyev_(&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
161 
162  if (return_value)
163  mooseError("error in lapack eigen solve");
164 }
165 
166 template <>
167 void
170 {
171  mooseError("Eigen solves with AD types is not supported.");
172 }
173 
174 template <typename T>
175 void
177  ColumnMajorMatrixTempl<T> & eval_img,
178  ColumnMajorMatrixTempl<T> & evec_right,
179  ColumnMajorMatrixTempl<T> & evec_left) const
180 {
181  this->checkSquareness();
182 
183  ColumnMajorMatrixTempl<T> a(*this);
184 
185  char jobvl = 'V';
186  char jobvr = 'V';
187  PetscBLASInt n = _n_rows;
188  PetscBLASInt return_value = 0;
189 
190  eval_real._n_rows = _n_rows;
191  eval_real._n_cols = 1;
192  eval_real._n_entries = _n_rows;
193  eval_real._values.resize(_n_rows);
194 
195  eval_img._n_rows = _n_rows;
196  eval_img._n_cols = 1;
197  eval_img._n_entries = _n_rows;
198  eval_img._values.resize(_n_rows);
199 
200  T * a_data = a.rawData();
201  T * eval_r = eval_real.rawData();
202  T * eval_i = eval_img.rawData();
203  T * evec_ri = evec_right.rawData();
204  T * evec_le = evec_left.rawData();
205 
206  PetscBLASInt buffer_size = n * 64;
207  std::vector<T> buffer(buffer_size);
208 
209  LAPACKgeev_(&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 
224  if (return_value)
225  mooseError("error in lapack eigen solve");
226 }
227 
228 template <>
229 void
234 {
235  mooseError("Eigen solves with AD types is not supported.");
236 }
237 
238 template <typename T>
239 void
241 {
242  this->checkSquareness();
243 
244  ColumnMajorMatrixTempl<T> a(*this);
245  ColumnMajorMatrixTempl<T> evals_real(_n_rows, 1), evals_img(_n_rows, 1),
246  evals_real2(_n_rows, _n_cols);
247  ColumnMajorMatrixTempl<T> evec_right(_n_rows, _n_cols), evec_left(_n_rows, _n_cols);
248  ColumnMajorMatrixTempl<T> evec_right_inverse(_n_rows, _n_cols);
249 
250  a.eigenNonsym(evals_real, evals_img, evec_right, evec_left);
251 
252  for (unsigned int i = 0; i < _n_rows; i++)
253  evals_real2(i, i) = std::exp(evals_real(i, 0));
254 
255  evec_right.inverse(evec_right_inverse);
256 
257  z = evec_right * evals_real2 * evec_right_inverse;
258 }
259 
260 template <typename T>
261 void
263 {
264  this->checkSquareness();
265  this->checkShapeEquality(invA);
266 
267  PetscBLASInt n = _n_rows;
268  PetscBLASInt return_value = 0;
269 
270  invA = *this;
271 
272  std::vector<PetscBLASInt> ipiv(n);
273  T * invA_data = invA.rawData();
274 
275  PetscBLASInt buffer_size = n * 64;
276  std::vector<T> buffer(buffer_size);
277 
278  LAPACKgetrf_(&n, &n, invA_data, &n, &ipiv[0], &return_value);
279 
280  LAPACKgetri_(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
281 
282  if (return_value)
283  mooseException("Error in LAPACK matrix-inverse calculation");
284 }
285 
286 template <typename T>
287 void
289 {
290  if (_n_rows != _n_cols)
291  mooseError("ColumnMajorMatrix error: Unable to perform the operation on a non-square matrix.");
292 }
293 
294 template <typename T>
295 void
297 {
298  if (_n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
299  mooseError("ColumnMajorMatrix error: Unable to perform the operation on matrices of different "
300  "shapes.");
301 }
302 
303 template <>
304 void
306 {
307  mooseError("Inverse solves with AD types is not supported for the ColumnMajorMatrix class.");
308 }
309 
310 template <typename T>
313 {
314  ColumnMajorMatrixTempl<T> & s = (*this);
315 
316  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
317 
318  for (unsigned int j = 0; j < _n_cols; j++)
319  for (unsigned int i = 0; i < _n_rows; i++)
320  ret_matrix(i, j) = std::abs(s(i, j));
321 
322  return ret_matrix;
323 }
324 
325 template <typename T>
326 inline T
328 {
329  return std::sqrt(doubleContraction(*this));
330 }
331 
332 template class ColumnMajorMatrixTempl<Real>;
333 template class ColumnMajorMatrixTempl<DualReal>;
void inverse(ColumnMajorMatrixTempl< T > &invA) const
Returns inverse of a general matrix.
This class defines a Tensor that can change its shape.
auto exp(const T &)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
ColumnMajorMatrixTempl< T > abs()
Returns a matrix that is the absolute value of the matrix this was called on.
T norm()
The Euclidean norm of the matrix.
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:148
ColumnMajorMatrixTempl< T > kronecker(const ColumnMajorMatrixTempl< T > &rhs) const
Kronecker Product.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
void eigen(ColumnMajorMatrixTempl< T > &eval, ColumnMajorMatrixTempl< T > &evec) const
Returns eigen system solve for a symmetric real matrix.
void checkSquareness() const
Check if matrix is square.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
T * rawData()
Returns a reference to the raw data pointer.
void eigenNonsym(ColumnMajorMatrixTempl< T > &eval_real, ColumnMajorMatrixTempl< T > &eval_img, ColumnMajorMatrixTempl< T > &evec_right, ColumnMajorMatrixTempl< T > &eve_left) const
Returns eigen system solve for a non-symmetric real matrix.
void checkShapeEquality(const ColumnMajorMatrixTempl< T > &rhs) const
Check if matrices are of same shape.
ColumnMajorMatrixTempl(const unsigned int rows=Moose::dim, const unsigned int cols=Moose::dim)
Constructor that sets an initial number of entries and shape.
std::vector< T > _values
IntRange< T > make_range(T beg, T end)
ColumnMajorMatrixTempl< T > & operator=(const TypeTensor< T > &rhs)
Sets the values in this tensor to the values on the RHS.
void exp(ColumnMajorMatrixTempl< T > &z) const
Returns matrix that is the exponential of the matrix this was called on.