larft.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_LARFT_HH
#define TLAPACK_LARFT_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/blas/gemm.hpp"
#include "tlapack/blas/gemv.hpp"
#include "tlapack/blas/trmv.hpp"
 
namespace tlapack {
 
template <TLAPACK_DIRECTION direction_t,
          TLAPACK_STOREV storage_t,
          TLAPACK_SMATRIX matrixV_t,
          TLAPACK_VECTOR vector_t,
          TLAPACK_SMATRIX matrixT_t>
int larft(direction_t direction,
          storage_t storeMode,
          const matrixV_t& V,
          const vector_t& tau,
          matrixT_t& T)
{
    // data traits
    using scalar_t = type_t<matrixT_t>;
    using real_t = real_type<scalar_t>;
    using tau_t = type_t<vector_t>;
    using idx_t = size_type<matrixV_t>;
 
    // using
    using range = pair<idx_t, idx_t>;
 
    // constants
    const real_t one(1);
    const real_t zero(0);
    const idx_t n = (storeMode == StoreV::Columnwise) ? nrows(V) : ncols(V);
    const idx_t k = size(tau);
 
    // check arguments
    tlapack_check_false(direction != Direction::Backward &&
                        direction != Direction::Forward);
    tlapack_check_false(storeMode != StoreV::Columnwise &&
                        storeMode != StoreV::Rowwise);
    tlapack_check_false(
        k > ((storeMode == StoreV::Columnwise) ? ncols(V) : nrows(V)));
    tlapack_check_false(nrows(T) < k || ncols(T) < k);
 
    // Quick return
    if (n == 0 || k == 0) return 0;
 
    if (direction == Direction::Forward) {
        // First iteration:
        T(0, 0) = tau[0];
 
        // Remaining iterations:
        for (idx_t i = 1; i < k; ++i) {
            // Column vector t := T(0:i,i)
            auto t = slice(T, range{0, i}, i);
 
            if (tau[i] == tau_t(0)) {
                // H(i) =  I
                for (idx_t j = 0; j < i; ++j)
                    t[j] = zero;
            }
 
            else {
                // General case
                if (storeMode == StoreV::Columnwise) {
                    // t := - tau[i] conj( V(i,0:i) )
                    for (idx_t j = 0; j < i; ++j)
                        t[j] = -tau[i] * conj(V(i, j));
 
                    // t := t - tau[i] V(i+1:n,0:i)^H V(i+1:n,i)
                    if (i + 1 < n) {
                        gemv(CONJ_TRANS, -tau[i],
                             slice(V, range{i + 1, n}, range{0, i}),
                             slice(V, range{i + 1, n}, i), one, t);
                    }
                }
                else {
                    // t := - tau[i] V(0:i,i)
                    for (idx_t j = 0; j < i; ++j)
                        t[j] = -tau[i] * V(j, i);
 
                    // t := t - tau[i] V(0:i,i:n) V(i,i+1:n)^H
                    if (i + 1 < n) {
                        auto Ti = slice(T, range{0, i}, range{i, i + 1});
                        gemm(NO_TRANS, CONJ_TRANS, -tau[i],
                             slice(V, range{0, i}, range{i + 1, n}),
                             slice(V, range{i, i + 1}, range{i + 1, n}), one,
                             Ti);
                    }
                }
 
                // t := T(0:i,0:i) * t
                trmv(UPPER_TRIANGLE, NO_TRANS, NON_UNIT_DIAG,
                     slice(T, range{0, i}, range{0, i}), t);
            }
 
            // Update diagonal
            T(i, i) = tau[i];
        }
    }
    else  // direct==Direction::Backward
    {
        // First iteration:
        T(k - 1, k - 1) = tau[k - 1];
 
        // Remaining iterations:
        for (idx_t i = k - 2; i != idx_t(-1); --i) {
            // Column vector t := T(0:i,i)
            auto t = slice(T, range{i + 1, k}, i);
 
            if (tau[i] == tau_t(0)) {
                // H(i) =  I
                for (idx_t j = 0; j < k - i - 1; ++j)
                    t[j] = zero;
            }
 
            else {
                // General case
                if (storeMode == StoreV::Columnwise) {
                    // t := - tau[i] conj(V(n-k+i,i+1:k))
                    for (idx_t j = 0; j < k - i - 1; ++j)
                        t[j] = -tau[i] * conj(V(n - k + i, j + i + 1));
 
                    // t := t - tau[i] V(0:n-k+i,i+1:k)^H V(0:n-k+i,i)
                    gemv(CONJ_TRANS, -tau[i],
                         slice(V, range{0, n - k + i}, range{i + 1, k}),
                         slice(V, range{0, n - k + i}, i), one, t);
                }
                else {
                    // t := - tau[i] V(i+1:k,n-k+i)
                    for (idx_t j = 0; j < k - i - 1; ++j)
                        t[j] = -tau[i] * V(j + i + 1, n - k + i);
 
                    // t := t - tau[i] V(i+1:k,0:n-k+i) V(i,0:n-k+i)^H
                    auto Ti = slice(T, range{i + 1, k}, range{i, i + 1});
                    gemm(NO_TRANS, CONJ_TRANS, -tau[i],
                         slice(V, range{i + 1, k}, range{0, n - k + i}),
                         slice(V, range{i, i + 1}, range{0, n - k + i}), one,
                         Ti);
                }
 
                // t := T(i+1:k,i+1:k) * t
                trmv(LOWER_TRIANGLE, NO_TRANS, NON_UNIT_DIAG,
                     slice(T, range{i + 1, k}, range{i + 1, k}), t);
            }
 
            // Update diagonal
            T(i, i) = tau[i];
        }
    }
    return 0;
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_LARFT_HH