LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
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 "cblas.h"
9 #include "cblas_test.h"
10 
11 void F77_dgemv(CBLAS_INT *layout, char *transp, CBLAS_INT *m, CBLAS_INT *n, double *alpha,
12  double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx, double *beta,
13  double *y, CBLAS_INT *incy ) {
14 
15  double *A;
16  CBLAS_INT i,j,LDA;
17  CBLAS_TRANSPOSE trans;
18 
19  get_transpose_type(transp, &trans);
20  if (*layout == 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 (*layout == 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(CBLAS_INT *layout, CBLAS_INT *m, CBLAS_INT *n, double *alpha, double *x, CBLAS_INT *incx,
39  double *y, CBLAS_INT *incy, double *a, CBLAS_INT *lda ) {
40 
41  double *A;
42  CBLAS_INT i,j,LDA;
43 
44  if (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
64  CBLAS_INT *n, double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx) {
65  double *A;
66  CBLAS_INT i,j,LDA;
67  CBLAS_TRANSPOSE trans;
68  CBLAS_UPLO uplo;
69  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 (*layout == 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 (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
92  CBLAS_INT *n, double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx ) {
93  double *A;
94  CBLAS_INT i,j,LDA;
95  CBLAS_TRANSPOSE trans;
96  CBLAS_UPLO uplo;
97  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 (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *a,
116  CBLAS_INT *lda, double *x, CBLAS_INT *incx, double *beta, double *y,
117  CBLAS_INT *incy) {
118  double *A;
119  CBLAS_INT i,j,LDA;
120  CBLAS_UPLO uplo;
121 
122  get_uplo_type(uplow,&uplo);
123 
124  if (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *x,
140  CBLAS_INT *incx, double *a, CBLAS_INT *lda) {
141  double *A;
142  CBLAS_INT i,j,LDA;
143  CBLAS_UPLO uplo;
144 
145  get_uplo_type(uplow,&uplo);
146 
147  if (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *x,
164  CBLAS_INT *incx, double *y, CBLAS_INT *incy, double *a, CBLAS_INT *lda) {
165  double *A;
166  CBLAS_INT i,j,LDA;
167  CBLAS_UPLO uplo;
168 
169  get_uplo_type(uplow,&uplo);
170 
171  if (*layout == 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(CBLAS_INT *layout, char *transp, CBLAS_INT *m, CBLAS_INT *n, CBLAS_INT *kl, CBLAS_INT *ku,
188  double *alpha, double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx,
189  double *beta, double *y, CBLAS_INT *incy ) {
190 
191  double *A;
192  CBLAS_INT i,irow,j,jcol,LDA;
193  CBLAS_TRANSPOSE trans;
194 
195  get_transpose_type(transp, &trans);
196 
197  if (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
226  CBLAS_INT *n, CBLAS_INT *k, double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx) {
227  double *A;
228  CBLAS_INT irow, jcol, i, j, LDA;
229  CBLAS_TRANSPOSE trans;
230  CBLAS_UPLO uplo;
231  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 (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
272  CBLAS_INT *n, CBLAS_INT *k, double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx) {
273  double *A;
274  CBLAS_INT irow, jcol, i, j, LDA;
275  CBLAS_TRANSPOSE trans;
276  CBLAS_UPLO uplo;
277  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 (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, CBLAS_INT *k, double *alpha,
318  double *a, CBLAS_INT *lda, double *x, CBLAS_INT *incx, double *beta,
319  double *y, CBLAS_INT *incy) {
320  double *A;
321  CBLAS_INT i,j,irow,jcol,LDA;
322  CBLAS_UPLO uplo;
323 
324  get_uplo_type(uplow,&uplo);
325 
326  if (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *ap,
363  double *x, CBLAS_INT *incx, double *beta, double *y, CBLAS_INT *incy) {
364  double *A,*AP;
365  CBLAS_INT i,j,k,LDA;
366  CBLAS_UPLO uplo;
367 
368  get_uplo_type(uplow,&uplo);
369 
370  if (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
401  CBLAS_INT *n, double *ap, double *x, CBLAS_INT *incx) {
402  double *A, *AP;
403  CBLAS_INT i, j, k, LDA;
404  CBLAS_TRANSPOSE trans;
405  CBLAS_UPLO uplo;
406  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 (*layout == 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(CBLAS_INT *layout, char *uplow, char *transp, char *diagn,
441  CBLAS_INT *n, double *ap, double *x, CBLAS_INT *incx) {
442  double *A, *AP;
443  CBLAS_INT i, j, k, LDA;
444  CBLAS_TRANSPOSE trans;
445  CBLAS_UPLO uplo;
446  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 (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *x,
482  CBLAS_INT *incx, double *ap ){
483  double *A, *AP;
484  CBLAS_INT i,j,k,LDA;
485  CBLAS_UPLO uplo;
486 
487  get_uplo_type(uplow,&uplo);
488 
489  if (*layout == 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(CBLAS_INT *layout, char *uplow, CBLAS_INT *n, double *alpha, double *x,
534  CBLAS_INT *incx, double *y, CBLAS_INT *incy, double *ap ){
535  double *A, *AP;
536  CBLAS_INT i,j,k,LDA;
537  CBLAS_UPLO uplo;
538 
539  get_uplo_type(uplow,&uplo);
540 
541  if (*layout == 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 }