1/*
2 *     Written by D.P. Manley, Digital Equipment Corporation.
3 *     Prefixed "C_" to BLAS routines and their declarations.
4 *
5 *     Modified by T. H. Do, 2/19/98, SGI/CRAY Research.
6 */
7#include <stdlib.h>
8#include "cblas.h"
9#include "cblas_test.h"
10#define  TEST_COL_MJR	0
11#define  TEST_ROW_MJR	1
12#define  UNDEFINED     -1
13
14void F77_dgemm(int *order, char *transpa, char *transpb, int *m, int *n,
15              int *k, double *alpha, double *a, int *lda, double *b, int *ldb,
16              double *beta, double *c, int *ldc ) {
17
18  double *A, *B, *C;
19  int i,j,LDA, LDB, LDC;
20  enum CBLAS_TRANSPOSE transa, transb;
21
22  get_transpose_type(transpa, &transa);
23  get_transpose_type(transpb, &transb);
24
25  if (*order == TEST_ROW_MJR) {
26     if (transa == CblasNoTrans) {
27        LDA = *k+1;
28        A = (double *)malloc( (*m)*LDA*sizeof( double ) );
29        for( i=0; i<*m; i++ )
30           for( j=0; j<*k; j++ )
31              A[i*LDA+j]=a[j*(*lda)+i];
32     }
33     else {
34        LDA = *m+1;
35        A   = ( double* )malloc( LDA*(*k)*sizeof( double ) );
36        for( i=0; i<*k; i++ )
37           for( j=0; j<*m; j++ )
38              A[i*LDA+j]=a[j*(*lda)+i];
39     }
40     if (transb == CblasNoTrans) {
41        LDB = *n+1;
42        B   = ( double* )malloc( (*k)*LDB*sizeof( double ) );
43        for( i=0; i<*k; i++ )
44           for( j=0; j<*n; j++ )
45              B[i*LDB+j]=b[j*(*ldb)+i];
46     }
47     else {
48        LDB = *k+1;
49        B   = ( double* )malloc( LDB*(*n)*sizeof( double ) );
50        for( i=0; i<*n; i++ )
51           for( j=0; j<*k; j++ )
52              B[i*LDB+j]=b[j*(*ldb)+i];
53     }
54     LDC = *n+1;
55     C   = ( double* )malloc( (*m)*LDC*sizeof( double ) );
56     for( j=0; j<*n; j++ )
57        for( i=0; i<*m; i++ )
58           C[i*LDC+j]=c[j*(*ldc)+i];
59
60     cblas_dgemm( CblasRowMajor, transa, transb, *m, *n, *k, *alpha, A, LDA,
61                  B, LDB, *beta, C, LDC );
62     for( j=0; j<*n; j++ )
63        for( i=0; i<*m; i++ )
64           c[j*(*ldc)+i]=C[i*LDC+j];
65     free(A);
66     free(B);
67     free(C);
68  }
69  else if (*order == TEST_COL_MJR)
70     cblas_dgemm( CblasColMajor, transa, transb, *m, *n, *k, *alpha, a, *lda,
71                  b, *ldb, *beta, c, *ldc );
72  else
73     cblas_dgemm( UNDEFINED, transa, transb, *m, *n, *k, *alpha, a, *lda,
74                  b, *ldb, *beta, c, *ldc );
75}
76void F77_dsymm(int *order, char *rtlf, char *uplow, int *m, int *n,
77              double *alpha, double *a, int *lda, double *b, int *ldb,
78              double *beta, double *c, int *ldc ) {
79
80  double *A, *B, *C;
81  int i,j,LDA, LDB, LDC;
82  enum CBLAS_UPLO uplo;
83  enum CBLAS_SIDE side;
84
85  get_uplo_type(uplow,&uplo);
86  get_side_type(rtlf,&side);
87
88  if (*order == TEST_ROW_MJR) {
89     if (side == CblasLeft) {
90        LDA = *m+1;
91        A   = ( double* )malloc( (*m)*LDA*sizeof( double ) );
92        for( i=0; i<*m; i++ )
93           for( j=0; j<*m; j++ )
94              A[i*LDA+j]=a[j*(*lda)+i];
95     }
96     else{
97        LDA = *n+1;
98        A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
99        for( i=0; i<*n; i++ )
100           for( j=0; j<*n; j++ )
101              A[i*LDA+j]=a[j*(*lda)+i];
102     }
103     LDB = *n+1;
104     B   = ( double* )malloc( (*m)*LDB*sizeof( double ) );
105     for( i=0; i<*m; i++ )
106        for( j=0; j<*n; j++ )
107           B[i*LDB+j]=b[j*(*ldb)+i];
108     LDC = *n+1;
109     C   = ( double* )malloc( (*m)*LDC*sizeof( double ) );
110     for( j=0; j<*n; j++ )
111        for( i=0; i<*m; i++ )
112           C[i*LDC+j]=c[j*(*ldc)+i];
113     cblas_dsymm( CblasRowMajor, side, uplo, *m, *n, *alpha, A, LDA, B, LDB,
114                  *beta, C, LDC );
115     for( j=0; j<*n; j++ )
116        for( i=0; i<*m; i++ )
117           c[j*(*ldc)+i]=C[i*LDC+j];
118     free(A);
119     free(B);
120     free(C);
121  }
122  else if (*order == TEST_COL_MJR)
123     cblas_dsymm( CblasColMajor, side, uplo, *m, *n, *alpha, a, *lda, b, *ldb,
124                  *beta, c, *ldc );
125  else
126     cblas_dsymm( UNDEFINED, side, uplo, *m, *n, *alpha, a, *lda, b, *ldb,
127                  *beta, c, *ldc );
128}
129
130void F77_dsyrk(int *order, char *uplow, char *transp, int *n, int *k,
131              double *alpha, double *a, int *lda,
132              double *beta, double *c, int *ldc ) {
133
134  int i,j,LDA,LDC;
135  double *A, *C;
136  enum CBLAS_UPLO uplo;
137  enum CBLAS_TRANSPOSE trans;
138
139  get_uplo_type(uplow,&uplo);
140  get_transpose_type(transp,&trans);
141
142  if (*order == TEST_ROW_MJR) {
143     if (trans == CblasNoTrans) {
144        LDA = *k+1;
145        A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
146        for( i=0; i<*n; i++ )
147           for( j=0; j<*k; j++ )
148              A[i*LDA+j]=a[j*(*lda)+i];
149     }
150     else{
151        LDA = *n+1;
152        A   = ( double* )malloc( (*k)*LDA*sizeof( double ) );
153        for( i=0; i<*k; i++ )
154           for( j=0; j<*n; j++ )
155              A[i*LDA+j]=a[j*(*lda)+i];
156     }
157     LDC = *n+1;
158     C   = ( double* )malloc( (*n)*LDC*sizeof( double ) );
159     for( i=0; i<*n; i++ )
160        for( j=0; j<*n; j++ )
161           C[i*LDC+j]=c[j*(*ldc)+i];
162     cblas_dsyrk(CblasRowMajor, uplo, trans, *n, *k, *alpha, A, LDA, *beta,
163	         C, LDC );
164     for( j=0; j<*n; j++ )
165        for( i=0; i<*n; i++ )
166           c[j*(*ldc)+i]=C[i*LDC+j];
167     free(A);
168     free(C);
169  }
170  else if (*order == TEST_COL_MJR)
171     cblas_dsyrk(CblasColMajor, uplo, trans, *n, *k, *alpha, a, *lda, *beta,
172	         c, *ldc );
173  else
174     cblas_dsyrk(UNDEFINED, uplo, trans, *n, *k, *alpha, a, *lda, *beta,
175	         c, *ldc );
176}
177
178void F77_dsyr2k(int *order, char *uplow, char *transp, int *n, int *k,
179               double *alpha, double *a, int *lda, double *b, int *ldb,
180               double *beta, double *c, int *ldc ) {
181  int i,j,LDA,LDB,LDC;
182  double *A, *B, *C;
183  enum CBLAS_UPLO uplo;
184  enum CBLAS_TRANSPOSE trans;
185
186  get_uplo_type(uplow,&uplo);
187  get_transpose_type(transp,&trans);
188
189  if (*order == TEST_ROW_MJR) {
190     if (trans == CblasNoTrans) {
191        LDA = *k+1;
192        LDB = *k+1;
193        A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
194        B   = ( double* )malloc( (*n)*LDB*sizeof( double ) );
195        for( i=0; i<*n; i++ )
196           for( j=0; j<*k; j++ ) {
197              A[i*LDA+j]=a[j*(*lda)+i];
198              B[i*LDB+j]=b[j*(*ldb)+i];
199           }
200     }
201     else {
202        LDA = *n+1;
203        LDB = *n+1;
204        A   = ( double* )malloc( LDA*(*k)*sizeof( double ) );
205        B   = ( double* )malloc( LDB*(*k)*sizeof( double ) );
206        for( i=0; i<*k; i++ )
207           for( j=0; j<*n; j++ ){
208              A[i*LDA+j]=a[j*(*lda)+i];
209              B[i*LDB+j]=b[j*(*ldb)+i];
210           }
211     }
212     LDC = *n+1;
213     C   = ( double* )malloc( (*n)*LDC*sizeof( double ) );
214     for( i=0; i<*n; i++ )
215        for( j=0; j<*n; j++ )
216           C[i*LDC+j]=c[j*(*ldc)+i];
217     cblas_dsyr2k(CblasRowMajor, uplo, trans, *n, *k, *alpha, A, LDA,
218		  B, LDB, *beta, C, LDC );
219     for( j=0; j<*n; j++ )
220        for( i=0; i<*n; i++ )
221           c[j*(*ldc)+i]=C[i*LDC+j];
222     free(A);
223     free(B);
224     free(C);
225  }
226  else if (*order == TEST_COL_MJR)
227     cblas_dsyr2k(CblasColMajor, uplo, trans, *n, *k, *alpha, a, *lda,
228		   b, *ldb, *beta, c, *ldc );
229  else
230     cblas_dsyr2k(UNDEFINED, uplo, trans, *n, *k, *alpha, a, *lda,
231		   b, *ldb, *beta, c, *ldc );
232}
233void F77_dtrmm(int *order, char *rtlf, char *uplow, char *transp, char *diagn,
234              int *m, int *n, double *alpha, double *a, int *lda, double *b,
235              int *ldb) {
236  int i,j,LDA,LDB;
237  double *A, *B;
238  enum CBLAS_SIDE side;
239  enum CBLAS_DIAG diag;
240  enum CBLAS_UPLO uplo;
241  enum CBLAS_TRANSPOSE trans;
242
243  get_uplo_type(uplow,&uplo);
244  get_transpose_type(transp,&trans);
245  get_diag_type(diagn,&diag);
246  get_side_type(rtlf,&side);
247
248  if (*order == TEST_ROW_MJR) {
249     if (side == CblasLeft) {
250        LDA = *m+1;
251        A   = ( double* )malloc( (*m)*LDA*sizeof( double ) );
252        for( i=0; i<*m; i++ )
253           for( j=0; j<*m; j++ )
254              A[i*LDA+j]=a[j*(*lda)+i];
255     }
256     else{
257        LDA = *n+1;
258        A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
259        for( i=0; i<*n; i++ )
260           for( j=0; j<*n; j++ )
261              A[i*LDA+j]=a[j*(*lda)+i];
262     }
263     LDB = *n+1;
264     B   = ( double* )malloc( (*m)*LDB*sizeof( double ) );
265     for( i=0; i<*m; i++ )
266        for( j=0; j<*n; j++ )
267           B[i*LDB+j]=b[j*(*ldb)+i];
268     cblas_dtrmm(CblasRowMajor, side, uplo, trans, diag, *m, *n, *alpha,
269		 A, LDA, B, LDB );
270     for( j=0; j<*n; j++ )
271        for( i=0; i<*m; i++ )
272           b[j*(*ldb)+i]=B[i*LDB+j];
273     free(A);
274     free(B);
275  }
276  else if (*order == TEST_COL_MJR)
277     cblas_dtrmm(CblasColMajor, side, uplo, trans, diag, *m, *n, *alpha,
278		   a, *lda, b, *ldb);
279  else
280     cblas_dtrmm(UNDEFINED, side, uplo, trans, diag, *m, *n, *alpha,
281		   a, *lda, b, *ldb);
282}
283
284void F77_dtrsm(int *order, char *rtlf, char *uplow, char *transp, char *diagn,
285              int *m, int *n, double *alpha, double *a, int *lda, double *b,
286              int *ldb) {
287  int i,j,LDA,LDB;
288  double *A, *B;
289  enum CBLAS_SIDE side;
290  enum CBLAS_DIAG diag;
291  enum CBLAS_UPLO uplo;
292  enum CBLAS_TRANSPOSE trans;
293
294  get_uplo_type(uplow,&uplo);
295  get_transpose_type(transp,&trans);
296  get_diag_type(diagn,&diag);
297  get_side_type(rtlf,&side);
298
299  if (*order == TEST_ROW_MJR) {
300     if (side == CblasLeft) {
301        LDA = *m+1;
302        A   = ( double* )malloc( (*m)*LDA*sizeof( double ) );
303        for( i=0; i<*m; i++ )
304           for( j=0; j<*m; j++ )
305              A[i*LDA+j]=a[j*(*lda)+i];
306     }
307     else{
308        LDA = *n+1;
309        A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
310        for( i=0; i<*n; i++ )
311           for( j=0; j<*n; j++ )
312              A[i*LDA+j]=a[j*(*lda)+i];
313     }
314     LDB = *n+1;
315     B   = ( double* )malloc( (*m)*LDB*sizeof( double ) );
316     for( i=0; i<*m; i++ )
317        for( j=0; j<*n; j++ )
318           B[i*LDB+j]=b[j*(*ldb)+i];
319     cblas_dtrsm(CblasRowMajor, side, uplo, trans, diag, *m, *n, *alpha,
320		 A, LDA, B, LDB );
321     for( j=0; j<*n; j++ )
322        for( i=0; i<*m; i++ )
323           b[j*(*ldb)+i]=B[i*LDB+j];
324     free(A);
325     free(B);
326  }
327  else if (*order == TEST_COL_MJR)
328     cblas_dtrsm(CblasColMajor, side, uplo, trans, diag, *m, *n, *alpha,
329		   a, *lda, b, *ldb);
330  else
331     cblas_dtrsm(UNDEFINED, side, uplo, trans, diag, *m, *n, *alpha,
332		   a, *lda, b, *ldb);
333}
334