gemv.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_LEGACY_GEMV_HH
#define TLAPACK_LEGACY_GEMV_HH
 
#include "tlapack/blas/gemv.hpp"
#include "tlapack/legacy_api/base/types.hpp"
#include "tlapack/legacy_api/base/utils.hpp"
 
namespace tlapack {
namespace legacy {
 
    template <typename TA, typename TX, typename TY>
    void gemv(Layout layout,
              Op trans,
              idx_t m,
              idx_t n,
              scalar_type<TA, TX, TY> alpha,
              TA const* A,
              idx_t lda,
              TX const* x,
              int_t incx,
              scalar_type<TA, TX, TY> beta,
              TY* y,
              int_t incy)
    {
        using internal::create_matrix;
        using scalar_t = scalar_type<TA, TX, TY>;
 
        // check arguments
        tlapack_check_false(layout != Layout::ColMajor &&
                            layout != Layout::RowMajor);
        tlapack_check_false(trans != Op::NoTrans && trans != Op::Trans &&
                            trans != Op::ConjTrans);
        tlapack_check_false(m < 0);
        tlapack_check_false(n < 0);
        tlapack_check_false(lda < ((layout == Layout::ColMajor) ? m : n));
        tlapack_check_false(incx == 0);
        tlapack_check_false(incy == 0);
 
        // quick return
        if (m == 0 || n == 0 ||
            ((alpha == scalar_t(0)) && (beta == scalar_t(1))))
            return;
 
        // Transpose if Row Major
        bool doConj = false;
        if (layout == Layout::RowMajor) {
            // A => A^T; A^T => A; A^H => A & doConj
            std::swap(m, n);
            if (trans == Op::NoTrans)
                trans = Op::Trans;
            else {
                if (trans == Op::ConjTrans) doConj = true;
                trans = Op::NoTrans;
            }
        }
 
        // Initialize indexes
        const idx_t lenx = ((trans == Op::NoTrans) ? n : m);
        const idx_t leny = ((trans == Op::NoTrans) ? m : n);
 
        // Conjugate if A is row-major and initially trans is Op::ConjTrans
        if (doConj) {
            TX* x2 = const_cast<TX*>(x);
            for (idx_t i = 0; i < lenx; ++i)
                x2[i * abs(incx)] = conj(x2[i * abs(incx)]);
            for (idx_t i = 0; i < leny; ++i)
                y[i * abs(incy)] = conj(y[i * abs(incy)]);
        }
 
        // Matrix views
        const auto A_ = create_matrix<TA>((TA*)A, m, n, lda);
 
        if (alpha == scalar_t(0)) {
            tlapack_expr_with_vector(
                y_, TY, leny, y, incy,
                if (beta == scalar_t(0)) for (idx_t i = 0; i < leny; ++i)
                    y_[i] = TY(0);
                else for (idx_t i = 0; i < leny; ++i) y_[i] *=
                (doConj) ? conj(beta) : beta);
        }
        else {
            tlapack_expr_with_2vectors(
                x_, TX, lenx, x, incx, y_, TY, leny, y, incy,
                if (beta == scalar_t(0))
                    gemv(trans, (doConj) ? conj(alpha) : alpha, A_, x_, y_);
                else gemv(trans, (doConj) ? conj(alpha) : alpha, A_, x_,
                          (doConj) ? conj(beta) : beta, y_));
        }
 
        // Conjugate if A is row-major and initially trans is Op::ConjTrans
        if (doConj) {
            TX* x2 = const_cast<TX*>(x);
            for (idx_t i = 0; i < lenx; ++i)
                x2[i * abs(incx)] = conj(x2[i * abs(incx)]);
            for (idx_t i = 0; i < leny; ++i)
                y[i * abs(incy)] = conj(y[i * abs(incy)]);
        }
    }
 
}  // namespace legacy
}  // namespace tlapack
 
#endif  //  #ifndef TLAPACK_LEGACY_GEMV_HH