MatrixTools.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7
8 #include "MatrixTools.h"
9
10 // MOOSE includes
11 #include "Conversion.h"
12 #include "MooseError.h"
13 #include "MooseException.h"
14
15 // PETSc includes
16 #include "petscblaslapack.h"
17
18 #if PETSC_VERSION_LESS_THAN(3, 5, 0)
19 extern "C" void FORTRAN_CALL(dgetri)(...); // matrix inversion routine from LAPACK
20 #endif
21
22 namespace MatrixTools
23 {
24 void
25 inverse(const std::vector<std::vector<Real>> & m, std::vector<std::vector<Real>> & m_inv)
26 {
27  unsigned int n = m.size();
28
29  // check the matrix m exists and is square
30  if (n == 0)
31  throw MooseException("Input matrix empty during matrix inversion.");
32  if (n != m_inv.size() || n != m[0].size() || n != m_inv[0].size())
33  throw MooseException("Input and output matrix are not same size square matrices.");
34
35  // build the vectorial representation
36  std::vector<PetscScalar> A;
37  for (const auto & rowvec : m)
38  for (const auto & matrix_entry : rowvec)
39  A.push_back(matrix_entry);
40
41  inverse(A, n);
42
43  // build the inverse
44  unsigned int i = 0;
45  for (auto & rowvec : m_inv)
46  for (auto & inv_entry : rowvec)
47  inv_entry = A[i++];
48 }
49
50 void
51 inverse(std::vector<PetscScalar> & A, unsigned int n)
52 {
53  mooseAssert(n >= 1, "MatrixTools::inverse - n (leading dimension) needs to be positive");
54  mooseAssert(n <= std::numeric_limits<int>::max(),
55  "MatrixTools::inverse - n (leading dimension) too large");
56
57  std::vector<PetscBLASInt> ipiv(n);
58  std::vector<PetscScalar> buffer(n * 64);
59
60  // Following does a LU decomposition of "square matrix A"
61  // upon return "A = P*L*U" if return_value == 0
62  // Here I use quotes because A is actually an array of length n^2, not a matrix of size n-by-n
63  int return_value;
64  LAPACKgetrf_(reinterpret_cast<int *>(&n),
65  reinterpret_cast<int *>(&n),
66  &A[0],
67  reinterpret_cast<int *>(&n),
68  &ipiv[0],
69  &return_value);
70
71  if (return_value != 0)
72  throw MooseException(
73  return_value < 0
74  ? "Argument " + Moose::stringify(-return_value) +
75  " was invalid during LU factorization in MatrixTools::inverse."
76  : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
77  " was exactly zero during LU factorization in MatrixTools::inverse.");
78
79  // get the inverse of A
80  int buffer_size = buffer.size();
81 #if PETSC_VERSION_LESS_THAN(3, 5, 0)
82  FORTRAN_CALL(dgetri)
83  (reinterpret_cast<int *>(&n),
84  &A[0],
85  reinterpret_cast<int *>(&n),
86  &ipiv[0],
87  &buffer[0],
88  &buffer_size,
89  &return_value);
90 #else
91  LAPACKgetri_(reinterpret_cast<int *>(&n),
92  &A[0],
93  reinterpret_cast<int *>(&n),
94  &ipiv[0],
95  &buffer[0],
96  &buffer_size,
97  &return_value);
98 #endif
99
100  if (return_value != 0)
101  throw MooseException(return_value < 0
102  ? "Argument " + Moose::stringify(-return_value) +
103  " was invalid during invert in MatrixTools::inverse."
104  : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
105  " was exactly zero during invert in MatrixTools::inverse.");
106 }
107 }
void inverse(const std::vector< std::vector< Real >> &m, std::vector< std::vector< Real >> &m_inv)
Inverse the dense square matrix m using LAPACK routines.
Definition: MatrixTools.C:25
void FORTRAN_CALL() dgetri(...)
static PetscErrorCode Mat * A
PetscInt m
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:66
Provides a way for users to bail out of the current solve.
PetscInt n