inv.hpp

Go to the documentation of this file.

// Copyright (c) 2019--present, The Simons Foundation
// This file is part of TRIQS/nda and is licensed under the Apache License, Version 2.0.
// SPDX-License-Identifier: Apache-2.0
// See LICENSE in the root of this distribution for details.

 
#pragma once
 
#include "../basic_array.hpp"
#include "../basic_functions.hpp"
#include "../blas/tools.hpp"
#include "../concepts.hpp"
#include "../exceptions.hpp"
#include "../lapack/getrf.hpp"
#include "../lapack/getri.hpp"
#include "../lapack/getrs.hpp"
#include "../layout/policies.hpp"
#include "../macros.hpp"
#include "../matrix_functions.hpp"
#include "../mem/address_space.hpp"
#include "../mem/policies.hpp"
#include "../traits.hpp"
 
#include <utility>
 
namespace nda::linalg {

 
  namespace detail {
 
    // Compute the inverse of a 2x2 matrix in place.
    void inv_in_place_2d(MemoryMatrix auto &&m) { // NOLINT (temporary views are allowed here)
      using value_t = get_value_t<decltype(m)>;
 
      // calculate the determinant of the matrix
      auto const det = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0));
      if (det == value_t{0.0}) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv_in_place: Matrix is not invertible";
      auto const detinv = value_t{1.0} / det;
 
      // multiply the adjoint by the inverse determinant
      std::swap(m(0, 0), m(1, 1));
      m(0, 0) *= +detinv;
      m(1, 1) *= +detinv;
      m(1, 0) *= -detinv;
      m(0, 1) *= -detinv;
    }
 
    // Compute the inverse of a 3x3 matrix in place.
    void inv_in_place_3d(MemoryMatrix auto &&m) { // NOLINT (temporary views are allowed here)
      using value_t = get_value_t<decltype(m)>;
      EXPECTS(is_matrix_square(m) and m.extent(0) == 3);
 
      // calculate the adjoint of the matrix
      auto adj  = stack_array<get_value_t<decltype(m)>, 3, 3>();
      adj(0, 0) = +m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1);
      adj(1, 0) = -m(1, 0) * m(2, 2) + m(1, 2) * m(2, 0);
      adj(2, 0) = +m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0);
      adj(0, 1) = -m(0, 1) * m(2, 2) + m(0, 2) * m(2, 1);
      adj(1, 1) = +m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0);
      adj(2, 1) = -m(0, 0) * m(2, 1) + m(0, 1) * m(2, 0);
      adj(0, 2) = +m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1);
      adj(1, 2) = -m(0, 0) * m(1, 2) + m(0, 2) * m(1, 0);
      adj(2, 2) = +m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
 
      // calculate the determinant of the matrix
      auto const det = m(0, 0) * adj(0, 0) + m(0, 1) * adj(1, 0) + m(0, 2) * adj(2, 0);
      if (det == value_t{0.0}) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv_in_place: Matrix is not invertible";
      auto const detinv = value_t{1.0} / det;
 
      // multiply the adjoint by the inverse determinant
      m = detinv * adj;
    }
 
  } // namespace detail

  template <blas_lapack::BlasArray<2> M>
    requires(get_algebra<M> == 'M' and mem::have_host_compatible_addr_space<M>)
  void inv_in_place(M &&m) { // NOLINT (temporary views are allowed here)
    using value_t = get_value_t<M>;
    EXPECTS(is_matrix_square(m));
 
    // use optimized routines for small matrices, otherwise use LAPACK routines
    auto const dim = m.shape()[0];
    if (dim == 1) {
      if (m(0, 0) == value_t{0.0}) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv_in_place: Matrix is not invertible";
      m(0, 0) = value_t{1.0} / m(0, 0);
    } else if (dim == 2) {
      detail::inv_in_place_2d(m);
    } else if (dim == 3) {
      detail::inv_in_place_3d(m);
    } else if (dim > 3) {
      // LU factorization with getrf
      auto ipiv = vector<int, heap<mem::get_addr_space<M>>>(dim);
      int info  = lapack::getrf(m, ipiv);
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv_in_place: getrf routine failed: info = " << info;
 
      // calculate the inverse with getri
      info = lapack::getri(m, ipiv);
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv_in_place: getri routine failed: info = " << info;
    }
  }

  template <Matrix M>
    requires(get_algebra<M> == 'M' and is_blas_lapack_v<get_value_t<M>>)
  auto inv(M const &m) {
    EXPECTS(is_matrix_square(m));
 
    // make a copy of the input matrix/view
    auto m_copy = make_regular(m);
 
    // for device compatible address spaces, we use getrf and getrs, otherwise we call inv_in_place
    if constexpr (mem::have_device_compatible_addr_space<M>) {
      // LU factorization with getrf
      auto ipiv = vector<int, heap<mem::get_addr_space<M>>>(m_copy.extent(0));
      int info  = lapack::getrf(m_copy, ipiv);
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv: getrf routine failed: info = " << info;
 
      // calculate the inverse with getrs and the identity matrix
      auto B = matrix<get_value_t<M>, F_layout, heap<mem::get_addr_space<M>>>{transpose(eye<get_value_t<M>>(m_copy.extent(0)))};
      info   = lapack::getrs(m_copy, B, ipiv);
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::inv: getrs routine failed: info = " << info;
      return B;
    } else {
      inv_in_place(m_copy);
      return m_copy;
    }
  }

  CLEF_MAKE_FNT_LAZY(inv)
 
  
 
} // namespace nda::linalg