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, 1/23/98, SGI/CRAY Research.
6 */
7#include <stdlib.h>
8#include "cblas.h"
9#include "cblas_test.h"
10
11void F77_dgemv(int *order, char *transp, int *m, int *n, double *alpha,
12	       double *a, int *lda, double *x, int *incx, double *beta,
13	       double *y, int *incy ) {
14
15  double *A;
16  int i,j,LDA;
17  enum CBLAS_TRANSPOSE trans;
18
19  get_transpose_type(transp, &trans);
20  if (*order == TEST_ROW_MJR) {
21     LDA = *n+1;
22     A   = ( double* )malloc( (*m)*LDA*sizeof( double ) );
23     for( i=0; i<*m; i++ )
24        for( j=0; j<*n; j++ )
25           A[ LDA*i+j ]=a[ (*lda)*j+i ];
26     cblas_dgemv( CblasRowMajor, trans,
27		  *m, *n, *alpha, A, LDA, x, *incx, *beta, y, *incy );
28     free(A);
29  }
30  else if (*order == TEST_COL_MJR)
31     cblas_dgemv( CblasColMajor, trans,
32		  *m, *n, *alpha, a, *lda, x, *incx, *beta, y, *incy );
33  else
34     cblas_dgemv( UNDEFINED, trans,
35		  *m, *n, *alpha, a, *lda, x, *incx, *beta, y, *incy );
36}
37
38void F77_dger(int *order, int *m, int *n, double *alpha, double *x, int *incx,
39	     double *y, int *incy, double *a, int *lda ) {
40
41  double *A;
42  int i,j,LDA;
43
44  if (*order == TEST_ROW_MJR) {
45     LDA = *n+1;
46     A   = ( double* )malloc( (*m)*LDA*sizeof( double ) );
47
48     for( i=0; i<*m; i++ ) {
49       for( j=0; j<*n; j++ )
50         A[ LDA*i+j ]=a[ (*lda)*j+i ];
51     }
52
53     cblas_dger(CblasRowMajor, *m, *n, *alpha, x, *incx, y, *incy, A, LDA );
54     for( i=0; i<*m; i++ )
55       for( j=0; j<*n; j++ )
56         a[ (*lda)*j+i ]=A[ LDA*i+j ];
57     free(A);
58  }
59  else
60     cblas_dger( CblasColMajor, *m, *n, *alpha, x, *incx, y, *incy, a, *lda );
61}
62
63void F77_dtrmv(int *order, char *uplow, char *transp, char *diagn,
64	      int *n, double *a, int *lda, double *x, int *incx) {
65  double *A;
66  int i,j,LDA;
67  enum CBLAS_TRANSPOSE trans;
68  enum CBLAS_UPLO uplo;
69  enum CBLAS_DIAG diag;
70
71  get_transpose_type(transp,&trans);
72  get_uplo_type(uplow,&uplo);
73  get_diag_type(diagn,&diag);
74
75  if (*order == TEST_ROW_MJR) {
76     LDA = *n+1;
77     A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
78     for( i=0; i<*n; i++ )
79       for( j=0; j<*n; j++ )
80         A[ LDA*i+j ]=a[ (*lda)*j+i ];
81     cblas_dtrmv(CblasRowMajor, uplo, trans, diag, *n, A, LDA, x, *incx);
82     free(A);
83  }
84  else if (*order == TEST_COL_MJR)
85     cblas_dtrmv(CblasColMajor, uplo, trans, diag, *n, a, *lda, x, *incx);
86  else {
87     cblas_dtrmv(UNDEFINED, uplo, trans, diag, *n, a, *lda, x, *incx);
88  }
89}
90
91void F77_dtrsv(int *order, char *uplow, char *transp, char *diagn,
92	       int *n, double *a, int *lda, double *x, int *incx ) {
93  double *A;
94  int i,j,LDA;
95  enum CBLAS_TRANSPOSE trans;
96  enum CBLAS_UPLO uplo;
97  enum CBLAS_DIAG diag;
98
99  get_transpose_type(transp,&trans);
100  get_uplo_type(uplow,&uplo);
101  get_diag_type(diagn,&diag);
102
103  if (*order == TEST_ROW_MJR) {
104     LDA = *n+1;
105     A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
106     for( i=0; i<*n; i++ )
107        for( j=0; j<*n; j++ )
108           A[ LDA*i+j ]=a[ (*lda)*j+i ];
109     cblas_dtrsv(CblasRowMajor, uplo, trans, diag, *n, A, LDA, x, *incx );
110     free(A);
111   }
112   else
113     cblas_dtrsv(CblasColMajor, uplo, trans, diag, *n, a, *lda, x, *incx );
114}
115void F77_dsymv(int *order, char *uplow, int *n, double *alpha, double *a,
116	      int *lda, double *x, int *incx, double *beta, double *y,
117	      int *incy) {
118  double *A;
119  int i,j,LDA;
120  enum CBLAS_UPLO uplo;
121
122  get_uplo_type(uplow,&uplo);
123
124  if (*order == TEST_ROW_MJR) {
125     LDA = *n+1;
126     A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
127     for( i=0; i<*n; i++ )
128        for( j=0; j<*n; j++ )
129           A[ LDA*i+j ]=a[ (*lda)*j+i ];
130     cblas_dsymv(CblasRowMajor, uplo, *n, *alpha, A, LDA, x, *incx,
131		 *beta, y, *incy );
132     free(A);
133   }
134   else
135     cblas_dsymv(CblasColMajor, uplo, *n, *alpha, a, *lda, x, *incx,
136		 *beta, y, *incy );
137}
138
139void F77_dsyr(int *order, char *uplow, int *n, double *alpha, double *x,
140	     int *incx, double *a, int *lda) {
141  double *A;
142  int i,j,LDA;
143  enum CBLAS_UPLO uplo;
144
145  get_uplo_type(uplow,&uplo);
146
147  if (*order == TEST_ROW_MJR) {
148     LDA = *n+1;
149     A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
150     for( i=0; i<*n; i++ )
151        for( j=0; j<*n; j++ )
152           A[ LDA*i+j ]=a[ (*lda)*j+i ];
153     cblas_dsyr(CblasRowMajor, uplo, *n, *alpha, x, *incx, A, LDA);
154     for( i=0; i<*n; i++ )
155       for( j=0; j<*n; j++ )
156         a[ (*lda)*j+i ]=A[ LDA*i+j ];
157     free(A);
158   }
159   else
160     cblas_dsyr(CblasColMajor, uplo, *n, *alpha, x, *incx, a, *lda);
161}
162
163void F77_dsyr2(int *order, char *uplow, int *n, double *alpha, double *x,
164	     int *incx, double *y, int *incy, double *a, int *lda) {
165  double *A;
166  int i,j,LDA;
167  enum CBLAS_UPLO uplo;
168
169  get_uplo_type(uplow,&uplo);
170
171  if (*order == TEST_ROW_MJR) {
172     LDA = *n+1;
173     A   = ( double* )malloc( (*n)*LDA*sizeof( double ) );
174     for( i=0; i<*n; i++ )
175        for( j=0; j<*n; j++ )
176           A[ LDA*i+j ]=a[ (*lda)*j+i ];
177     cblas_dsyr2(CblasRowMajor, uplo, *n, *alpha, x, *incx, y, *incy, A, LDA);
178     for( i=0; i<*n; i++ )
179       for( j=0; j<*n; j++ )
180         a[ (*lda)*j+i ]=A[ LDA*i+j ];
181     free(A);
182   }
183   else
184     cblas_dsyr2(CblasColMajor, uplo, *n, *alpha, x, *incx, y, *incy, a, *lda);
185}
186
187void F77_dgbmv(int *order, char *transp, int *m, int *n, int *kl, int *ku,
188	       double *alpha, double *a, int *lda, double *x, int *incx,
189	       double *beta, double *y, int *incy ) {
190
191  double *A;
192  int i,irow,j,jcol,LDA;
193  enum CBLAS_TRANSPOSE trans;
194
195  get_transpose_type(transp, &trans);
196
197  if (*order == TEST_ROW_MJR) {
198     LDA = *ku+*kl+2;
199     A   = ( double* )malloc( (*n+*kl)*LDA*sizeof( double ) );
200     for( i=0; i<*ku; i++ ){
201        irow=*ku+*kl-i;
202        jcol=(*ku)-i;
203        for( j=jcol; j<*n; j++ )
204           A[ LDA*(j-jcol)+irow ]=a[ (*lda)*j+i ];
205     }
206     i=*ku;
207     irow=*ku+*kl-i;
208     for( j=0; j<*n; j++ )
209        A[ LDA*j+irow ]=a[ (*lda)*j+i ];
210     for( i=*ku+1; i<*ku+*kl+1; i++ ){
211        irow=*ku+*kl-i;
212        jcol=i-(*ku);
213        for( j=jcol; j<(*n+*kl); j++ )
214           A[ LDA*j+irow ]=a[ (*lda)*(j-jcol)+i ];
215     }
216     cblas_dgbmv( CblasRowMajor, trans, *m, *n, *kl, *ku, *alpha,
217		  A, LDA, x, *incx, *beta, y, *incy );
218     free(A);
219  }
220  else
221     cblas_dgbmv( CblasColMajor, trans, *m, *n, *kl, *ku, *alpha,
222		  a, *lda, x, *incx, *beta, y, *incy );
223}
224
225void F77_dtbmv(int *order, char *uplow, char *transp, char *diagn,
226	      int *n, int *k, double *a, int *lda, double *x, int *incx) {
227  double *A;
228  int irow, jcol, i, j, LDA;
229  enum CBLAS_TRANSPOSE trans;
230  enum CBLAS_UPLO uplo;
231  enum CBLAS_DIAG diag;
232
233  get_transpose_type(transp,&trans);
234  get_uplo_type(uplow,&uplo);
235  get_diag_type(diagn,&diag);
236
237  if (*order == TEST_ROW_MJR) {
238     LDA = *k+1;
239     A = ( double* )malloc( (*n+*k)*LDA*sizeof( double ) );
240     if (uplo == CblasUpper) {
241        for( i=0; i<*k; i++ ){
242           irow=*k-i;
243           jcol=(*k)-i;
244           for( j=jcol; j<*n; j++ )
245              A[ LDA*(j-jcol)+irow ]=a[ (*lda)*j+i ];
246        }
247        i=*k;
248        irow=*k-i;
249        for( j=0; j<*n; j++ )
250           A[ LDA*j+irow ]=a[ (*lda)*j+i ];
251     }
252     else {
253       i=0;
254       irow=*k-i;
255       for( j=0; j<*n; j++ )
256          A[ LDA*j+irow ]=a[ (*lda)*j+i ];
257       for( i=1; i<*k+1; i++ ){
258          irow=*k-i;
259          jcol=i;
260          for( j=jcol; j<(*n+*k); j++ )
261             A[ LDA*j+irow ]=a[ (*lda)*(j-jcol)+i ];
262       }
263     }
264     cblas_dtbmv(CblasRowMajor, uplo, trans, diag, *n, *k, A, LDA, x, *incx);
265     free(A);
266   }
267   else
268     cblas_dtbmv(CblasColMajor, uplo, trans, diag, *n, *k, a, *lda, x, *incx);
269}
270
271void F77_dtbsv(int *order, char *uplow, char *transp, char *diagn,
272	      int *n, int *k, double *a, int *lda, double *x, int *incx) {
273  double *A;
274  int irow, jcol, i, j, LDA;
275  enum CBLAS_TRANSPOSE trans;
276  enum CBLAS_UPLO uplo;
277  enum CBLAS_DIAG diag;
278
279  get_transpose_type(transp,&trans);
280  get_uplo_type(uplow,&uplo);
281  get_diag_type(diagn,&diag);
282
283  if (*order == TEST_ROW_MJR) {
284     LDA = *k+1;
285     A = ( double* )malloc( (*n+*k)*LDA*sizeof( double ) );
286     if (uplo == CblasUpper) {
287        for( i=0; i<*k; i++ ){
288        irow=*k-i;
289        jcol=(*k)-i;
290        for( j=jcol; j<*n; j++ )
291           A[ LDA*(j-jcol)+irow ]=a[ (*lda)*j+i ];
292        }
293        i=*k;
294        irow=*k-i;
295        for( j=0; j<*n; j++ )
296           A[ LDA*j+irow ]=a[ (*lda)*j+i ];
297     }
298     else {
299        i=0;
300        irow=*k-i;
301        for( j=0; j<*n; j++ )
302           A[ LDA*j+irow ]=a[ (*lda)*j+i ];
303        for( i=1; i<*k+1; i++ ){
304           irow=*k-i;
305           jcol=i;
306           for( j=jcol; j<(*n+*k); j++ )
307              A[ LDA*j+irow ]=a[ (*lda)*(j-jcol)+i ];
308        }
309     }
310     cblas_dtbsv(CblasRowMajor, uplo, trans, diag, *n, *k, A, LDA, x, *incx);
311     free(A);
312  }
313  else
314     cblas_dtbsv(CblasColMajor, uplo, trans, diag, *n, *k, a, *lda, x, *incx);
315}
316
317void F77_dsbmv(int *order, char *uplow, int *n, int *k, double *alpha,
318	      double *a, int *lda, double *x, int *incx, double *beta,
319	      double *y, int *incy) {
320  double *A;
321  int i,j,irow,jcol,LDA;
322  enum CBLAS_UPLO uplo;
323
324  get_uplo_type(uplow,&uplo);
325
326  if (*order == TEST_ROW_MJR) {
327     LDA = *k+1;
328     A   = ( double* )malloc( (*n+*k)*LDA*sizeof( double ) );
329     if (uplo == CblasUpper) {
330        for( i=0; i<*k; i++ ){
331           irow=*k-i;
332           jcol=(*k)-i;
333           for( j=jcol; j<*n; j++ )
334        A[ LDA*(j-jcol)+irow ]=a[ (*lda)*j+i ];
335        }
336        i=*k;
337        irow=*k-i;
338        for( j=0; j<*n; j++ )
339           A[ LDA*j+irow ]=a[ (*lda)*j+i ];
340     }
341     else {
342        i=0;
343        irow=*k-i;
344        for( j=0; j<*n; j++ )
345           A[ LDA*j+irow ]=a[ (*lda)*j+i ];
346        for( i=1; i<*k+1; i++ ){
347           irow=*k-i;
348           jcol=i;
349           for( j=jcol; j<(*n+*k); j++ )
350              A[ LDA*j+irow ]=a[ (*lda)*(j-jcol)+i ];
351        }
352     }
353     cblas_dsbmv(CblasRowMajor, uplo, *n, *k, *alpha, A, LDA, x, *incx,
354		 *beta, y, *incy );
355     free(A);
356   }
357   else
358     cblas_dsbmv(CblasColMajor, uplo, *n, *k, *alpha, a, *lda, x, *incx,
359		 *beta, y, *incy );
360}
361
362void F77_dspmv(int *order, char *uplow, int *n, double *alpha, double *ap,
363	      double *x, int *incx, double *beta, double *y, int *incy) {
364  double *A,*AP;
365  int i,j,k,LDA;
366  enum CBLAS_UPLO uplo;
367
368  get_uplo_type(uplow,&uplo);
369
370  if (*order == TEST_ROW_MJR) {
371     LDA = *n;
372     A   = ( double* )malloc( LDA*LDA*sizeof( double ) );
373     AP  = ( double* )malloc( (((LDA+1)*LDA)/2)*sizeof( double ) );
374     if (uplo == CblasUpper) {
375        for( j=0, k=0; j<*n; j++ )
376           for( i=0; i<j+1; i++, k++ )
377              A[ LDA*i+j ]=ap[ k ];
378        for( i=0, k=0; i<*n; i++ )
379           for( j=i; j<*n; j++, k++ )
380              AP[ k ]=A[ LDA*i+j ];
381     }
382     else {
383        for( j=0, k=0; j<*n; j++ )
384           for( i=j; i<*n; i++, k++ )
385              A[ LDA*i+j ]=ap[ k ];
386        for( i=0, k=0; i<*n; i++ )
387           for( j=0; j<i+1; j++, k++ )
388              AP[ k ]=A[ LDA*i+j ];
389     }
390     cblas_dspmv( CblasRowMajor, uplo, *n, *alpha, AP, x, *incx, *beta, y,
391		  *incy );
392     free(A);
393     free(AP);
394  }
395  else
396     cblas_dspmv( CblasColMajor, uplo, *n, *alpha, ap, x, *incx, *beta, y,
397		  *incy );
398}
399
400void F77_dtpmv(int *order, char *uplow, char *transp, char *diagn,
401	      int *n, double *ap, double *x, int *incx) {
402  double *A, *AP;
403  int i, j, k, LDA;
404  enum CBLAS_TRANSPOSE trans;
405  enum CBLAS_UPLO uplo;
406  enum CBLAS_DIAG diag;
407
408  get_transpose_type(transp,&trans);
409  get_uplo_type(uplow,&uplo);
410  get_diag_type(diagn,&diag);
411
412  if (*order == TEST_ROW_MJR) {
413     LDA = *n;
414     A   = ( double* )malloc( LDA*LDA*sizeof( double ) );
415     AP  = ( double* )malloc( (((LDA+1)*LDA)/2)*sizeof( double ) );
416     if (uplo == CblasUpper) {
417        for( j=0, k=0; j<*n; j++ )
418           for( i=0; i<j+1; i++, k++ )
419              A[ LDA*i+j ]=ap[ k ];
420        for( i=0, k=0; i<*n; i++ )
421           for( j=i; j<*n; j++, k++ )
422              AP[ k ]=A[ LDA*i+j ];
423     }
424     else {
425        for( j=0, k=0; j<*n; j++ )
426           for( i=j; i<*n; i++, k++ )
427              A[ LDA*i+j ]=ap[ k ];
428        for( i=0, k=0; i<*n; i++ )
429           for( j=0; j<i+1; j++, k++ )
430              AP[ k ]=A[ LDA*i+j ];
431     }
432     cblas_dtpmv( CblasRowMajor, uplo, trans, diag, *n, AP, x, *incx );
433     free(A);
434     free(AP);
435  }
436  else
437     cblas_dtpmv( CblasColMajor, uplo, trans, diag, *n, ap, x, *incx );
438}
439
440void F77_dtpsv(int *order, char *uplow, char *transp, char *diagn,
441	      int *n, double *ap, double *x, int *incx) {
442  double *A, *AP;
443  int i, j, k, LDA;
444  enum CBLAS_TRANSPOSE trans;
445  enum CBLAS_UPLO uplo;
446  enum CBLAS_DIAG diag;
447
448  get_transpose_type(transp,&trans);
449  get_uplo_type(uplow,&uplo);
450  get_diag_type(diagn,&diag);
451
452  if (*order == TEST_ROW_MJR) {
453     LDA = *n;
454     A   = ( double* )malloc( LDA*LDA*sizeof( double ) );
455     AP  = ( double* )malloc( (((LDA+1)*LDA)/2)*sizeof( double ) );
456     if (uplo == CblasUpper) {
457        for( j=0, k=0; j<*n; j++ )
458           for( i=0; i<j+1; i++, k++ )
459              A[ LDA*i+j ]=ap[ k ];
460        for( i=0, k=0; i<*n; i++ )
461           for( j=i; j<*n; j++, k++ )
462              AP[ k ]=A[ LDA*i+j ];
463
464     }
465     else {
466        for( j=0, k=0; j<*n; j++ )
467           for( i=j; i<*n; i++, k++ )
468              A[ LDA*i+j ]=ap[ k ];
469        for( i=0, k=0; i<*n; i++ )
470           for( j=0; j<i+1; j++, k++ )
471              AP[ k ]=A[ LDA*i+j ];
472     }
473     cblas_dtpsv( CblasRowMajor, uplo, trans, diag, *n, AP, x, *incx );
474     free(A);
475     free(AP);
476  }
477  else
478     cblas_dtpsv( CblasColMajor, uplo, trans, diag, *n, ap, x, *incx );
479}
480
481void F77_dspr(int *order, char *uplow, int *n, double *alpha, double *x,
482	     int *incx, double *ap ){
483  double *A, *AP;
484  int i,j,k,LDA;
485  enum CBLAS_UPLO uplo;
486
487  get_uplo_type(uplow,&uplo);
488
489  if (*order == TEST_ROW_MJR) {
490     LDA = *n;
491     A   = ( double* )malloc( LDA*LDA*sizeof( double ) );
492     AP  = ( double* )malloc( (((LDA+1)*LDA)/2)*sizeof( double ) );
493     if (uplo == CblasUpper) {
494        for( j=0, k=0; j<*n; j++ )
495           for( i=0; i<j+1; i++, k++ )
496              A[ LDA*i+j ]=ap[ k ];
497        for( i=0, k=0; i<*n; i++ )
498           for( j=i; j<*n; j++, k++ )
499              AP[ k ]=A[ LDA*i+j ];
500     }
501     else {
502        for( j=0, k=0; j<*n; j++ )
503           for( i=j; i<*n; i++, k++ )
504              A[ LDA*i+j ]=ap[ k ];
505        for( i=0, k=0; i<*n; i++ )
506           for( j=0; j<i+1; j++, k++ )
507              AP[ k ]=A[ LDA*i+j ];
508     }
509     cblas_dspr( CblasRowMajor, uplo, *n, *alpha, x, *incx, AP );
510     if (uplo == CblasUpper) {
511        for( i=0, k=0; i<*n; i++ )
512           for( j=i; j<*n; j++, k++ )
513              A[ LDA*i+j ]=AP[ k ];
514        for( j=0, k=0; j<*n; j++ )
515           for( i=0; i<j+1; i++, k++ )
516              ap[ k ]=A[ LDA*i+j ];
517     }
518     else {
519        for( i=0, k=0; i<*n; i++ )
520           for( j=0; j<i+1; j++, k++ )
521              A[ LDA*i+j ]=AP[ k ];
522        for( j=0, k=0; j<*n; j++ )
523           for( i=j; i<*n; i++, k++ )
524              ap[ k ]=A[ LDA*i+j ];
525     }
526     free(A);
527     free(AP);
528  }
529  else
530     cblas_dspr( CblasColMajor, uplo, *n, *alpha, x, *incx, ap );
531}
532
533void F77_dspr2(int *order, char *uplow, int *n, double *alpha, double *x,
534	     int *incx, double *y, int *incy, double *ap ){
535  double *A, *AP;
536  int i,j,k,LDA;
537  enum CBLAS_UPLO uplo;
538
539  get_uplo_type(uplow,&uplo);
540
541  if (*order == TEST_ROW_MJR) {
542     LDA = *n;
543     A   = ( double* )malloc( LDA*LDA*sizeof( double ) );
544     AP  = ( double* )malloc( (((LDA+1)*LDA)/2)*sizeof( double ) );
545     if (uplo == CblasUpper) {
546        for( j=0, k=0; j<*n; j++ )
547           for( i=0; i<j+1; i++, k++ )
548              A[ LDA*i+j ]=ap[ k ];
549        for( i=0, k=0; i<*n; i++ )
550           for( j=i; j<*n; j++, k++ )
551              AP[ k ]=A[ LDA*i+j ];
552     }
553     else {
554        for( j=0, k=0; j<*n; j++ )
555           for( i=j; i<*n; i++, k++ )
556              A[ LDA*i+j ]=ap[ k ];
557        for( i=0, k=0; i<*n; i++ )
558           for( j=0; j<i+1; j++, k++ )
559              AP[ k ]=A[ LDA*i+j ];
560     }
561     cblas_dspr2( CblasRowMajor, uplo, *n, *alpha, x, *incx, y, *incy, AP );
562     if (uplo == CblasUpper) {
563        for( i=0, k=0; i<*n; i++ )
564           for( j=i; j<*n; j++, k++ )
565              A[ LDA*i+j ]=AP[ k ];
566        for( j=0, k=0; j<*n; j++ )
567           for( i=0; i<j+1; i++, k++ )
568              ap[ k ]=A[ LDA*i+j ];
569     }
570     else {
571        for( i=0, k=0; i<*n; i++ )
572           for( j=0; j<i+1; j++, k++ )
573              A[ LDA*i+j ]=AP[ k ];
574        for( j=0, k=0; j<*n; j++ )
575           for( i=j; i<*n; i++, k++ )
576              ap[ k ]=A[ LDA*i+j ];
577     }
578     free(A);
579     free(AP);
580  }
581  else
582     cblas_dspr2( CblasColMajor, uplo, *n, *alpha, x, *incx, y, *incy, ap );
583}
584