MatrixTools Namespace Reference

Functions
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. More...
void	inverse (std::vector< PetscScalar > &A, unsigned int n)
	Inverts the dense "matrix" A using LAPACK routines. More...

Function Documentation

inverse() [1/2]

void MatrixTools::inverse	(	const std::vector< std::vector< Real >> &	m,
		std::vector< std::vector< Real >> &	m_inv
	)

Inverse the dense square matrix m using LAPACK routines.

If you need to invert a matrix "in place", make the two arguments the same @ param m The matrix to invert @ param[out] m_inv The inverse of m, which must be of the same size as m. @ return if zero then the inversion was successful. Otherwise m was not square, contained illegal entries or was singular

Definition at line 23 of file MatrixTools.C.

Referenced by RankFourTensorTempl< T >::invSymm(), Moose::SlepcSupport::mooseSlepcEPSMonitor(), and EigenProblem::outputInverseEigenvalue().

{
  unsigned int n = m.size();

  // check the matrix m exists and is square
  if (n == 0)
    throw MooseException("Input matrix empty during matrix inversion.");
  if (n != m_inv.size() || n != m[0].size() || n != m_inv[0].size())
    throw MooseException("Input and output matrix are not same size square matrices.");

  // build the vectorial representation
  std::vector<PetscScalar> A;
  for (const auto & rowvec : m)
    for (const auto & matrix_entry : rowvec)
      A.push_back(matrix_entry);

  inverse(A, n);

  // build the inverse
  unsigned int i = 0;
  for (auto & rowvec : m_inv)
    for (auto & inv_entry : rowvec)
      inv_entry = A[i++];
}

inverse() [2/2]

void MatrixTools::inverse	(	std::vector< PetscScalar > &	A,
		unsigned int	n
	)

Inverts the dense "matrix" A using LAPACK routines.

Parameters

A	upon input this is a row vector representing a square matrix of size sqrt(n)*sqrt(n). Upon output it is the inverse (as a row-vector)
n	size of the vector A

Returns

if zero then inversion was successful. Otherwise A contained illegal entries or was singular

Definition at line 49 of file MatrixTools.C.

{
  mooseAssert(n >= 1, "MatrixTools::inverse - n (leading dimension) needs to be positive");
  mooseAssert(n <= std::numeric_limits<unsigned int>::max(),
              "MatrixTools::inverse - n (leading dimension) too large");

  std::vector<PetscBLASInt> ipiv(n);
  std::vector<PetscScalar> buffer(n * 64);

  // Following does a LU decomposition of "square matrix A"
  // upon return "A = P*L*U" if return_value == 0
  // Here I use quotes because A is actually an array of length n^2, not a matrix of size n-by-n
  PetscBLASInt return_value;
  LAPACKgetrf_(reinterpret_cast<PetscBLASInt *>(&n),
               reinterpret_cast<PetscBLASInt *>(&n),
               &A[0],
               reinterpret_cast<PetscBLASInt *>(&n),
               &ipiv[0],
               &return_value);

  if (return_value != 0)
    throw MooseException(
        return_value < 0
            ? "Argument " + Moose::stringify(-return_value) +
                  " was invalid during LU factorization in MatrixTools::inverse."
            : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
                  " was exactly zero during LU factorization in MatrixTools::inverse.");

  // get the inverse of A
  PetscBLASInt buffer_size = buffer.size();
  LAPACKgetri_(reinterpret_cast<PetscBLASInt *>(&n),
               &A[0],
               reinterpret_cast<PetscBLASInt *>(&n),
               &ipiv[0],
               &buffer[0],
               &buffer_size,
               &return_value);

  if (return_value != 0)
    throw MooseException(return_value < 0
                             ? "Argument " + Moose::stringify(-return_value) +
                                   " was invalid during invert in MatrixTools::inverse."
                             : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
                                   " was exactly zero during invert in MatrixTools::inverse.");
}