1 #include <cxxblas/cxxblas.cxx>
  2 #include <cxxblas/interface/aux.h>
  3 #include <complex>
  4 
  5 using cxxblas::StorageOrder;
  6 using cxxblas::ColMajor;
  7 using cxxblas::RowMajor;
  8 using cxxblas::Transpose;
  9 using cxxblas::NoTrans;
 10 using cxxblas::Trans;
 11 using cxxblas::Conj;
 12 using cxxblas::ConjTrans;
 13 using cxxblas::Side;
 14 using cxxblas::Left;
 15 using cxxblas::Right;
 16 using cxxblas::Diag;
 17 using cxxblas::Unit;
 18 using cxxblas::NonUnit;
 19 
 20 
 21 extern "C" {
 22 
 23     void xerbla_(const char* srname, int* info);
 24 
 25     enum CBLAS_ORDER        {CblasRowMajor=101, CblasColMajor=102};
 26     enum CBLAS_TRANSPOSE    {CblasNoTrans=111, CblasTrans=112,
 27                              CblasConjTrans=113, CblasConjNoTrans=114};
 28     enum CBLAS_SIDE         {CblasLeft=141, CblasRight=142};
 29     enum CBLAS_UPLO         {CblasUpper=121, CblasLower=122};
 30     enum CBLAS_DIAG         {CblasNonUnit=131, CblasUnit=132};
 31 
 32 #ifndef COMPLEX_FLOAT1
 33     typedef CBLAS_FLOAT                 CBLAS_ALPHA;
 34     typedef CBLAS_FLOAT                 CXXBLAS_FLOAT;
 35 #else
 36     typedef const CBLAS_FLOAT *         CBLAS_ALPHA;
 37     typedef std::complex<CBLAS_FLOAT>   CXXBLAS_FLOAT;
 38 #endif
 39 
 40     void
 41     CBLAS_NAME(enum CBLAS_ORDER _order,
 42                enum CBLAS_UPLO _upLoC, enum CBLAS_TRANSPOSE _transAB,
 43                CBLAS_INT n, CBLAS_INT k,
 44                CBLAS_ALPHA _alpha,
 45                const CBLAS_FLOAT  *_A, CBLAS_INT ldA,
 46                const CBLAS_FLOAT  *_B, CBLAS_INT ldB,
 47                CBLAS_FLOAT _beta,
 48                CBLAS_FLOAT  *_C, CBLAS_INT ldC)
 49 #ifdef CREATE_CBLAS
 50     {
 51         StorageOrder order = (_order==CblasColMajor) ? ColMajor
 52                                                      : RowMajor;
 53 
 54         StorageUpLo upLoC = (_upLoC==CblasUpper) ? Upper : Lower;
 55 
 56         Transpose transAB = NoTrans;
 57         CBLAS_INT numRowsAB = n;
 58         CBLAS_INT numColsAB = k;
 59 
 60         if (_transAB==CblasTrans) {
 61             transAB = Trans;
 62             numRowsAB = k;
 63             numColsAB = n;
 64         }
 65         if (_transAB==CblasConjTrans) {
 66             transAB = ConjTrans;
 67             numRowsAB = k;
 68             numColsAB = n;
 69         }
 70         if (_transAB==CblasConjNoTrans) {
 71             transAB = Conj;
 72             numRowsAB = n;
 73             numColsAB = k;
 74         }
 75 
 76         const CXXBLAS_FLOAT *A = reinterpret_cast<const CXXBLAS_FLOAT *>(_A);
 77         const CXXBLAS_FLOAT *B = reinterpret_cast<const CXXBLAS_FLOAT *>(_B);
 78         CXXBLAS_FLOAT *C = reinterpret_cast<CXXBLAS_FLOAT *>(_C);
 79 
 80 #   ifndef COMPLEX_FLOAT1
 81         CXXBLAS_FLOAT alpha = _alpha;
 82 #   else
 83         CXXBLAS_FLOAT alpha(_alpha[0], _alpha[1]);
 84 #   endif
 85         CBLAS_FLOAT beta = _beta;
 86 
 87         CXXBLAS_FLOAT *__C = C;
 88         CBLAS_INT      __ldC = ldC;
 89         StorageOrder   __order = order;
 90 
 91 #   ifdef TEST_ROW_MAJOR
 92         switchFullStorageOrder(order, numRowsAB, numColsAB, A, ldA);
 93         switchFullStorageOrder(order, numRowsAB, numColsAB, B, ldB);
 94 
 95         __order = (order==ColMajor) ? RowMajor : ColMajor;
 96         allocateFullStorage(__order, n, n, __C, __ldC);
 97         switchFullStorageOrder(order, n, n, C, ldC, __C, __ldC);
 98 #   endif
 99 
100         cxxblas::her2k(__order, upLoC, transAB,
101                       n, k,
102                       alpha,
103                       A, ldA,
104                       B, ldB,
105                       beta,
106                       __C, __ldC);
107 
108 #   ifdef TEST_ROW_MAJOR
109         releaseStorage(A);
110         releaseStorage(B);
111 
112         switchFullStorageOrder(__order, n, n, __C, __ldC, C, ldC);
113         releaseStorage(__C);
114 #   endif
115     }
116 #else
117     ;
118 #endif // CREATE_CBLAS
119 
120 
121 #ifdef CREATE_BLAS
122     void
123     BLAS_NAME(const char *_upLoC, const char *_transAB,
124               const CBLAS_INT *_n, const CBLAS_INT *_k,
125               const CBLAS_FLOAT *_alpha,
126               const CBLAS_FLOAT *A, const CBLAS_INT *_ldA,
127               const CBLAS_FLOAT *B, const CBLAS_INT *_ldB,
128               const CBLAS_FLOAT *_beta,
129               CBLAS_FLOAT *C, const CBLAS_INT *_ldC)
130     {
131         CBLAS_UPLO          upLoC;
132         CBLAS_TRANSPOSE     transAB;
133         CBLAS_INT           numRowsAB  = -1;
134 
135         bool checkUpLoC = false;
136         bool checkTransAB = false;
137 
138         CBLAS_INT n       = *_n;
139         CBLAS_INT k       = *_k;
140         CBLAS_INT ldA     = *_ldA;
141         CBLAS_INT ldB     = *_ldB;
142         CBLAS_INT ldC     = *_ldC;
143 
144         if ((*_upLoC=='L') || (*_upLoC=='l')) {
145             upLoC = CblasLower;
146             checkUpLoC = true;
147         }
148         if ((*_upLoC=='U') || (*_upLoC=='u')) {
149             upLoC = CblasUpper;
150             checkUpLoC = true;
151         }
152 
153         if ((*_transAB=='N') || (*_transAB=='n')) {
154             transAB = CblasNoTrans;
155             numRowsAB = n;
156             checkTransAB = true;
157         }
158         if ((*_transAB=='T') || (*_transAB=='t')) {
159             transAB = CblasTrans;
160             numRowsAB = k;
161             checkTransAB = false;  // 'T' is illegal
162         }
163         if ((*_transAB=='C') || (*_transAB=='c')) {
164             transAB = CblasConjTrans;
165             numRowsAB = k;
166             checkTransAB = true;
167         }
168         if ((*_transAB=='R') || (*_transAB=='r')) {
169             transAB = CblasConjNoTrans;
170             numRowsAB = n;
171             checkTransAB = false// 'R' is illegal
172         }
173 
174 
175 #   ifndef COMPLEX_FLOAT1
176         CBLAS_ALPHA alpha  = *_alpha;
177 #   else
178         CBLAS_ALPHA alpha = _alpha;
179 #   endif
180         CBLAS_FLOAT beta  = *_beta;
181 
182         CBLAS_INT info = 0;
183         if (ldC<std::max(1, n)) {
184             info = 12;
185         }
186         if (ldB<std::max(1, numRowsAB)) {
187             info = 9;
188         }
189         if (ldA<std::max(1, numRowsAB)) {
190             info = 7;
191         }
192         if (k<0) {
193             info = 4;
194         }
195         if (n<0) {
196             info = 3;
197         }
198         if (!checkTransAB) {
199             info = 2;
200         }
201         if (!checkUpLoC) {
202             info = 1;
203         }
204         if (info!=0) {
205             char blasName[6];
206             strncpy(blasName, BLAS_NAME_STR, 6);
207             for (int i=0; i<6; ++i) {
208                 blasName[i] = std::toupper(blasName[i]);
209             }
210             xerbla_(blasName, &info);
211             return;
212         }
213 
214         // the blas interface calls the cblas interface
215         // so any blas-test will also test the cblas-interface
216         CBLAS_NAME(CblasColMajor, upLoC, transAB,
217                    n, k,
218                    alpha,
219                    A, ldA,
220                    B, ldB,
221                    beta,
222                    C, ldC);
223     }
224 #endif // CREATE_BLAS
225 
226 // extern "C"