MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
taddensemat.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef TADDENSEMATRIX_H
13#define TADDENSEMATRIX_H
14
15#include "mfem.hpp"
16#include "tadvector.hpp"
17
18namespace mfem
19{
20/// Templated dense matrix data type.
21/** The main goal of the TAutoDiffDenseMatrix class is to serve as a data
22 container for representing dense matrices in classes, methods, and functions
23 utilized with automatic differentiation (AD). The functionality/interface is
24 copied from the standard MFEM dense matrix mfem::DenseMatrix. The basic idea
25 is to utilize the templated vector class in combination with AD during the
26 development phase. The AD parts can be replaced with optimized code once the
27 initial development of the application is complete. The common interface
28 between TAutoDiffDenseMatrix and DenseMatrix will ease the transition from
29 AD to hand-optimized code as it does not require a change in the interface
30 or the code structure. TAutoDiffDenseMatrix is intended to be utilized for
31 dense serial matrices. The objects can be combined with TAutoDiffVector or
32 standard Vector.*/
33template<typename dtype>
35{
36private:
37 int height; ///< Dimension of the output / number of rows in the matrix.
38 int width; ///< Dimension of the input / number of columns in the matrix.
39 dtype *data;
40 int capacity; // zero or negative capacity means we do not own the data.
41
42public:
43 /// Get the height (size of output) of the Operator. Synonym with NumRows().
44 inline int Height() const { return height; }
45 /** @brief Get the number of rows (size of output) of the Operator. Synonym
46 with Height(). */
47 inline int NumRows() const { return height; }
48
49 /// Get the width (size of input) of the Operator. Synonym with NumCols().
50 inline int Width() const { return width; }
51 /** @brief Get the number of columns (size of input) of the Operator. Synonym
52 with Width(). */
53 inline int NumCols() const { return width; }
54
55 /** Default constructor for TAutoDiffDenseMatrix.
56 Sets data = NULL and height = width = 0. */
58 {
59 data = nullptr;
60 capacity = 0;
61 height = 0;
62 width = 0;
63 }
64
65 /// Copy constructor
66 template<typename idtype>
68 {
69 height = m.GetHeight();
70 width = m.GetWidth();
71 const int hw = height * width;
72 if (hw > 0)
73 {
74 idtype *mdata = m.Data();
75 MFEM_ASSERT(mdata, "invalid source matrix");
76 data = new dtype[hw];
77 capacity = hw;
78 for (int i = 0; i < hw; i++)
79 {
80 data[i] = mdata[i];
81 }
82 }
83 else
84 {
85 data = nullptr;
86 capacity = 0;
87 width = 0;
88 height = 0;
89 }
90 }
91 /// Copy constructor using standard DenseMatrix
93 {
94 height = m.Height();
95 width = m.Width();
96 const int hw = height * width;
97 if (hw > 0)
98 {
99 real_t *mdata = m.Data();
100 MFEM_ASSERT(mdata, "invalid source matrix");
101 data = new dtype[hw];
102 capacity = hw;
103 for (int i = 0; i < hw; i++)
104 {
105 data[i] = mdata[i];
106 }
107 }
108 else
109 {
110 data = nullptr;
111 capacity = 0;
112 width = 0;
113 height = 0;
114 }
115 }
116
117 /// Creates square matrix of size s.
119 {
120 MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
121 height = s;
122 width = s;
123 capacity = s * s;
124 if (capacity > 0)
125 {
126 data = new dtype[capacity](); // init with zeroes
127 }
128 else
129 {
130 data = nullptr;
131 }
132 }
133
134 /// Creates rectangular matrix of size m x n.
136 {
137 MFEM_ASSERT(m >= 0 && n >= 0,
138 "invalid DenseMatrix size: " << m << " x " << n);
139 height = m;
140 width = n;
141 capacity = m * n;
142 if (capacity > 0)
143 {
144 data = new dtype[capacity](); // init with zeroes
145 }
146 else
147 {
148 data = nullptr;
149 }
150 }
151
153 {
154 height = mat.Width();
155 width = mat.Height();
156 capacity = height * width;
157 if (capacity > 0)
158 {
159 data = new dtype[capacity];
160
161 for (int i = 0; i < height; i++)
162 {
163 for (int j = 0; j < width; j++)
164 {
165 (*this)(i, j) = mat(j, i);
166 }
167 }
168 }
169 else
170 {
171 data = nullptr;
172 }
173 }
174
175 /// Change the size of the DenseMatrix to s x s.
176 void SetSize(int s) { SetSize(s, s); }
177
178 /// Change the size of the DenseMatrix to h x w.
179 void SetSize(int h, int w)
180 {
181 MFEM_ASSERT(h >= 0 && w >= 0,
182 "invalid DenseMatrix size: " << h << " x " << w);
183 if (Height() == h && Width() == w)
184 {
185 return;
186 }
187 height = h;
188 width = w;
189 const int hw = h * w;
190 if (hw > std::abs(capacity))
191 {
192 if (capacity > 0)
193 {
194 delete[] data;
195 }
196 capacity = hw;
197 data = new dtype[hw](); // init with zeroes
198 }
199 }
200
201 /// Returns the matrix data array.
202 inline dtype *Data() const { return data; }
203 /// Returns the matrix data array.
204 inline dtype *GetData() const { return data; }
205
206 inline bool OwnsData() const { return (capacity > 0); }
207
208 /// Returns reference to a_{ij}.
209 dtype &operator()(int i, int j)
210 {
211 MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
212 return data[i + j * height];
213 }
214
215 const dtype &operator()(int i, int j) const
216 {
217 MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
218 return data[i + j * height];
219 }
220
221 dtype &Elem(int i, int j) { return (*this)(i, j); }
222
223 const dtype &Elem(int i, int j) const { return (*this)(i, j); }
224
225 void Mult(const dtype *x, dtype *y) const
226 {
227 if (width == 0)
228 {
229 for (int row = 0; row < height; row++)
230 {
231 y[row] = 0.0;
232 }
233 return;
234 }
235 dtype *d_col = data;
236 dtype x_col = x[0];
237 for (int row = 0; row < height; row++)
238 {
239 y[row] = x_col * d_col[row];
240 }
241 d_col += height;
242 for (int col = 1; col < width; col++)
243 {
244 x_col = x[col];
245 for (int row = 0; row < height; row++)
246 {
247 y[row] += x_col * d_col[row];
248 }
249 d_col += height;
250 }
251 }
252
254 {
255 MFEM_ASSERT(height == y.Size() && width == x.Size(),
256 "incompatible dimensions");
257
258 Mult((const dtype *) x, (dtype *) y);
259 }
260
262 {
263 MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
264 "incompatible dimensions");
265
266 const int hw = height * width;
267 dtype a = 0.0;
268 for (int i = 0; i < hw; i++)
269 {
270 a += data[i] * m.data[i];
271 }
272
273 return a;
274 }
275
276 void MultTranspose(const dtype *x, dtype *y) const
277 {
278 dtype *d_col = data;
279 for (int col = 0; col < width; col++)
280 {
281 real_t y_col = 0.0;
282 for (int row = 0; row < height; row++)
283 {
284 y_col += x[row] * d_col[row];
285 }
286 y[col] = y_col;
287 d_col += height;
288 }
289 }
290
292 TAutoDiffVector<dtype> &y) const
293 {
294 MFEM_ASSERT(height == x.Size() && width == y.Size(),
295 "incompatible dimensions");
296
297 MultTranspose((const dtype *) x, (dtype *) y);
298 }
299
300 void Randomize(int seed)
301 {
302 // static unsigned int seed = time(0);
303
304 if (seed == 0)
305 {
306 seed = (int) time(0);
307 }
308
309 // srand(seed++);
310 srand((unsigned) seed);
311
312 for (int i = 0; i < capacity; i++)
313 {
314 data[i] = (dtype)(rand_real());
315 }
316 }
317
318 void RandomizeDiag(int seed)
319 {
320 // static unsigned int seed = time(0);
321
322 if (seed == 0)
323 {
324 seed = (int) time(0);
325 }
326
327 // srand(seed++);
328 srand((unsigned) seed);
329
330 for (int i = 0; i < std::min(height, width); i++)
331 {
332 Elem(i, i) = (dtype)(rand_real());
333 }
334 }
335
336 /// Creates n x n diagonal matrix with diagonal elements c
337 void Diag(dtype c, int n)
338 {
339 SetSize(n);
340
341 const int N = n * n;
342 for (int i = 0; i < N; i++)
343 {
344 data[i] = (dtype) 0.0;
345 }
346 for (int i = 0; i < n; i++)
347 {
348 data[i * (n + 1)] = c;
349 }
350 }
351 /// Creates n x n diagonal matrix with diagonal given by diag
352 template<typename itype>
353 void Diag(itype *diag, int n)
354 {
355 SetSize(n);
356
357 int i, N = n * n;
358 for (i = 0; i < N; i++)
359 {
360 data[i] = 0.0;
361 }
362 for (i = 0; i < n; i++)
363 {
364 data[i * (n + 1)] = (dtype) diag[i];
365 }
366 }
367
368 /// (*this) = (*this)^t
370 {
371 int i, j;
372 dtype t;
373
374 if (Width() == Height())
375 {
376 for (i = 0; i < Height(); i++)
377 for (j = i + 1; j < Width(); j++)
378 {
379 t = (*this)(i, j);
380 (*this)(i, j) = (*this)(j, i);
381 (*this)(j, i) = t;
382 }
383 }
384 else
385 {
386 TAutoDiffDenseMatrix<dtype> T(*this, 't');
387 (*this) = T;
388 }
389 }
390 /// (*this) = A^t
391 template<typename itype>
393 {
394 SetSize(A.Width(), A.Height());
395
396 for (int i = 0; i < Height(); i++)
397 for (int j = 0; j < Width(); j++)
398 {
399 (*this)(i, j) = (dtype) A(j, i);
400 }
401 }
402
403 /// (*this) = 1/2 ((*this) + (*this)^t)
405 {
406#ifdef MFEM_DEBUG
407 if (Width() != Height())
408 {
409 mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
410 }
411#endif
412
413 for (int i = 0; i < Height(); i++)
414 for (int j = 0; j < i; j++)
415 {
416 dtype a = 0.5 * ((*this)(i, j) + (*this)(j, i));
417 (*this)(j, i) = (*this)(i, j) = a;
418 }
419 }
420
421 void Lump()
422 {
423 for (int i = 0; i < Height(); i++)
424 {
425 dtype L = 0.0;
426 for (int j = 0; j < Width(); j++)
427 {
428 L += (*this)(i, j);
429 (*this)(i, j) = (dtype) 0.0;
430 }
431 (*this)(i, i) = L;
432 }
433 }
434};
435
436template<typename dtype>
439{
440#ifdef MFEM_DEBUG
441 if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
442 {
443 mfem_error("CalcAdjugate(...)");
444 }
445 if (a.Width() != adja.Height() || a.Height() != adja.Width())
446 {
447 mfem_error("CalcAdjugate(...)");
448 }
449#endif
450
451 if (a.Width() < a.Height())
452 {
453 const dtype *d = a.Data();
454 dtype *ad = adja.Data();
455 if (a.Width() == 1)
456 {
457 // N x 1, N = 2,3
458 ad[0] = d[0];
459 ad[1] = d[1];
460 if (a.Height() == 3)
461 {
462 ad[2] = d[2];
463 }
464 }
465 else
466 {
467 // 3 x 2
468 real_t e, g, f;
469 e = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
470 g = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
471 f = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
472
473 ad[0] = d[0] * g - d[3] * f;
474 ad[1] = d[3] * e - d[0] * f;
475 ad[2] = d[1] * g - d[4] * f;
476 ad[3] = d[4] * e - d[1] * f;
477 ad[4] = d[2] * g - d[5] * f;
478 ad[5] = d[5] * e - d[2] * f;
479 }
480 return;
481 }
482
483 if (a.Width() == 1)
484 {
485 adja(0, 0) = (dtype) 1.0;
486 }
487 else if (a.Width() == 2)
488 {
489 adja(0, 0) = a(1, 1);
490 adja(0, 1) = -a(0, 1);
491 adja(1, 0) = -a(1, 0);
492 adja(1, 1) = a(0, 0);
493 }
494 else
495 {
496 adja(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
497 adja(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
498 adja(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
499
500 adja(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
501 adja(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
502 adja(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
503
504 adja(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
505 adja(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
506 adja(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
507 }
508}
509
510} // namespace mfem
511
512#endif
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:111
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Templated dense matrix data type.
TAutoDiffDenseMatrix(int m, int n)
Creates rectangular matrix of size m x n.
dtype * GetData() const
Returns the matrix data array.
void SetSize(int h, int w)
Change the size of the DenseMatrix to h x w.
dtype & Elem(int i, int j)
TAutoDiffDenseMatrix(const TAutoDiffDenseMatrix< idtype > &m)
Copy constructor.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void Diag(itype *diag, int n)
Creates n x n diagonal matrix with diagonal given by diag.
TAutoDiffDenseMatrix(const DenseMatrix &m)
Copy constructor using standard DenseMatrix.
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
dtype * Data() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
dtype operator*(const TAutoDiffDenseMatrix< dtype > &m) const
const dtype & Elem(int i, int j) const
void MultTranspose(const dtype *x, dtype *y) const
TAutoDiffDenseMatrix(const TAutoDiffDenseMatrix< dtype > &mat, char ch)
void Diag(dtype c, int n)
Creates n x n diagonal matrix with diagonal elements c.
const dtype & operator()(int i, int j) const
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void Mult(const dtype *x, dtype *y) const
TAutoDiffDenseMatrix(int s)
Creates square matrix of size s.
void Transpose()
(*this) = (*this)^t
void Transpose(const TAutoDiffDenseMatrix< itype > &A)
(*this) = A^t
void Mult(const TAutoDiffVector< dtype > &x, TAutoDiffVector< dtype > &y) const
void MultTranspose(const TAutoDiffVector< dtype > &x, TAutoDiffVector< dtype > &y) const
dtype & operator()(int i, int j)
Returns reference to a_{ij}.
Templated vector data type.
Definition tadvector.hpp:41
int Size() const
Returns the size of the vector.
IntegrationType itype
Definition ex38.cpp:46
real_t a
Definition lissajous.cpp:41
void mfem_error(const char *msg)
Definition error.cpp:154
real_t rand_real()
Generate a random real_t number in the interval [0,1) using rand().
Definition vector.hpp:59
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
RefCoord t[3]
RefCoord s[3]