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_LEVEL2_GER_TCC
 34 #define CXXBLAS_LEVEL2_GER_TCC 1
 35 
 36 #include <complex>
 37 #include <cxxblas/level1/level1.h>
 38 
 39 namespace cxxblas {
 40 
 41 template <typename IndexType, typename ALPHA, typename VX, typename VY,
 42           typename MA>
 43 void
 44 gerc_generic(StorageOrder order, Transpose conjugateA,
 45              IndexType m, IndexType n,
 46              const ALPHA &alpha,
 47              const VX *x, IndexType incX,
 48              const VY *y, IndexType incY,
 49              MA *A, IndexType ldA)
 50 {
 51     if (order==ColMajor) {
 52         conjugateA = Transpose(conjugateA^Conj);
 53         gerc_generic(RowMajor, conjugateA, n, m,
 54                      alpha, y, incY, x, incX,
 55                      A, ldA);
 56         return;
 57     }
 58 
 59     #ifdef CXXBLAS_USE_XERBLA
 60         // insert error check here
 61     #endif
 62     if (conjugateA==Conj) {
 63         for (IndexType i=0, iX=0; i<m; ++i, iX+=incX) {
 64             axpy_generic(n, alpha*conjugate(x[iX]),
 65                             y, incY,
 66                             A+i*ldA, IndexType(1));
 67         }
 68     } else {
 69         for (IndexType i=0, iX=0; i<m; ++i, iX+=incX) {
 70             acxpy_generic(n, alpha*x[iX], y, incY, A+i*ldA, IndexType(1));
 71         }
 72     }
 73 }
 74 
 75 template <typename IndexType, typename ALPHA, typename VX, typename VY,
 76           typename MA>
 77 void
 78 geru_generic(StorageOrder order,
 79              IndexType m, IndexType n,
 80              const ALPHA &alpha,
 81              const VX *x, IndexType incX,
 82              const VY *y, IndexType incY,
 83              MA *A, IndexType ldA)
 84 {
 85     CXXBLAS_DEBUG_OUT("geru_generic");
 86     if (order==ColMajor) {
 87         geru_generic(RowMajor, n, m,
 88                      alpha, y, incY, x, incX,
 89                      A, ldA);
 90         return;
 91     }
 92     #ifdef CXXBLAS_USE_XERBLA
 93         // insert error check here
 94     #endif
 95     for (IndexType i=0, iX=0; i<m; ++i, iX+=incX) {
 96         axpy_generic(n, alpha*x[iX], y, incY, A+i*ldA, IndexType(1));
 97     }
 98 }
 99 
100 //------------------------------------------------------------------------------
101 
102 template <typename IndexType, typename ALPHA, typename VX, typename VY,
103           typename MA>
104 void
105 gerc(StorageOrder order,
106      IndexType m, IndexType n,
107      const ALPHA &alpha,
108      const VX *x, IndexType incX,
109      const VY *y, IndexType incY,
110      MA *A, IndexType ldA)
111 {
112     CXXBLAS_DEBUG_OUT("gerc_generic");
113 
114     if (incX<0) {
115         x -= incX*(m-1);
116     }
117     if (incY<0) {
118         y -= incY*(n-1);
119     }
120     gerc_generic(order, NoTrans, m, n, alpha, x, incX, y, incY, A, ldA);
121 }
122 
123 template <typename IndexType, typename ALPHA, typename VX, typename VY,
124           typename MA>
125 void
126 ger(StorageOrder order,
127     IndexType m, IndexType n,
128     const ALPHA &alpha,
129     const VX *x, IndexType incX,
130     const VY *y, IndexType incY,
131     MA *A, IndexType ldA)
132 {
133     gerc(order, m, n, alpha, x, incX, y, incY, A, ldA);
134 }
135 
136 template <typename IndexType, typename ALPHA, typename VX, typename VY,
137           typename MA>
138 void
139 geru(StorageOrder order,
140      IndexType m, IndexType n,
141      const ALPHA &alpha,
142      const VX *x, IndexType incX,
143      const VY *y, IndexType incY,
144      MA *A, IndexType ldA)
145 {
146     if (incX<0) {
147         x -= incX*(m-1);
148     }
149     if (incY<0) {
150         y -= incY*(n-1);
151     }
152     geru_generic(order, m, n, alpha, x, incX, y, incY, A, ldA);
153 }
154 
155 #ifdef HAVE_CBLAS
156 
157 // sger
158 template <typename IndexType>
159 typename If<IndexType>::isBlasCompatibleInteger
160 ger(StorageOrder order,
161     IndexType m, IndexType n,
162     const float &alpha,
163     const float *x, IndexType incX,
164     const float *y, IndexType incY,
165     float *A, IndexType ldA)
166 {
167     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_sger");
168 
169     cblas_sger(CBLAS::getCblasType(order),
170                m, n,
171                alpha,
172                x, incX,
173                y, incY,
174                A, ldA);
175 }
176 
177 // dger
178 template <typename IndexType>
179 typename If<IndexType>::isBlasCompatibleInteger
180 ger(StorageOrder order,
181     IndexType m, IndexType n,
182     const double &alpha,
183     const double *x, IndexType incX,
184     const double *y, IndexType incY,
185     double *A, IndexType ldA)
186 {
187     CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_dger");
188 
189     cblas_dger(CBLAS::getCblasType(order),
190                m, n,
191                alpha,
192                x, incX,
193                y, incY,
194                A, ldA);
195 }
196 
197 #endif // HAVE_CBLAS
198 
199 // namespace cxxblas
200 
201 #endif // CXXBLAS_LEVEL2_GER_TCC