unmql.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_UNMQL_HH
#define TLAPACK_UNMQL_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/lapack/unmq.hpp"
 
namespace tlapack {
 
struct UnmqlOpts {
    size_t nb = 32;  
};
 
template <class T,
          TLAPACK_SMATRIX matrixA_t,
          TLAPACK_SMATRIX matrixC_t,
          TLAPACK_SVECTOR tau_t,
          TLAPACK_SIDE side_t,
          TLAPACK_OP trans_t>
constexpr WorkInfo unmql_worksize(side_t side,
                                  trans_t trans,
                                  const matrixA_t& A,
                                  const tau_t& tau,
                                  const matrixC_t& C,
                                  const UnmqlOpts& opts = {})
{
    using idx_t = size_type<matrixC_t>;
    using matrixT_t = matrix_type<matrixA_t, tau_t>;
    using range = pair<idx_t, idx_t>;
 
    // Constants
    const idx_t k = size(tau);
    const idx_t nb = min<idx_t>(opts.nb, k);
 
    // Local workspace sizes
    WorkInfo workinfo =
        (is_same_v<T, type_t<matrixT_t>>) ? WorkInfo(nb, nb) : WorkInfo(0);
 
    // larfb:
    {
        // Constants
        const idx_t m = nrows(C);
        const idx_t n = ncols(C);
        const idx_t nA = (side == Side::Left) ? m : n;
 
        // Empty matrices
        auto&& V = slice(A, range{0, nA}, range{0, nb});
        auto&& matrixT = slice(A, range{0, nb}, range{0, nb});
 
        // Internal workspace queries
        workinfo += larfb_worksize<T>(side, trans, BACKWARD, COLUMNWISE_STORAGE,
                                      V, matrixT, C);
    }
 
    return workinfo;
}
 
template <TLAPACK_SMATRIX matrixA_t,
          TLAPACK_SMATRIX matrixC_t,
          TLAPACK_SVECTOR tau_t,
          TLAPACK_SIDE side_t,
          TLAPACK_OP trans_t>
int unmql(side_t side,
          trans_t trans,
          const matrixA_t& A,
          const tau_t& tau,
          matrixC_t& C,
          const UnmqlOpts& opts = {})
{
    return unmq(side, trans, BACKWARD, COLUMNWISE_STORAGE, A, tau, C,
                UnmqOpts{opts.nb});
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_UNMQL_HH