15 #include "libmesh/petsc_macro.h" 19 #include <petscblaslapack.h> 23 : _n_rows(rows), _n_cols(cols), _n_entries(rows * cols), _values(rows * cols, 0.0)
30 : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
37 : _n_rows(LIBMESH_DIM),
39 _n_entries(LIBMESH_DIM * LIBMESH_DIM),
40 _values(LIBMESH_DIM * LIBMESH_DIM)
44 (*this)(i, j) = rhs(i, j);
49 : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
56 : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
63 const TypeVector<T> & col2,
64 const TypeVector<T> & col3)
65 : _n_rows(LIBMESH_DIM),
67 _n_entries(LIBMESH_DIM * LIBMESH_DIM),
68 _values(LIBMESH_DIM * LIBMESH_DIM)
70 unsigned int entry = 0;
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);
102 if (_n_rows != rhs.m() || _n_cols != rhs.n())
103 mooseError(
"ColumnMajorMatrix and DenseMatrix should be of the same shape.");
107 _n_entries = rhs.m() * rhs.n();
108 _values.resize(rhs.m() * rhs.n());
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);
117 template <
typename T>
121 if (_n_rows != rhs.size() || _n_cols != 1)
122 mooseError(
"ColumnMajorMatrix and DenseVector should be of the same shape.");
124 _n_rows = rhs.size();
126 _n_entries = rhs.size();
127 _values.resize(rhs.size());
129 for (
unsigned int i = 0; i < _n_rows; ++i)
135 template <
typename T>
140 this->checkSquareness();
144 PetscBLASInt n = _n_rows;
145 PetscBLASInt return_value = 0;
154 T * eval_data = eval.
rawData();
155 T * evec_data = evec.
rawData();
157 PetscBLASInt buffer_size = n * 64;
158 std::vector<T> buffer(buffer_size);
160 LAPACKsyev_(&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
171 mooseError(
"Eigen solves with AD types is not supported.");
174 template <
typename T>
181 this->checkSquareness();
187 PetscBLASInt n = _n_rows;
188 PetscBLASInt return_value = 0;
193 eval_real.
_values.resize(_n_rows);
198 eval_img.
_values.resize(_n_rows);
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();
206 PetscBLASInt buffer_size = n * 64;
207 std::vector<T> buffer(buffer_size);
235 mooseError(
"Eigen solves with AD types is not supported.");
238 template <
typename T>
242 this->checkSquareness();
246 evals_real2(_n_rows, _n_cols);
250 a.
eigenNonsym(evals_real, evals_img, evec_right, evec_left);
252 for (
unsigned int i = 0; i < _n_rows; i++)
253 evals_real2(i, i) =
std::exp(evals_real(i, 0));
255 evec_right.inverse(evec_right_inverse);
257 z = evec_right * evals_real2 * evec_right_inverse;
260 template <
typename T>
264 this->checkSquareness();
265 this->checkShapeEquality(invA);
267 PetscBLASInt n = _n_rows;
268 PetscBLASInt return_value = 0;
272 std::vector<PetscBLASInt> ipiv(n);
273 T * invA_data = invA.
rawData();
275 PetscBLASInt buffer_size = n * 64;
276 std::vector<T> buffer(buffer_size);
278 LAPACKgetrf_(&n, &n, invA_data, &n, &ipiv[0], &return_value);
280 LAPACKgetri_(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
283 mooseException(
"Error in LAPACK matrix-inverse calculation");
286 template <
typename T>
290 if (_n_rows != _n_cols)
291 mooseError(
"ColumnMajorMatrix error: Unable to perform the operation on a non-square matrix.");
294 template <
typename T>
299 mooseError(
"ColumnMajorMatrix error: Unable to perform the operation on matrices of different " 307 mooseError(
"Inverse solves with AD types is not supported for the ColumnMajorMatrix class.");
310 template <
typename T>
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));
325 template <
typename T>
329 return std::sqrt(doubleContraction(*
this));
void inverse(ColumnMajorMatrixTempl< T > &invA) const
Returns inverse of a general matrix.
This class defines a Tensor that can change its shape.
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
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.
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.
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.