MFEM  v4.5.2
Finite element discretization library
taddensemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
18 namespace 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.*/
33 template<typename dtype>
35 {
36 private:
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 
42 public:
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  double *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.
118  explicit TAutoDiffDenseMatrix(int 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.
135  TAutoDiffDenseMatrix(int m, int 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  double 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  const double max = (double) (RAND_MAX) + 1.;
304 
305  if (seed == 0)
306  {
307  seed = (int) time(0);
308  }
309 
310  // srand(seed++);
311  srand((unsigned) seed);
312 
313  for (int i = 0; i < capacity; i++)
314  {
315  data[i] = (dtype)(std::abs(rand() / max));
316  }
317  }
318 
319  void RandomizeDiag(int seed)
320  {
321  // static unsigned int seed = time(0);
322  const double max = (double) (RAND_MAX) + 1.;
323 
324  if (seed == 0)
325  {
326  seed = (int) time(0);
327  }
328 
329  // srand(seed++);
330  srand((unsigned) seed);
331 
332  for (int i = 0; i < std::min(height, width); i++)
333  {
334  Elem(i, i) = (dtype)(std::abs(rand() / max));
335  }
336  }
337 
338  /// Creates n x n diagonal matrix with diagonal elements c
339  void Diag(dtype c, int n)
340  {
341  SetSize(n);
342 
343  const int N = n * n;
344  for (int i = 0; i < N; i++)
345  {
346  data[i] = (dtype) 0.0;
347  }
348  for (int i = 0; i < n; i++)
349  {
350  data[i * (n + 1)] = c;
351  }
352  }
353  /// Creates n x n diagonal matrix with diagonal given by diag
354  template<typename itype>
355  void Diag(itype *diag, int n)
356  {
357  SetSize(n);
358 
359  int i, N = n * n;
360  for (i = 0; i < N; i++)
361  {
362  data[i] = 0.0;
363  }
364  for (i = 0; i < n; i++)
365  {
366  data[i * (n + 1)] = (dtype) diag[i];
367  }
368  }
369 
370  /// (*this) = (*this)^t
371  void Transpose()
372  {
373  int i, j;
374  dtype t;
375 
376  if (Width() == Height())
377  {
378  for (i = 0; i < Height(); i++)
379  for (j = i + 1; j < Width(); j++)
380  {
381  t = (*this)(i, j);
382  (*this)(i, j) = (*this)(j, i);
383  (*this)(j, i) = t;
384  }
385  }
386  else
387  {
388  TAutoDiffDenseMatrix<dtype> T(*this, 't');
389  (*this) = T;
390  }
391  }
392  /// (*this) = A^t
393  template<typename itype>
395  {
396  SetSize(A.Width(), A.Height());
397 
398  for (int i = 0; i < Height(); i++)
399  for (int j = 0; j < Width(); j++)
400  {
401  (*this)(i, j) = (dtype) A(j, i);
402  }
403  }
404 
405  /// (*this) = 1/2 ((*this) + (*this)^t)
406  void Symmetrize()
407  {
408 #ifdef MFEM_DEBUG
409  if (Width() != Height())
410  {
411  mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
412  }
413 #endif
414 
415  for (int i = 0; i < Height(); i++)
416  for (int j = 0; j < i; j++)
417  {
418  dtype a = 0.5 * ((*this)(i, j) + (*this)(j, i));
419  (*this)(j, i) = (*this)(i, j) = a;
420  }
421  }
422 
423  void Lump()
424  {
425  for (int i = 0; i < Height(); i++)
426  {
427  dtype L = 0.0;
428  for (int j = 0; j < Width(); j++)
429  {
430  L += (*this)(i, j);
431  (*this)(i, j) = (dtype) 0.0;
432  }
433  (*this)(i, i) = L;
434  }
435  }
436 };
437 
438 template<typename dtype>
441 {
442 #ifdef MFEM_DEBUG
443  if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
444  {
445  mfem_error("CalcAdjugate(...)");
446  }
447  if (a.Width() != adja.Height() || a.Height() != adja.Width())
448  {
449  mfem_error("CalcAdjugate(...)");
450  }
451 #endif
452 
453  if (a.Width() < a.Height())
454  {
455  const dtype *d = a.Data();
456  dtype *ad = adja.Data();
457  if (a.Width() == 1)
458  {
459  // N x 1, N = 2,3
460  ad[0] = d[0];
461  ad[1] = d[1];
462  if (a.Height() == 3)
463  {
464  ad[2] = d[2];
465  }
466  }
467  else
468  {
469  // 3 x 2
470  double e, g, f;
471  e = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
472  g = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
473  f = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
474 
475  ad[0] = d[0] * g - d[3] * f;
476  ad[1] = d[3] * e - d[0] * f;
477  ad[2] = d[1] * g - d[4] * f;
478  ad[3] = d[4] * e - d[1] * f;
479  ad[4] = d[2] * g - d[5] * f;
480  ad[5] = d[5] * e - d[2] * f;
481  }
482  return;
483  }
484 
485  if (a.Width() == 1)
486  {
487  adja(0, 0) = (dtype) 1.0;
488  }
489  else if (a.Width() == 2)
490  {
491  adja(0, 0) = a(1, 1);
492  adja(0, 1) = -a(0, 1);
493  adja(1, 0) = -a(1, 0);
494  adja(1, 1) = a(0, 0);
495  }
496  else
497  {
498  adja(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
499  adja(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
500  adja(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
501 
502  adja(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
503  adja(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
504  adja(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
505 
506  adja(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
507  adja(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
508  adja(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
509  }
510 }
511 
512 } // namespace mfem
513 
514 #endif
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: taddensemat.hpp:53
const dtype & Elem(int i, int j) const
void SetSize(int h, int w)
Change the size of the DenseMatrix to h x w.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2479
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void MultTranspose(const dtype *x, dtype *y) const
const dtype & operator()(int i, int j) const
void MultTranspose(const TAutoDiffVector< dtype > &x, TAutoDiffVector< dtype > &y) const
TAutoDiffDenseMatrix(const TAutoDiffDenseMatrix< dtype > &mat, char ch)
dtype * Data() const
Returns the matrix data array.
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
dtype operator*(const TAutoDiffDenseMatrix< dtype > &m) const
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
void Mult(const dtype *x, dtype *y) const
void Diag(dtype c, int n)
Creates n x n diagonal matrix with diagonal elements c.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: taddensemat.hpp:50
dtype * GetData() const
Returns the matrix data array.
dtype & Elem(int i, int j)
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: taddensemat.hpp:47
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: taddensemat.hpp:44
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
TAutoDiffDenseMatrix(int s)
Creates square matrix of size s.
int Size() const
Returns the size of the vector.
Definition: tadvector.hpp:246
double a
Definition: lissajous.cpp:41
TAutoDiffDenseMatrix(int m, int n)
Creates rectangular matrix of size m x n.
TAutoDiffDenseMatrix(const TAutoDiffDenseMatrix< idtype > &m)
Copy constructor.
Definition: taddensemat.hpp:67
RefCoord t[3]
Templated vector data type.
Definition: tadvector.hpp:40
void Diag(itype *diag, int n)
Creates n x n diagonal matrix with diagonal given by diag.
void Transpose(const TAutoDiffDenseMatrix< itype > &A)
(*this) = A^t
void Transpose()
(*this) = (*this)^t
void RandomizeDiag(int seed)
RefCoord s[3]
Templated dense matrix data type.
Definition: taddensemat.hpp:34
TAutoDiffDenseMatrix(const DenseMatrix &m)
Copy constructor using standard DenseMatrix.
Definition: taddensemat.hpp:92
void Mult(const TAutoDiffVector< dtype > &x, TAutoDiffVector< dtype > &y) const
dtype & operator()(int i, int j)
Returns reference to a_{ij}.