unmq.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_UNMQ_HH
#define TLAPACK_UNMQ_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/lapack/larfb.hpp"
#include "tlapack/lapack/larft.hpp"
 
namespace tlapack {
 
struct UnmqOpts {
    size_t nb = 32;  
};
 
template <class T,
          TLAPACK_SMATRIX matrixV_t,
          TLAPACK_SMATRIX matrixC_t,
          TLAPACK_SVECTOR vector_t,
          TLAPACK_SIDE side_t,
          TLAPACK_OP trans_t,
          TLAPACK_DIRECTION direction_t,
          TLAPACK_STOREV storage_t>
constexpr WorkInfo unmq_worksize(side_t side,
                                 trans_t trans,
                                 direction_t direction,
                                 storage_t storeMode,
                                 const matrixV_t& V,
                                 const vector_t& tau,
                                 const matrixC_t& C,
                                 const UnmqOpts& opts = {})
{
    using idx_t = size_type<matrixC_t>;
    using work_t = matrix_type<matrixV_t, vector_t>;
    using range = pair<idx_t, idx_t>;
 
    // Constants
    const idx_t m = nrows(C);
    const idx_t n = ncols(C);
    const idx_t k = size(tau);
    const idx_t nQ = (side == Side::Left) ? m : n;
    const idx_t nb = min((idx_t)opts.nb, k);
 
    // Local workspace sizes
    WorkInfo workinfo =
        (is_same_v<T, type_t<work_t>>) ? WorkInfo(nb, nb) : WorkInfo(0);
 
    auto&& Vi = (storeMode == StoreV::Columnwise)
                    ? slice(V, range{0, nQ}, range{0, nb})
                    : slice(V, range{0, nb}, range{0, nQ});
    auto&& matrixTi = slice(V, range{0, nb}, range{0, nb});
 
    // larfb:
    workinfo += larfb_worksize<T>(side, NO_TRANS, direction, storeMode, Vi,
                                  matrixTi, C);
 
    return workinfo;
}
 
template <TLAPACK_SMATRIX matrixV_t,
          TLAPACK_SMATRIX matrixC_t,
          TLAPACK_SVECTOR vector_t,
          TLAPACK_SIDE side_t,
          TLAPACK_OP trans_t,
          TLAPACK_DIRECTION direction_t,
          TLAPACK_STOREV storage_t,
          TLAPACK_WORKSPACE work_t>
int unmq_work(side_t side,
              trans_t trans,
              direction_t direction,
              storage_t storeMode,
              const matrixV_t& V,
              const vector_t& tau,
              matrixC_t& C,
              work_t& work,
              const UnmqOpts& opts = {})
{
    using idx_t = size_type<matrixC_t>;
    using range = pair<idx_t, idx_t>;
 
    // constants
    const idx_t m = nrows(C);
    const idx_t n = ncols(C);
    const idx_t k = size(tau);
    const idx_t nQ = (side == Side::Left) ? m : n;
    const idx_t nb = min((idx_t)opts.nb, k);
 
    // check arguments
    tlapack_check_false(side != Side::Left && side != Side::Right);
    tlapack_check_false(
        trans != Op::NoTrans && trans != Op::ConjTrans &&
        ((trans != Op::Trans) || is_complex<type_t<matrixV_t>>));
    tlapack_check_false(direction != Direction::Backward &&
                        direction != Direction::Forward);
    tlapack_check((storeMode == StoreV::Columnwise)
                      ? ((ncols(V) == k) && (nrows(V) == nQ))
                      : ((nrows(V) == k) && (ncols(V) == nQ)));
 
    // quick return
    if (m <= 0 || n <= 0 || k <= 0) return 0;
 
    // Matrix matrixT
    auto [matrixT, work1] = reshape(work, nb, nb);
 
    // const expressions
    const bool positiveIncLeft =
        (storeMode == StoreV::Columnwise)
            ? ((direction == Direction::Backward) ? (trans == Op::NoTrans)
                                                  : (trans != Op::NoTrans))
            : ((direction == Direction::Forward) ? (trans == Op::NoTrans)
                                                 : (trans != Op::NoTrans));
    const bool positiveInc =
        (side == Side::Left) ? positiveIncLeft : !positiveIncLeft;
    const idx_t i0 = (positiveInc) ? 0 : ((k - 1) / nb) * nb;
    const idx_t iN = (positiveInc) ? ((k - 1) / nb + 1) * nb : -nb;
    const idx_t inc = (positiveInc) ? nb : -nb;
 
    if (storeMode == StoreV::Columnwise) {
        for (idx_t i = i0; i != iN; i += inc) {
            const idx_t ib = min(nb, k - i);
            const auto rangev = (direction == Direction::Forward)
                                    ? range{i, nQ}
                                    : range{0, nQ - k + i + ib};
            const auto Vi = slice(V, rangev, range{i, i + ib});
            const auto taui = slice(tau, range{i, i + ib});
            auto matrixTi = slice(matrixT, range{0, ib}, range{0, ib});
 
            // Form the triangular factor of the block reflector
            larft(direction, COLUMNWISE_STORAGE, Vi, taui, matrixTi);
 
            // H or H**H is applied to either C[i:m,0:n] or C[0:m,i:n]
            auto Ci = (side == Side::Left) ? slice(C, rangev, range{0, n})
                                           : slice(C, range{0, m}, rangev);
            larfb_work(side, trans, direction, COLUMNWISE_STORAGE, Vi, matrixTi,
                       Ci, work1);
        }
    }
    else {
        for (idx_t i = i0; i != iN; i += inc) {
            const idx_t ib = min(nb, k - i);
            const auto rangev = (direction == Direction::Forward)
                                    ? range{i, nQ}
                                    : range{0, nQ - k + i + ib};
            const auto Vi = slice(V, range{i, i + ib}, rangev);
            const auto taui = slice(tau, range{i, i + ib});
            auto matrixTi = slice(matrixT, range{0, ib}, range{0, ib});
 
            // Form the triangular factor of the block reflector
            larft(direction, ROWWISE_STORAGE, Vi, taui, matrixTi);
 
            // H or H**H is applied to either C[i:m,0:n] or C[0:m,i:n]
            auto Ci = (side == Side::Left) ? slice(C, rangev, range{0, n})
                                           : slice(C, range{0, m}, rangev);
            larfb_work(side,
                       (trans == Op::NoTrans) ? Op::ConjTrans : Op::NoTrans,
                       direction, ROWWISE_STORAGE, Vi, matrixTi, Ci, work1);
        }
    }
 
    return 0;
}
 
template <TLAPACK_SMATRIX matrixV_t,
          TLAPACK_SMATRIX matrixC_t,
          TLAPACK_SVECTOR vector_t,
          TLAPACK_SIDE side_t,
          TLAPACK_OP trans_t,
          TLAPACK_DIRECTION direction_t,
          TLAPACK_STOREV storage_t>
int unmq(side_t side,
         trans_t trans,
         direction_t direction,
         storage_t storeMode,
         const matrixV_t& V,
         const vector_t& tau,
         matrixC_t& C,
         const UnmqOpts& opts = {})
{
    using idx_t = size_type<matrixC_t>;
    using work_t = matrix_type<matrixV_t, vector_t>;
    using T = type_t<work_t>;
 
    // Functor
    Create<work_t> new_matrix;
 
    // constants
    const idx_t m = nrows(C);
    const idx_t n = ncols(C);
    const idx_t k = size(tau);
 
    // quick return
    if (m <= 0 || n <= 0 || k <= 0) return 0;
 
    // Allocates workspace
    WorkInfo workinfo =
        unmq_worksize<T>(side, trans, direction, storeMode, V, tau, C, opts);
    std::vector<T> work_;
    auto work = new_matrix(work_, workinfo.m, workinfo.n);
 
    return unmq_work(side, trans, direction, storeMode, V, tau, C, work, opts);
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_UNMQ_HH