trsv.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2017-2021, University of Tennessee. All rights reserved.
// 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_BLAS_TRSV_HH
#define TLAPACK_BLAS_TRSV_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/lapack/conjugate.hpp"
 
namespace tlapack {
 
template <
    TLAPACK_MATRIX matrixA_t,
    TLAPACK_VECTOR vectorX_t,
    class T = type_t<vectorX_t>,
    disable_if_allow_optblas_t<pair<matrixA_t, T>, pair<vectorX_t, T> > = 0>
void trsv(Uplo uplo, Op trans, Diag diag, const matrixA_t& A, vectorX_t& x)
{
    // data traits
    using TA = type_t<matrixA_t>;
    using TX = type_t<vectorX_t>;
    using idx_t = size_type<matrixA_t>;
 
    // constants
    const idx_t n = nrows(A);
    const bool nonunit = (diag == Diag::NonUnit);
 
    // check arguments
    tlapack_check_false(uplo != Uplo::Lower && uplo != Uplo::Upper);
    tlapack_check_false(trans != Op::NoTrans && trans != Op::Trans &&
                        trans != Op::ConjTrans && trans != Op::Conj);
    tlapack_check_false(diag != Diag::NonUnit && diag != Diag::Unit);
    tlapack_check_false(nrows(A) != ncols(A));
    tlapack_check_false(size(x) != n);
 
    if (trans == Op::NoTrans) {
        // Form x := A^{-1} * x
        if (uplo == Uplo::Upper) {
            // upper
            for (idx_t j = n - 1; j != idx_t(-1); --j) {
                // note: NOT skipping if x[j] is zero, for consistent NAN
                // handling
                if (nonunit) {
                    x[j] /= A(j, j);
                }
                for (idx_t i = j - 1; i != idx_t(-1); --i) {
                    x[i] -= x[j] * A(i, j);
                }
            }
        }
        else {
            // lower
            for (idx_t j = 0; j < n; ++j) {
                // note: NOT skipping if x[j] is zero ...
                if (nonunit) {
                    x[j] /= A(j, j);
                }
                for (idx_t i = j + 1; i < n; ++i) {
                    x[i] -= x[j] * A(i, j);
                }
            }
        }
    }
    else if (trans == Op::Conj) {
        // Form x := A^{-1} * x
        if (uplo == Uplo::Upper) {
            // upper
            for (idx_t j = n - 1; j != idx_t(-1); --j) {
                // note: NOT skipping if x[j] is zero, for consistent NAN
                // handling
                if (nonunit) {
                    x[j] /= conj(A(j, j));
                }
                for (idx_t i = j - 1; i != idx_t(-1); --i) {
                    x[i] -= x[j] * conj(A(i, j));
                }
            }
        }
        else {
            // lower
            for (idx_t j = 0; j < n; ++j) {
                // note: NOT skipping if x[j] is zero ...
                if (nonunit) {
                    x[j] /= conj(A(j, j));
                }
                for (idx_t i = j + 1; i < n; ++i) {
                    x[i] -= x[j] * conj(A(i, j));
                }
            }
        }
    }
    else if (trans == Op::Trans) {
        // Form  x := A^{-T} * x
 
        using scalar_t = scalar_type<TA, TX>;
 
        if (uplo == Uplo::Upper) {
            // upper
            for (idx_t j = 0; j < n; ++j) {
                scalar_t tmp = x[j];
                for (idx_t i = 0; i < j; ++i) {
                    tmp -= A(i, j) * x[i];
                }
                if (nonunit) {
                    tmp /= A(j, j);
                }
                x[j] = tmp;
            }
        }
        else {
            // lower
            for (idx_t j = n - 1; j != idx_t(-1); --j) {
                scalar_t tmp = x[j];
                for (idx_t i = j + 1; i < n; ++i) {
                    tmp -= A(i, j) * x[i];
                }
                if (nonunit) {
                    tmp /= A(j, j);
                }
                x[j] = tmp;
            }
        }
    }
    else {
        // Form x := A^{-H} * x
        // same code as above A^{-T} * x case, except add conj()
 
        using scalar_t = scalar_type<TA, TX>;
 
        if (uplo == Uplo::Upper) {
            // upper
            for (idx_t j = 0; j < n; ++j) {
                scalar_t tmp = x[j];
                for (idx_t i = 0; i < j; ++i) {
                    tmp -= conj(A(i, j)) * x[i];
                }
                if (nonunit) {
                    tmp /= conj(A(j, j));
                }
                x[j] = tmp;
            }
        }
        else {
            // lower
            for (idx_t j = n - 1; j != idx_t(-1); --j) {
                scalar_t tmp = x[j];
                for (idx_t i = j + 1; i < n; ++i) {
                    tmp -= conj(A(i, j)) * x[i];
                }
                if (nonunit) {
                    tmp /= conj(A(j, j));
                }
                x[j] = tmp;
            }
        }
    }
}
 
#ifdef TLAPACK_USE_LAPACKPP
 
template <
    TLAPACK_LEGACY_MATRIX matrixA_t,
    TLAPACK_LEGACY_VECTOR vectorX_t,
    class T = type_t<vectorX_t>,
    enable_if_allow_optblas_t<pair<matrixA_t, T>, pair<vectorX_t, T> > = 0>
void trsv(Uplo uplo, Op trans, Diag diag, const matrixA_t& A, vectorX_t& x)
{
    // Legacy objects
    auto A_ = legacy_matrix(A);
    auto x_ = legacy_vector(x);
 
    // Constants to forward
    constexpr Layout L = layout<matrixA_t>;
    const auto& n = A_.n;
 
    if (trans != Op::Conj)
        ::blas::trsv((::blas::Layout)L, (::blas::Uplo)uplo, (::blas::Op)trans,
                     (::blas::Diag)diag, n, A_.ptr, A_.ldim, x_.ptr, x_.inc);
    else {
        conjugate(x);
        ::blas::trsv((::blas::Layout)L, (::blas::Uplo)uplo, ::blas::Op::NoTrans,
                     (::blas::Diag)diag, n, A_.ptr, A_.ldim, x_.ptr, x_.inc);
        conjugate(x);
    }
}
 
#endif
 
}  // namespace tlapack
 
#endif  //  #ifndef TLAPACK_BLAS_TRSV_HH