rot.hh

// Copyright (c) 2017-2023, University of Tennessee. All rights reserved.
// SPDX-License-Identifier: BSD-3-Clause
// This program 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 BLAS_ROT_HH
#define BLAS_ROT_HH
 
#include "blas/util.hh"
 
#include <limits>
 
namespace blas {
 
// =============================================================================
 
template <typename TX, typename TY>
void rot(
    int64_t n,
    TX *x, int64_t incx,
    TY *y, int64_t incy,
    blas::real_type<TX, TY>   c,
    blas::scalar_type<TX, TY> s )
{
    typedef scalar_type<TX, TY> scalar_t;
 
    // check arguments
    blas_error_if( n < 0 );      // standard BLAS returns, doesn't fail
    blas_error_if( incx == 0 );
    blas_error_if( incy == 0 );
 
    scalar_t zero( 0 );
 
    // quick return
    if (n == 0 || (c == 1 && s == zero))
        return;
 
    if (incx == 1 && incy == 1) {
        // unit stride
        for (int64_t i = 0; i < n; ++i) {
            scalar_t stmp = c*x[i] + s*y[i];
            y[i] = c*y[i] - conj(s)*x[i];
            x[i] = stmp;
        }
    }
    else {
        // non-unit stride
        int64_t ix = (incx > 0 ? 0 : (-n + 1)*incx);
        int64_t iy = (incy > 0 ? 0 : (-n + 1)*incy);
        for (int64_t i = 0; i < n; ++i) {
            scalar_t stmp = c*x[ix] + s*y[iy];
            y[iy] = c*y[iy] - conj(s)*x[ix];
            x[ix] = stmp;
            ix += incx;
            iy += incy;
        }
    }
}
 
}  // namespace blas
 
#endif        //  #ifndef BLAS_ROT_HH