det.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 "../concepts.hpp"
#include "../exceptions.hpp"
#include "../lapack/getrf.hpp"
#include "../matrix_functions.hpp"
#include "../mem/address_space.hpp"
#include "../mem/policies.hpp"
#include "../traits.hpp"
 
namespace nda::linalg {


  template <Matrix M>
    requires(get_algebra<M> == 'M' and nda::mem::have_host_compatible_addr_space<M> and is_blas_lapack_v<get_value_t<M>>)
  auto det_in_place(M &&m) { // NOLINT (temporary views are allowed here)
    EXPECTS(is_matrix_square(m));
 
    // use optimized routines for small matrices, otherwise use LAPACK routine
    auto const dim = m.shape()[0];
    if (dim == 1) {
      return m(0, 0);
    } else if (dim == 2) {
      return m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
    } else if (dim == 3) {
      return m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) //
         - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0))   //
         + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));  //
    } else if (dim > 3) {
      // LU factorization with getrf
      auto ipiv = vector<int, sso<100>>(dim);
      int info  = nda::lapack::getrf(m, ipiv);
      if (info < 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::det_in_place: getrf routine failed: info = " << info;
 
      // calculate the determinant from the LU decomposition
      auto det    = get_value_t<M>{1};
      int n_flips = 0;
      for (int i = 0; i < dim; i++) {
        det *= m(i, i);
        // count the number of column interchanges performed by getrf
        if (ipiv(i) != i + 1) ++n_flips;
      }
 
      return ((n_flips % 2 == 1) ? -det : det);
    } else {
      // empty matrix
      return get_value_t<M>{1};
    }
  }

  template <Matrix M>
    requires(get_algebra<M> == 'M' and nda::mem::have_host_compatible_addr_space<M> and is_blas_lapack_v<get_value_t<M>>)
  auto det(M const &m) {
    EXPECTS(is_matrix_square(m));
    auto m_copy = make_regular(m);
    return det_in_place(m_copy);
  }

 
} // namespace nda::linalg
 
namespace nda::clef {
 
  // Make nda::linalg::det visible to lazy function calls.
  using nda::linalg::det;

  CLEF_MAKE_FNT_LAZY(det)
 
} // namespace nda::clef