ctest/c_dblas2.c

   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 "common.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 }