her2.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_BLAS_HER2_HH
#define TLAPACK_BLAS_HER2_HH
 
#include "tlapack/base/utils.hpp"
 
namespace tlapack {
 
template <TLAPACK_MATRIX matrixA_t,
          TLAPACK_VECTOR vectorX_t,
          TLAPACK_VECTOR vectorY_t,
          TLAPACK_SCALAR alpha_t,
          class T = type_t<matrixA_t>,
          disable_if_allow_optblas_t<pair<alpha_t, T>,
                                     pair<matrixA_t, T>,
                                     pair<vectorX_t, T>,
                                     pair<vectorY_t, T> > = 0>
void her2(Uplo uplo,
          const alpha_t& alpha,
          const vectorX_t& x,
          const vectorY_t& y,
          matrixA_t& A)
{
    // data traits
    using idx_t = size_type<matrixA_t>;
    using TX = type_t<vectorX_t>;
    using TY = type_t<vectorY_t>;
 
    // constants
    const idx_t n = nrows(A);
 
    // check arguments
    tlapack_check_false(uplo != Uplo::Lower && uplo != Uplo::Upper);
    tlapack_check_false(size(x) != n);
    tlapack_check_false(size(y) != n);
    tlapack_check_false(ncols(A) != n);
 
    if (uplo == Uplo::Upper) {
        for (idx_t j = 0; j < n; ++j) {
            const scalar_type<alpha_t, TY> tmp1 = alpha * conj(y[j]);
            const scalar_type<alpha_t, TX> tmp2 = conj(alpha * x[j]);
            for (idx_t i = 0; i < j; ++i)
                A(i, j) += x[i] * tmp1 + y[i] * tmp2;
            A(j, j) = real(A(j, j)) + real(x[j]) * real(tmp1) -
                      imag(x[j]) * imag(tmp1) + real(y[j]) * real(tmp2) -
                      imag(y[j]) * imag(tmp2);
        }
    }
    else {
        for (idx_t j = 0; j < n; ++j) {
            const scalar_type<alpha_t, TY> tmp1 = alpha * conj(y[j]);
            const scalar_type<alpha_t, TX> tmp2 = conj(alpha * x[j]);
            A(j, j) = real(A(j, j)) + real(x[j]) * real(tmp1) -
                      imag(x[j]) * imag(tmp1) + real(y[j]) * real(tmp2) -
                      imag(y[j]) * imag(tmp2);
            for (idx_t i = j + 1; i < n; ++i)
                A(i, j) += x[i] * tmp1 + y[i] * tmp2;
        }
    }
}
 
#ifdef TLAPACK_USE_LAPACKPP
 
template <TLAPACK_LEGACY_MATRIX matrixA_t,
          TLAPACK_LEGACY_VECTOR vectorX_t,
          TLAPACK_LEGACY_VECTOR vectorY_t,
          TLAPACK_SCALAR alpha_t,
          class T = type_t<matrixA_t>,
          enable_if_allow_optblas_t<pair<alpha_t, T>,
                                    pair<matrixA_t, T>,
                                    pair<vectorX_t, T>,
                                    pair<vectorY_t, T> > = 0>
void her2(Uplo uplo,
          const alpha_t alpha,
          const vectorX_t& x,
          const vectorY_t& y,
          matrixA_t& A)
{
    // Legacy objects
    auto A_ = legacy_matrix(A);
    auto x_ = legacy_vector(x);
    auto y_ = legacy_vector(y);
 
    // Constants to forward
    constexpr Layout L = layout<matrixA_t>;
    const auto& n = A_.n;
 
    return ::blas::her2((::blas::Layout)L, (::blas::Uplo)uplo, n, alpha, x_.ptr,
                        x_.inc, y_.ptr, y_.inc, A_.ptr, A_.ldim);
}
 
#endif
 
}  // namespace tlapack
 
#endif  //  #ifndef TLAPACK_BLAS_HER2_HH