1 /*
  2  *   Copyright (c) 2009, Michael Lehn
  3  *
  4  *   All rights reserved.
  5  *
  6  *   Redistribution and use in source and binary forms, with or without
  7  *   modification, are permitted provided that the following conditions
  8  *   are met:
  9  *
 10  *   1) Redistributions of source code must retain the above copyright
 11  *      notice, this list of conditions and the following disclaimer.
 12  *   2) Redistributions in binary form must reproduce the above copyright
 13  *      notice, this list of conditions and the following disclaimer in
 14  *      the documentation and/or other materials provided with the
 15  *      distribution.
 16  *   3) Neither the name of the FLENS development group nor the names of
 17  *      its contributors may be used to endorse or promote products derived
 18  *      from this software without specific prior written permission.
 19  *
 20  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 21  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 22  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 23  *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 24  *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 25  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 26  *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 27  *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 28  *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 29  *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 30  *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 31  */
 32 
 33 #ifndef CXXBLAS_LEVEL1_AXPY_TCC
 34 #define CXXBLAS_LEVEL1_AXPY_TCC 1
 35 
 36 #include <cstdio>
 37 #include <cxxblas/aux/aux.h>
 38 
 39 namespace cxxblas {
 40 
 41 template <typename IndexType, typename ALPHA, typename X, typename Y>
 42 void
 43 axpy_generic(IndexType n, const ALPHA &alpha, const X *x,
 44              IndexType incX, Y *y, IndexType incY)
 45 {
 46     CXXBLAS_DEBUG_OUT("axpy_generic");
 47 
 48     for (IndexType i=0, iX=0, iY=0; i<n; ++i, iX+=incX, iY+=incY) {
 49         y[iY] += alpha*x[iX];
 50     }
 51 }
 52 
 53 template <typename IndexType, typename ALPHA, typename X, typename Y>
 54 void
 55 axpy(IndexType n, const ALPHA &alpha, const X *x,
 56      IndexType incX, Y *y, IndexType incY)
 57 {
 58     if (incX<0) {
 59         x -= incX*(n-1);
 60     }
 61     if (incY<0) {
 62         y -= incY*(n-1);
 63     }
 64     axpy_generic(n, alpha, x, incX, y, incY);
 65 }
 66 
 67 #ifdef HAVE_CBLAS
 68 // saxpy
 69 template <typename IndexType>
 70 typename If<IndexType>::isBlasCompatibleInteger
 71 axpy(IndexType n, const float &alpha, const float *x, IndexType incX,
 72      float *y, IndexType incY)
 73 {
 74     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_saxpy");
 75 
 76     cblas_saxpy(n, alpha, x, incX, y, incY);
 77 }
 78 
 79 // daxpy
 80 template <typename IndexType>
 81 typename If<IndexType>::isBlasCompatibleInteger
 82 axpy(IndexType n, const double &alpha, const double *x, IndexType incX,
 83      double *y, IndexType incY)
 84 {
 85     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_daxpy");
 86 
 87     cblas_daxpy(n, alpha, x, incX, y, incY);
 88 }
 89 
 90 // caxpy
 91 template <typename IndexType>
 92 typename If<IndexType>::isBlasCompatibleInteger
 93 axpy(IndexType n, const ComplexFloat &alpha,
 94      const ComplexFloat *x, IndexType incX,
 95      ComplexFloat *y, IndexType incY)
 96 {
 97     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_caxpy");
 98 
 99     cblas_caxpy(n, reinterpret_cast<const float *>(&alpha),
100                    reinterpret_cast<const float *>(x), incX,
101                    reinterpret_cast<const float *>(y), incY);
102 }
103 
104 // zaxpy
105 template <typename IndexType>
106 typename If<IndexType>::isBlasCompatibleInteger
107 axpy(IndexType n, const ComplexDouble &alpha,
108      const ComplexDouble *x, IndexType incX,
109      ComplexDouble *y, IndexType incY)
110 {
111     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_zaxpy");
112 
113     cblas_zaxpy(n, reinterpret_cast<const double *>(&alpha),
114                    reinterpret_cast<const double *>(x), incX,
115                    reinterpret_cast<const double *>(y), incY);
116 }
117 
118 #endif // HAVE_CBLAS
119 
120 // namespace cxxblas
121 
122 #endif // CXXBLAS_LEVEL1_AXPY_TCC