gerqf.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_GERQF_HH
#define TLAPACK_GERQF_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/lapack/gerq2.hpp"
#include "tlapack/lapack/larfb.hpp"
#include "tlapack/lapack/larft.hpp"
 
namespace tlapack {
struct GerqfOpts {
    size_t nb = 32;  
};
 
template <class T, TLAPACK_SMATRIX A_t, TLAPACK_SVECTOR tau_t>
constexpr WorkInfo gerqf_worksize(const A_t& A,
                                  const tau_t& tau,
                                  const GerqfOpts& opts = {})
{
    using idx_t = size_type<A_t>;
    using range = pair<idx_t, idx_t>;
    using work_t = matrix_type<A_t, tau_t>;
 
    // constants
    const idx_t m = nrows(A);
    const idx_t n = ncols(A);
    const idx_t k = min(m, n);
    const idx_t nb = min((idx_t)opts.nb, k);
 
    auto&& A11 = rows(A, range(0, nb));
    auto&& tauw1 = slice(tau, range(0, nb));
    WorkInfo workinfo = gerq2_worksize<T>(A11, tauw1);
 
    if (m > nb) {
        auto&& TT1 = slice(A, range(0, nb), range(0, nb));
        auto&& A12 = slice(A, range(nb, m), range(0, n));
        workinfo.minMax(larfb_worksize<T>(RIGHT_SIDE, NO_TRANS, BACKWARD,
                                          ROWWISE_STORAGE, A11, TT1, A12));
        if constexpr (is_same_v<T, type_t<work_t>>)
            workinfo += WorkInfo(nb, nb);
    }
 
    return workinfo;
}
 
template <TLAPACK_SMATRIX A_t, TLAPACK_SVECTOR tau_t, TLAPACK_WORKSPACE work_t>
int gerqf_work(A_t& A, tau_t& tau, work_t& work, const GerqfOpts& opts = {})
{
    using idx_t = size_type<A_t>;
    using range = pair<idx_t, idx_t>;
 
    // constants
    const idx_t m = nrows(A);
    const idx_t n = ncols(A);
    const idx_t k = min(m, n);
    const idx_t nb = min((idx_t)opts.nb, k);
 
    // check arguments
    tlapack_check((idx_t)size(tau) >= k);
 
    // Matrix TT
    auto [TT, work2] = (m > nb) ? reshape(work, nb, nb) : reshape(work, 0, 0);
 
    // Main computational loop
    for (idx_t j2 = 0; j2 < k; j2 += nb) {
        const idx_t ib = min(nb, k - j2);
        const idx_t j = m - j2 - ib;
 
        // Compute the RQ factorization of the current block A(j:j+ib,0:n-j2)
        auto A11 = slice(A, range(j, j + ib), range(0, n - j2));
        auto tauw1 = slice(tau, range(k - j2 - ib, k - j2));
 
        gerq2_work(A11, tauw1, work);
 
        if (j > 0) {
            // Form the triangular factor of the block reflector
            auto TT1 = slice(TT, range(0, ib), range(0, ib));
            larft(BACKWARD, ROWWISE_STORAGE, A11, tauw1, TT1);
 
            // Apply H to A(0:j,0:n-j2) from the right
            auto A12 = slice(A, range(0, j), range(0, n - j2));
            larfb_work(RIGHT_SIDE, NO_TRANS, BACKWARD, ROWWISE_STORAGE, A11,
                       TT1, A12, work2);
        }
    }
 
    return 0;
}
 
template <TLAPACK_SMATRIX A_t, TLAPACK_SVECTOR tau_t>
int gerqf(A_t& A, tau_t& tau, const GerqfOpts& opts = {})
{
    using work_t = matrix_type<A_t, tau_t>;
    using T = type_t<work_t>;
    Create<work_t> new_matrix;
 
    // Allocate or get workspace
    WorkInfo workinfo = gerqf_worksize<T>(A, tau, opts);
    std::vector<T> work_;
    auto work = new_matrix(work_, workinfo.m, workinfo.n);
 
    return gerqf_work(A, tau, work, opts);
}
}  // namespace tlapack
#endif  // TLAPACK_GERQF_HH