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             : #include "MatrixTools.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Conversion.h"
      14             : #include "MooseError.h"
      15             : #include "MooseException.h"
      16             : 
      17             : // PETSc includes
      18             : #include "petscblaslapack.h"
      19             : 
      20             : namespace MatrixTools
      21             : {
      22             : void
      23          10 : inverse(const std::vector<std::vector<Real>> & m, std::vector<std::vector<Real>> & m_inv)
      24             : {
      25          10 :   unsigned int n = m.size();
      26             : 
      27             :   // check the matrix m exists and is square
      28          10 :   if (n == 0)
      29           2 :     throw MooseException("Input matrix empty during matrix inversion.");
      30           8 :   if (n != m_inv.size() || n != m[0].size() || n != m_inv[0].size())
      31           2 :     throw MooseException("Input and output matrix are not same size square matrices.");
      32             : 
      33             :   // build the vectorial representation
      34           6 :   std::vector<PetscScalar> A;
      35          22 :   for (const auto & rowvec : m)
      36          60 :     for (const auto & matrix_entry : rowvec)
      37          44 :       A.push_back(matrix_entry);
      38             : 
      39           6 :   inverse(A, n);
      40             : 
      41             :   // build the inverse
      42           4 :   unsigned int i = 0;
      43          14 :   for (auto & rowvec : m_inv)
      44          36 :     for (auto & inv_entry : rowvec)
      45          26 :       inv_entry = A[i++];
      46           6 : }
      47             : 
      48             : void
      49          18 : inverse(std::vector<PetscScalar> & A, unsigned int n)
      50             : {
      51             :   mooseAssert(n >= 1, "MatrixTools::inverse - n (leading dimension) needs to be positive");
      52             :   mooseAssert(n <= std::numeric_limits<unsigned int>::max(),
      53             :               "MatrixTools::inverse - n (leading dimension) too large");
      54             : 
      55          36 :   std::vector<PetscBLASInt> ipiv(n);
      56          18 :   std::vector<PetscScalar> buffer(n * 64);
      57             : 
      58             :   // Following does a LU decomposition of "square matrix A"
      59             :   // upon return "A = P*L*U" if return_value == 0
      60             :   // Here I use quotes because A is actually an array of length n^2, not a matrix of size n-by-n
      61             :   PetscBLASInt return_value;
      62          18 :   LAPACKgetrf_(reinterpret_cast<PetscBLASInt *>(&n),
      63             :                reinterpret_cast<PetscBLASInt *>(&n),
      64          18 :                &A[0],
      65             :                reinterpret_cast<PetscBLASInt *>(&n),
      66          18 :                &ipiv[0],
      67             :                &return_value);
      68             : 
      69          18 :   if (return_value != 0)
      70             :     throw MooseException(
      71           4 :         return_value < 0
      72          12 :             ? "Argument " + Moose::stringify(-return_value) +
      73             :                   " was invalid during LU factorization in MatrixTools::inverse."
      74           8 :             : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
      75          12 :                   " was exactly zero during LU factorization in MatrixTools::inverse.");
      76             : 
      77             :   // get the inverse of A
      78          14 :   PetscBLASInt buffer_size = buffer.size();
      79          14 :   LAPACKgetri_(reinterpret_cast<PetscBLASInt *>(&n),
      80          14 :                &A[0],
      81             :                reinterpret_cast<PetscBLASInt *>(&n),
      82          14 :                &ipiv[0],
      83          14 :                &buffer[0],
      84             :                &buffer_size,
      85             :                &return_value);
      86             : 
      87          14 :   if (return_value != 0)
      88           0 :     throw MooseException(return_value < 0
      89           0 :                              ? "Argument " + Moose::stringify(-return_value) +
      90             :                                    " was invalid during invert in MatrixTools::inverse."
      91           0 :                              : "Matrix on-diagonal entry " + Moose::stringify(return_value) +
      92           0 :                                    " was exactly zero during invert in MatrixTools::inverse.");
      93          22 : }
      94             : }