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