getri.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2021-2023, University of Colorado Denver. All rights reserved.
//
// This file is part of <T>LAPACK.
// <T>LAPACK is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
 
#ifndef TLAPACK_GETRI_HH
#define TLAPACK_GETRI_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/blas/swap.hpp"
#include "tlapack/lapack/getri_uili.hpp"
#include "tlapack/lapack/getri_uxli.hpp"
 
namespace tlapack {
 
enum class GetriVariant : char {
    UILI = 'D',  
    UXLI = 'C'   
};
 
struct GetriOpts {
    GetriVariant variant = GetriVariant::UILI;
};
 
template <class T, TLAPACK_SMATRIX matrix_t, TLAPACK_VECTOR piv_t>
constexpr WorkInfo getri_worksize(const matrix_t& A,
                                  const piv_t& piv,
                                  const GetriOpts& opts = {})
{
    if (opts.variant == GetriVariant::UXLI)
        return getri_uxli_worksize<T>(A, opts);
 
    return WorkInfo(0);
}
 
template <TLAPACK_SMATRIX matrix_t,
          TLAPACK_VECTOR piv_t,
          TLAPACK_WORKSPACE work_t>
int getri_work(matrix_t& A,
               const piv_t& piv,
               work_t& work,
               const GetriOpts& opts = {})
{
    using idx_t = size_type<matrix_t>;
 
    // Constants
    const idx_t n = ncols(A);
 
    // Call variant
    int info;
    if (opts.variant == GetriVariant::UXLI)
        info = getri_uxli_work(A, work);
    else
        info = getri_uili(A);
 
    // Return is matrix is not invertible
    if (info != 0) return info;
 
    // swap columns of X to find A^{-1} since A^{-1}=X P
    for (idx_t j = n; j-- > 0;) {
        if (piv[j] != j) {
            auto vect1 = tlapack::col(A, j);
            auto vect2 = tlapack::col(A, piv[j]);
            tlapack::swap(vect1, vect2);
        }
    }
 
    return 0;
}
 
template <TLAPACK_SMATRIX matrix_t, TLAPACK_VECTOR piv_t>
int getri(matrix_t& A, const piv_t& piv, const GetriOpts& opts = {})
{
    using idx_t = size_type<matrix_t>;
 
    // Constants
    const idx_t n = ncols(A);
 
    // Call variant
    int info;
    if (opts.variant == GetriVariant::UXLI)
        info = getri_uxli(A);
    else
        info = getri_uili(A);
 
    // Return is matrix is not invertible
    if (info != 0) return info;
 
    // swap columns of X to find A^{-1} since A^{-1}=X P
    for (idx_t j = n; j-- > 0;) {
        if (piv[j] != j) {
            auto vect1 = tlapack::col(A, j);
            auto vect2 = tlapack::col(A, piv[j]);
            tlapack::swap(vect1, vect2);
        }
    }
 
    return 0;
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_GETRI_HH