Home | History | Annotate | Download | only in testing
      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 
     11 void 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 
     38 void 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 
     63 void 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 
     91 void 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 }
    115 void 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 
    139 void 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 
    163 void 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 
    187 void 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 
    225 void 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 
    271 void 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 
    317 void 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 
    362 void 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 
    400 void 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 
    440 void 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 
    481 void 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 
    533 void 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