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 1006 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(unsigned int rows, unsigned int cols)
21 2012 : : _n_rows(rows), _n_cols(cols), _n_entries(rows * cols), _values(rows * cols, 0.0)
22 : {
23 1006 : _values.resize(rows * cols);
24 1006 : }
25 :
26 : template <typename T>
27 206 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const ColumnMajorMatrixTempl<T> & rhs)
28 206 : : _n_rows(LIBMESH_DIM), _n_cols(LIBMESH_DIM), _n_entries(_n_cols * _n_cols)
29 : {
30 206 : *this = rhs;
31 206 : }
32 :
33 : template <typename T>
34 2 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const TypeTensor<T> & rhs)
35 2 : : _n_rows(LIBMESH_DIM),
36 2 : _n_cols(LIBMESH_DIM),
37 2 : _n_entries(LIBMESH_DIM * LIBMESH_DIM),
38 4 : _values(LIBMESH_DIM * LIBMESH_DIM)
39 : {
40 8 : for (const auto j : make_range(Moose::dim))
41 24 : for (const auto i : make_range(Moose::dim))
42 18 : (*this)(i, j) = rhs(i, j);
43 2 : }
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 2 : ColumnMajorMatrixTempl<T>::ColumnMajorMatrixTempl(const TypeVector<T> & col1,
61 : const TypeVector<T> & col2,
62 : const TypeVector<T> & col3)
63 2 : : _n_rows(LIBMESH_DIM),
64 2 : _n_cols(LIBMESH_DIM),
65 2 : _n_entries(LIBMESH_DIM * LIBMESH_DIM),
66 4 : _values(LIBMESH_DIM * LIBMESH_DIM)
67 : {
68 2 : unsigned int entry = 0;
69 8 : for (const auto i : make_range(Moose::dim))
70 6 : _values[entry++] = col1(i);
71 :
72 8 : for (const auto i : make_range(Moose::dim))
73 6 : _values[entry++] = col2(i);
74 :
75 8 : for (const auto i : make_range(Moose::dim))
76 6 : _values[entry++] = col3(i);
77 2 : }
78 :
79 : template <typename T>
80 : ColumnMajorMatrixTempl<T>
81 4 : ColumnMajorMatrixTempl<T>::kronecker(const ColumnMajorMatrixTempl<T> & rhs) const
82 : {
83 4 : rhs.checkSquareness();
84 :
85 2 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows * rhs._n_rows, _n_cols * rhs._n_cols);
86 :
87 6 : for (unsigned int i = 0; i < _n_rows; i++)
88 12 : for (unsigned int j = 0; j < _n_cols; j++)
89 24 : for (unsigned int k = 0; k < rhs._n_rows; k++)
90 48 : for (unsigned int l = 0; l < rhs._n_cols; l++)
91 32 : ret_matrix(((i * _n_rows) + k), ((j * _n_cols) + l)) = (*this)(i, j) * rhs(k, l);
92 :
93 2 : return ret_matrix;
94 0 : }
95 :
96 : template <typename T>
97 : ColumnMajorMatrixTempl<T> &
98 2 : ColumnMajorMatrixTempl<T>::operator=(const DenseMatrix<T> & rhs)
99 : {
100 2 : if (_n_rows != rhs.m() || _n_cols != rhs.n())
101 2 : 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 2 : ColumnMajorMatrixTempl<T>::operator=(const DenseVector<T> & rhs)
118 : {
119 2 : if (_n_rows != rhs.size() || _n_cols != 1)
120 2 : 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 2 : ColumnMajorMatrixTempl<T>::eigen(ColumnMajorMatrixTempl<T> & eval,
136 : ColumnMajorMatrixTempl<T> & evec) const
137 : {
138 2 : this->checkSquareness();
139 :
140 2 : char jobz = 'V';
141 2 : char uplo = 'U';
142 2 : PetscBLASInt n = _n_rows;
143 2 : PetscBLASInt return_value = 0;
144 :
145 2 : eval._n_rows = _n_rows;
146 2 : eval._n_cols = 1;
147 2 : eval._n_entries = _n_rows;
148 2 : eval._values.resize(_n_rows);
149 :
150 2 : evec = *this;
151 :
152 2 : T * eval_data = eval.rawData();
153 2 : T * evec_data = evec.rawData();
154 :
155 2 : PetscBLASInt buffer_size = n * 64;
156 2 : std::vector<T> buffer(buffer_size);
157 :
158 2 : LAPACKsyev_(&jobz, &uplo, &n, evec_data, &n, eval_data, &buffer[0], &buffer_size, &return_value);
159 :
160 2 : if (return_value)
161 0 : mooseError("error in lapack eigen solve");
162 2 : }
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 6 : 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 6 : this->checkSquareness();
180 :
181 6 : ColumnMajorMatrixTempl<T> a(*this);
182 :
183 6 : char jobvl = 'V';
184 6 : char jobvr = 'V';
185 6 : PetscBLASInt n = _n_rows;
186 6 : PetscBLASInt return_value = 0;
187 :
188 6 : eval_real._n_rows = _n_rows;
189 6 : eval_real._n_cols = 1;
190 6 : eval_real._n_entries = _n_rows;
191 6 : eval_real._values.resize(_n_rows);
192 :
193 6 : eval_img._n_rows = _n_rows;
194 6 : eval_img._n_cols = 1;
195 6 : eval_img._n_entries = _n_rows;
196 6 : eval_img._values.resize(_n_rows);
197 :
198 6 : T * a_data = a.rawData();
199 6 : T * eval_r = eval_real.rawData();
200 6 : T * eval_i = eval_img.rawData();
201 6 : T * evec_ri = evec_right.rawData();
202 6 : T * evec_le = evec_left.rawData();
203 :
204 6 : PetscBLASInt buffer_size = n * 64;
205 6 : std::vector<T> buffer(buffer_size);
206 :
207 6 : 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 6 : &buffer[0],
219 : &buffer_size,
220 : &return_value);
221 :
222 6 : if (return_value)
223 0 : mooseError("error in lapack eigen solve");
224 6 : }
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 4 : ColumnMajorMatrixTempl<T>::exp(ColumnMajorMatrixTempl<T> & z) const
239 : {
240 : using std::exp;
241 4 : this->checkSquareness();
242 :
243 4 : ColumnMajorMatrixTempl<T> a(*this);
244 4 : ColumnMajorMatrixTempl<T> evals_real(_n_rows, 1), evals_img(_n_rows, 1),
245 4 : evals_real2(_n_rows, _n_cols);
246 4 : ColumnMajorMatrixTempl<T> evec_right(_n_rows, _n_cols), evec_left(_n_rows, _n_cols);
247 4 : ColumnMajorMatrixTempl<T> evec_right_inverse(_n_rows, _n_cols);
248 :
249 4 : a.eigenNonsym(evals_real, evals_img, evec_right, evec_left);
250 :
251 16 : for (unsigned int i = 0; i < _n_rows; i++)
252 12 : evals_real2(i, i) = exp(evals_real(i, 0));
253 :
254 4 : evec_right.inverse(evec_right_inverse);
255 :
256 4 : z = evec_right * evals_real2 * evec_right_inverse;
257 4 : }
258 :
259 : template <typename T>
260 : void
261 6 : ColumnMajorMatrixTempl<T>::inverse(ColumnMajorMatrixTempl<T> & invA) const
262 : {
263 6 : this->checkSquareness();
264 6 : this->checkShapeEquality(invA);
265 :
266 6 : PetscBLASInt n = _n_rows;
267 6 : PetscBLASInt return_value = 0;
268 :
269 6 : invA = *this;
270 :
271 6 : std::vector<PetscBLASInt> ipiv(n);
272 6 : T * invA_data = invA.rawData();
273 :
274 6 : PetscBLASInt buffer_size = n * 64;
275 6 : std::vector<T> buffer(buffer_size);
276 :
277 6 : LAPACKgetrf_(&n, &n, invA_data, &n, &ipiv[0], &return_value);
278 :
279 6 : LAPACKgetri_(&n, invA_data, &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
280 :
281 6 : if (return_value)
282 0 : mooseException("Error in LAPACK matrix-inverse calculation");
283 6 : }
284 :
285 : template <typename T>
286 : void
287 32 : ColumnMajorMatrixTempl<T>::checkSquareness() const
288 : {
289 32 : if (_n_rows != _n_cols)
290 6 : mooseError("ColumnMajorMatrix error: Unable to perform the operation on a non-square matrix.");
291 26 : }
292 :
293 : template <typename T>
294 : void
295 20 : ColumnMajorMatrixTempl<T>::checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const
296 : {
297 20 : if (_n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
298 4 : mooseError("ColumnMajorMatrix error: Unable to perform the operation on matrices of different "
299 : "shapes.");
300 16 : }
301 :
302 : template <>
303 : void
304 0 : ColumnMajorMatrixTempl<ADReal>::inverse(ColumnMajorMatrixTempl<ADReal> &) const
305 : {
306 0 : mooseError("Inverse solves with AD types is not supported for the ColumnMajorMatrix class.");
307 : }
308 :
309 : template <typename T>
310 : inline ColumnMajorMatrixTempl<T>
311 0 : ColumnMajorMatrixTempl<T>::abs()
312 : {
313 : using std::abs;
314 0 : ColumnMajorMatrixTempl<T> & s = (*this);
315 :
316 0 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
317 :
318 0 : for (unsigned int j = 0; j < _n_cols; j++)
319 0 : for (unsigned int i = 0; i < _n_rows; i++)
320 0 : ret_matrix(i, j) = abs(s(i, j));
321 :
322 0 : return ret_matrix;
323 0 : }
324 :
325 : template <typename T>
326 : inline T
327 2 : ColumnMajorMatrixTempl<T>::norm()
328 : {
329 : using std::sqrt;
330 2 : return sqrt(doubleContraction(*this));
331 : }
332 :
333 : template class ColumnMajorMatrixTempl<Real>;
334 : template class ColumnMajorMatrixTempl<ADReal>;
|
