Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 "libmesh/petsc_macro.h"
14 :
15 : // PETSc includes
16 : #include <petscsys.h>
17 : #include <petscblaslapack.h>
18 :
19 : template <typename T>
20 581 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(unsigned int rows, unsigned int cols)
21 581 : : _n_rows(rows), _n_cols(cols), _n_entries(rows * cols), _values(rows * cols, 0.0)
22 : {
23 581 : _values.resize(rows * cols);
24 581 : }
25 :
26 : template <typename T>
27 175 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const ColumnMajorMatrixTempl<T> & rhs)
28 175 : : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
29 : {
30 175 : *this = rhs;
31 175 : }
32 :
33 : template <typename T>
34 1 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const TypeTensor<T> & rhs)
35 1 : : _n_rows(LIBMESH_DIM),
36 1 : _n_cols(LIBMESH_DIM),
37 1 : _n_entries(LIBMESH_DIM * LIBMESH_DIM),
38 1 : _values(LIBMESH_DIM * LIBMESH_DIM)
39 : {
40 4 : for (const auto j : make_range(Moose::dim))
41 12 : for (const auto i : make_range(Moose::dim))
42 9 : (*this)(i, j) = rhs(i, j);
43 1 : }
44 :
45 : template <typename T>
46 0 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs)
47 0 : : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
48 : {
49 0 : *this = rhs;
50 0 : }
51 :
52 : template <typename T>
53 0 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const DenseVector<T> & rhs)
54 0 : : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
55 : {
56 0 : *this = rhs;
57 0 : }
58 :
59 : template <typename T>
60 1 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const TypeVector<T> & col1,
61 : const TypeVector<T> & col2,
62 : const TypeVector<T> & col3)
63 1 : : _n_rows(LIBMESH_DIM),
64 1 : _n_cols(LIBMESH_DIM),
65 1 : _n_entries(LIBMESH_DIM * LIBMESH_DIM),
66 1 : _values(LIBMESH_DIM * LIBMESH_DIM)
67 : {
68 1 : unsigned int entry = 0;
69 4 : for (const auto i : make_range(Moose::dim))
70 3 : _values[entry++] = col1(i);
71 :
72 4 : for (const auto i : make_range(Moose::dim))
73 3 : _values[entry++] = col2(i);
74 :
75 4 : for (const auto i : make_range(Moose::dim))
76 3 : _values[entry++] = col3(i);
77 1 : }
78 :
79 : template <typename T>
80 : ColumnMajorMatrixTempl<T>
81 2 : ColumnMajorMatrixTempl<T>::kronecker(const ColumnMajorMatrixTempl<T> & rhs) const
82 : {
83 2 : rhs.checkSquareness();
84 :
85 1 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows * rhs._n_rows, _n_cols * rhs._n_cols);
86 :
87 3 : for (unsigned int i = 0; i < _n_rows; i++)
88 6 : for (unsigned int j = 0; j < _n_cols; j++)
89 12 : for (unsigned int k = 0; k < rhs._n_rows; k++)
90 24 : for (unsigned int l = 0; l < rhs._n_cols; l++)
91 16 : ret_matrix(((i * _n_rows) + k), ((j * _n_cols) + l)) = (*this)(i, j) * rhs(k, l);
92 :
93 1 : return ret_matrix;
94 0 : }
95 :
96 : template <typename T>
97 : ColumnMajorMatrixTempl<T> &
98 1 : ColumnMajorMatrixTempl<T>::operator=(const DenseMatrix<T> & rhs)
99 : {
100 1 : if (_n_rows != rhs.m() || _n_cols != rhs.n())
101 1 : mooseError("ColumnMajorMatrix and DenseMatrix should be of the same shape.");
102 :
103 0 : _n_rows = rhs.m();
104 0 : _n_cols = rhs.n();
105 0 : _n_entries = rhs.m() * rhs.n();
106 0 : _values.resize(rhs.m() * rhs.n());
107 :
108 0 : for (unsigned int j = 0; j < _n_cols; ++j)
109 0 : for (unsigned int i = 0; i < _n_rows; ++i)
110 0 : (*this)(i, j) = rhs(i, j);
111 :
112 0 : return *this;
113 : }
114 :
115 : template <typename T>
116 : ColumnMajorMatrixTempl<T> &
117 1 : ColumnMajorMatrixTempl<T>::operator=(const DenseVector<T> & rhs)
118 : {
119 1 : if (_n_rows != rhs.size() || _n_cols != 1)
120 1 : mooseError("ColumnMajorMatrix and DenseVector should be of the same shape.");
121 :
122 0 : _n_rows = rhs.size();
123 0 : _n_cols = 1;
124 0 : _n_entries = rhs.size();
125 0 : _values.resize(rhs.size());
126 :
127 0 : for (unsigned int i = 0; i < _n_rows; ++i)
128 0 : (*this)(i) = rhs(i);
129 :
130 0 : return *this;
131 : }
132 :
133 : template <typename T>
134 : void
135 1 : ColumnMajorMatrixTempl<T>::eigen(ColumnMajorMatrixTempl<T> & eval,
136 : ColumnMajorMatrixTempl<T> & evec) const
137 : {
138 1 : this->checkSquareness();
139 :
140 1 : char jobz = 'V';
141 1 : char uplo = 'U';
142 1 : PetscBLASInt n = _n_rows;
143 1 : PetscBLASInt return_value = 0;
144 :
145 1 : eval._n_rows = _n_rows;
146 1 : eval._n_cols = 1;
147 1 : eval._n_entries = _n_rows;
148 1 : eval._values.resize(_n_rows);
149 :
150 1 : evec = *this;
151 :
152 1 : T * eval_data = eval.rawData();
153 1 : T * evec_data = evec.rawData();
154 :
155 1 : PetscBLASInt buffer_size = n * 64;
156 1 : std::vector<T> buffer(buffer_size);
157 :
158 1 : LAPACKsyev_(&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
159 :
160 1 : if (return_value)
161 0 : mooseError("error in lapack eigen solve");
162 1 : }
163 :
164 : template <>
165 : void
166 0 : ColumnMajorMatrixTempl<ADReal>::eigen(ColumnMajorMatrixTempl<ADReal> &,
167 : ColumnMajorMatrixTempl<ADReal> &) const
168 : {
169 0 : mooseError("Eigen solves with AD types is not supported.");
170 : }
171 :
172 : template <typename T>
173 : void
174 3 : ColumnMajorMatrixTempl<T>::eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
175 : ColumnMajorMatrixTempl<T> & eval_img,
176 : ColumnMajorMatrixTempl<T> & evec_right,
177 : ColumnMajorMatrixTempl<T> & evec_left) const
178 : {
179 3 : this->checkSquareness();
180 :
181 3 : ColumnMajorMatrixTempl<T> a(*this);
182 :
183 3 : char jobvl = 'V';
184 3 : char jobvr = 'V';
185 3 : PetscBLASInt n = _n_rows;
186 3 : PetscBLASInt return_value = 0;
187 :
188 3 : eval_real._n_rows = _n_rows;
189 3 : eval_real._n_cols = 1;
190 3 : eval_real._n_entries = _n_rows;
191 3 : eval_real._values.resize(_n_rows);
192 :
193 3 : eval_img._n_rows = _n_rows;
194 3 : eval_img._n_cols = 1;
195 3 : eval_img._n_entries = _n_rows;
196 3 : eval_img._values.resize(_n_rows);
197 :
198 3 : T * a_data = a.rawData();
199 3 : T * eval_r = eval_real.rawData();
200 3 : T * eval_i = eval_img.rawData();
201 3 : T * evec_ri = evec_right.rawData();
202 3 : T * evec_le = evec_left.rawData();
203 :
204 3 : PetscBLASInt buffer_size = n * 64;
205 3 : std::vector<T> buffer(buffer_size);
206 :
207 3 : LAPACKgeev_(&jobvl,
208 : &jobvr,
209 : &n,
210 : a_data,
211 : &n,
212 : eval_r,
213 : eval_i,
214 : evec_le,
215 : &n,
216 : evec_ri,
217 : &n,
218 3 : &buffer[0],
219 : &buffer_size,
220 : &return_value);
221 :
222 3 : if (return_value)
223 0 : mooseError("error in lapack eigen solve");
224 3 : }
225 :
226 : template <>
227 : void
228 0 : ColumnMajorMatrixTempl<ADReal>::eigenNonsym(ColumnMajorMatrixTempl<ADReal> &,
229 : ColumnMajorMatrixTempl<ADReal> &,
230 : ColumnMajorMatrixTempl<ADReal> &,
231 : ColumnMajorMatrixTempl<ADReal> &) const
232 : {
233 0 : mooseError("Eigen solves with AD types is not supported.");
234 : }
235 :
236 : template <typename T>
237 : void
238 2 : ColumnMajorMatrixTempl<T>::exp(ColumnMajorMatrixTempl<T> & z) const
239 : {
240 2 : this->checkSquareness();
241 :
242 2 : ColumnMajorMatrixTempl<T> a(*this);
243 2 : ColumnMajorMatrixTempl<T> evals_real(_n_rows, 1), evals_img(_n_rows, 1),
244 2 : evals_real2(_n_rows, _n_cols);
245 2 : ColumnMajorMatrixTempl<T> evec_right(_n_rows, _n_cols), evec_left(_n_rows, _n_cols);
246 2 : ColumnMajorMatrixTempl<T> evec_right_inverse(_n_rows, _n_cols);
247 :
248 2 : a.eigenNonsym(evals_real, evals_img, evec_right, evec_left);
249 :
250 8 : for (unsigned int i = 0; i < _n_rows; i++)
251 6 : evals_real2(i, i) = std::exp(evals_real(i, 0));
252 :
253 2 : evec_right.inverse(evec_right_inverse);
254 :
255 2 : z = evec_right * evals_real2 * evec_right_inverse;
256 2 : }
257 :
258 : template <typename T>
259 : void
260 3 : ColumnMajorMatrixTempl<T>::inverse(ColumnMajorMatrixTempl<T> & invA) const
261 : {
262 3 : this->checkSquareness();
263 3 : this->checkShapeEquality(invA);
264 :
265 3 : PetscBLASInt n = _n_rows;
266 3 : PetscBLASInt return_value = 0;
267 :
268 3 : invA = *this;
269 :
270 3 : std::vector<PetscBLASInt> ipiv(n);
271 3 : T * invA_data = invA.rawData();
272 :
273 3 : PetscBLASInt buffer_size = n * 64;
274 3 : std::vector<T> buffer(buffer_size);
275 :
276 3 : LAPACKgetrf_(&n, &n, invA_data, &n, &ipiv[0], &return_value);
277 :
278 3 : LAPACKgetri_(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
279 :
280 3 : if (return_value)
281 0 : mooseException("Error in LAPACK matrix-inverse calculation");
282 3 : }
283 :
284 : template <typename T>
285 : void
286 16 : ColumnMajorMatrixTempl<T>::checkSquareness() const
287 : {
288 16 : if (_n_rows != _n_cols)
289 3 : mooseError("ColumnMajorMatrix error: Unable to perform the operation on a non-square matrix.");
290 13 : }
291 :
292 : template <typename T>
293 : void
294 10 : ColumnMajorMatrixTempl<T>::checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const
295 : {
296 10 : if (_n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
297 2 : mooseError("ColumnMajorMatrix error: Unable to perform the operation on matrices of different "
298 : "shapes.");
299 8 : }
300 :
301 : template <>
302 : void
303 0 : ColumnMajorMatrixTempl<ADReal>::inverse(ColumnMajorMatrixTempl<ADReal> &) const
304 : {
305 0 : mooseError("Inverse solves with AD types is not supported for the ColumnMajorMatrix class.");
306 : }
307 :
308 : template <typename T>
309 : inline ColumnMajorMatrixTempl<T>
310 0 : ColumnMajorMatrixTempl<T>::abs()
311 : {
312 0 : ColumnMajorMatrixTempl<T> & s = (*this);
313 :
314 0 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
315 :
316 0 : for (unsigned int j = 0; j < _n_cols; j++)
317 0 : for (unsigned int i = 0; i < _n_rows; i++)
318 0 : ret_matrix(i, j) = std::abs(s(i, j));
319 :
320 0 : return ret_matrix;
321 0 : }
322 :
323 : template <typename T>
324 : inline T
325 1 : ColumnMajorMatrixTempl<T>::norm()
326 : {
327 1 : return std::sqrt(doubleContraction(*this));
328 : }
329 :
330 : template class ColumnMajorMatrixTempl<Real>;
331 : template class ColumnMajorMatrixTempl<ADReal>;
|
