diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..1f40305 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.vscode/* +build/* diff --git a/blas.h b/blas.h new file mode 100644 index 0000000..f7ca91c --- /dev/null +++ b/blas.h @@ -0,0 +1,33 @@ +#ifndef _BLAS_H +#define _BLAS_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include + +void caxpy ( int n, float complex ca, float complex cx[], int incx, float complex cy[], int incy ); +void ccopy ( int n, float complex cx[], int incx, float complex cy[], int incy ); +float complex cdotc ( int n, float complex cx[], int incx, float complex cy[], int incy ); +float complex cdotu ( int n, float complex cx[], int incx, float complex cy[], int incy ); +void cgemm(char transa, char transb, int m, int n, int k, complex float alpha, + complex float a[], int lda, complex float b[], int ldb, complex float beta, + complex float c[], int ldc); +void crotg ( float complex *ca, float complex cb, float *c, float complex *s ); +void cscal ( int n, float complex ca, float complex cx[], int incx ); +void csrot ( int n, float complex cx[], int incx, float complex cy[], int incy, float c, float s ); +void csscal ( int n, float sa, float complex cx[], int incx ); +void cswap ( int n, float complex cx[], int incx, float complex cy[], int incy ); +int icamax ( int n, float complex x[], int incx ); +float scasum ( int n, float complex x[], int incx ); +float scnrm2 ( int n, float complex x[], int incx ); + + + +#ifdef __cplusplus +} +#endif + + +#endif \ No newline at end of file diff --git a/src/caxpy.c b/src/caxpy.c index 8dfdba9..b2b2397 100644 --- a/src/caxpy.c +++ b/src/caxpy.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CAXPY * * =========== DOCUMENTATION =========== @@ -84,56 +89,51 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX CA - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - COMPLEX CX(*),CY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I,IX,IY -* .. -* .. External Functions .. - REAL SCABS1 - EXTERNAL SCABS1 -* .. - IF (N.LE.0) RETURN - IF (SCABS1(CA).EQ.0.0E+0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* - DO I = 1,N - CY(I) = CY(I) + CA*CX(I) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - CY(IY) = CY(IY) + CA*CX(IX) - IX = IX + INCX - IY = IY + INCY - END DO - END IF -* - RETURN -* -* End of CAXPY -* - END +*/ +void caxpy ( int n, float complex ca, float complex cx[], int incx, + float complex cy[], int incy ) { +// +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 +// + int i, ix, iy; + + if ( n <= 0 ) return; + + + if ( scabs1 ( ca ) == 0.0 ) return; + + if ( incx != 1 || incy != 1 ) { + if ( 0 <= incx ) { + ix = 0; + } + else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } + else { + iy = ( -n + 1 ) * incy; + } + + for (i = 0; i < n; i++) { + cy[iy] = cy[iy] + ca * cx[ix]; + ix = ix + incx; + iy = iy + incy; + } + } + else { + for (i = 0; i < n; i++) { + cy[i] = cy[i] + ca * cx[i]; + } + } + + return; + + // End of CAXPY + +} diff --git a/src/ccopy.c b/src/ccopy.c index 1b4999b..369f366 100644 --- a/src/ccopy.c +++ b/src/ccopy.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CCOPY * * =========== DOCUMENTATION =========== @@ -77,49 +82,46 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CCOPY(N,CX,INCX,CY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - COMPLEX CX(*),CY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I,IX,IY -* .. - IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* - DO I = 1,N - CY(I) = CX(I) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - CY(IY) = CX(IX) - IX = IX + INCX - IY = IY + INCY - END DO - END IF - RETURN -* -* End of CCOPY -* - END +*/ +void ccopy ( int n, float complex cx[], int incx, float complex cy[], + int incy ) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + + int i, ix, iy; + + if ( n <= 0 ) { + return; + } + + if ( incx != 1 || incy != 1 ) { + if ( 0 <= incx ) { + ix = 0; + } else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } else { + iy = ( -n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + cy[iy] = cx[ix]; + ix = ix + incx; + iy = iy + incy; + } + } else { + for ( i = 0; i < n; i++ ) { + cy[i] = cx[i]; + } + } + return; + + // End of CCOPY + +} diff --git a/src/cdotc.c b/src/cdotc.c index da29fdf..afc2274 100644 --- a/src/cdotc.c +++ b/src/cdotc.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CDOTC * * =========== DOCUMENTATION =========== @@ -79,56 +84,51 @@ *> \endverbatim *> * ===================================================================== - COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - COMPLEX CX(*),CY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - COMPLEX CTEMP - INTEGER I,IX,IY -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG -* .. - CTEMP = (0.0,0.0) - CDOTC = (0.0,0.0) - IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* - DO I = 1,N - CTEMP = CTEMP + CONJG(CX(I))*CY(I) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - CTEMP = CTEMP + CONJG(CX(IX))*CY(IY) - IX = IX + INCX - IY = IY + INCY - END DO - END IF - CDOTC = CTEMP - RETURN -* -* End of CDOTC -* - END +*/ +float complex cdotc ( int n, float complex cx[], int incx, + float complex cy[], int incy ) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + + + int i, ix, iy; + float complex value; + + value = 0.0; + + if ( n <= 0 ) { + return value; + } + + if ( incx == 1 && incy == 1 ) { + for ( i = 0; i < n; i++ ) { + value = value + ( ~cx[i] ) * cy[i]; + } + } else { + if ( 0 <= incx ) { + ix = 0; + } else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } else { + iy = ( -n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + value = value + ( ~cx[ix] ) * cy[iy]; + ix = ix + incx; + iy = iy + incy; + } + } + + return value; + + // End of CDOTC + +} diff --git a/src/cdotu.c b/src/cdotu.c index d8c21d1..35ff563 100644 --- a/src/cdotu.c +++ b/src/cdotu.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CDOTU * * =========== DOCUMENTATION =========== @@ -79,53 +84,51 @@ *> \endverbatim *> * ===================================================================== - COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - COMPLEX CX(*),CY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - COMPLEX CTEMP - INTEGER I,IX,IY -* .. - CTEMP = (0.0,0.0) - CDOTU = (0.0,0.0) - IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* - DO I = 1,N - CTEMP = CTEMP + CX(I)*CY(I) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - CTEMP = CTEMP + CX(IX)*CY(IY) - IX = IX + INCX - IY = IY + INCY - END DO - END IF - CDOTU = CTEMP - RETURN -* -* End of CDOTU -* - END +*/ +float complex cdotu ( int n, float complex cx[], int incx, + float complex cy[], int incy ) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + + + int i, ix, iy; + float complex value; + + value = 0.0; + + if ( n <= 0 ) { + return value; + } + + if ( incx == 1 && incy == 1 ) { + for ( i = 0; i < n; i++ ) { + value = value + cx[i] * cy[i]; + } + } else { + if ( 0 <= incx ) { + ix = 0; + } else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } else { + iy = ( -n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + value = value + cx[ix] * cy[iy]; + ix = ix + incx; + iy = iy + incy; + } + } + + return value; + + // End of CDOTU + +} diff --git a/src/cgemm.c b/src/cgemm.c index baefe21..42e31c4 100644 --- a/src/cgemm.c +++ b/src/cgemm.c @@ -1,3 +1,11 @@ +#include +#include +#include +#include + +#include "blas.h" + +/* *> \brief \b CGEMM * * =========== DOCUMENTATION =========== @@ -183,295 +191,319 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) -* -* -- Reference BLAS level3 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX ALPHA,BETA - INTEGER K,LDA,LDB,LDC,M,N - CHARACTER TRANSA,TRANSB -* .. -* .. Array Arguments .. - COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) -* .. -* -* ===================================================================== -* -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,MAX -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,J,L,NROWA,NROWB - LOGICAL CONJA,CONJB,NOTA,NOTB -* .. -* .. Parameters .. - COMPLEX ONE - PARAMETER (ONE= (1.0E+0,0.0E+0)) - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* -* Set NOTA and NOTB as true if A and B respectively are not -* conjugated or transposed, set CONJA and CONJB as true if A and -* B respectively are to be transposed but not conjugated and set -* NROWA and NROWB as the number of rows of A and B respectively. -* - NOTA = LSAME(TRANSA,'N') - NOTB = LSAME(TRANSB,'N') - CONJA = LSAME(TRANSA,'C') - CONJB = LSAME(TRANSB,'C') - IF (NOTA) THEN - NROWA = M - ELSE - NROWA = K - END IF - IF (NOTB) THEN - NROWB = K - ELSE - NROWB = N - END IF -* -* Test the input parameters. -* - INFO = 0 - IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND. - + (.NOT.LSAME(TRANSA,'T'))) THEN - INFO = 1 - ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND. - + (.NOT.LSAME(TRANSB,'T'))) THEN - INFO = 2 - ELSE IF (M.LT.0) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (K.LT.0) THEN - INFO = 5 - ELSE IF (LDA.LT.MAX(1,NROWA)) THEN - INFO = 8 - ELSE IF (LDB.LT.MAX(1,NROWB)) THEN - INFO = 10 - ELSE IF (LDC.LT.MAX(1,M)) THEN - INFO = 13 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CGEMM ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. - + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN -* -* And when alpha.eq.zero. -* - IF (ALPHA.EQ.ZERO) THEN - IF (BETA.EQ.ZERO) THEN - DO 20 J = 1,N - DO 10 I = 1,M - C(I,J) = ZERO - 10 CONTINUE - 20 CONTINUE - ELSE - DO 40 J = 1,N - DO 30 I = 1,M - C(I,J) = BETA*C(I,J) - 30 CONTINUE - 40 CONTINUE - END IF - RETURN - END IF -* -* Start the operations. -* - IF (NOTB) THEN - IF (NOTA) THEN -* -* Form C := alpha*A*B + beta*C. -* - DO 90 J = 1,N - IF (BETA.EQ.ZERO) THEN - DO 50 I = 1,M - C(I,J) = ZERO - 50 CONTINUE - ELSE IF (BETA.NE.ONE) THEN - DO 60 I = 1,M - C(I,J) = BETA*C(I,J) - 60 CONTINUE - END IF - DO 80 L = 1,K - TEMP = ALPHA*B(L,J) - DO 70 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 70 CONTINUE - 80 CONTINUE - 90 CONTINUE - ELSE IF (CONJA) THEN -* -* Form C := alpha*A**H*B + beta*C. -* - DO 120 J = 1,N - DO 110 I = 1,M - TEMP = ZERO - DO 100 L = 1,K - TEMP = TEMP + CONJG(A(L,I))*B(L,J) - 100 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 110 CONTINUE - 120 CONTINUE - ELSE -* -* Form C := alpha*A**T*B + beta*C -* - DO 150 J = 1,N - DO 140 I = 1,M - TEMP = ZERO - DO 130 L = 1,K - TEMP = TEMP + A(L,I)*B(L,J) - 130 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 140 CONTINUE - 150 CONTINUE - END IF - ELSE IF (NOTA) THEN - IF (CONJB) THEN -* -* Form C := alpha*A*B**H + beta*C. -* - DO 200 J = 1,N - IF (BETA.EQ.ZERO) THEN - DO 160 I = 1,M - C(I,J) = ZERO - 160 CONTINUE - ELSE IF (BETA.NE.ONE) THEN - DO 170 I = 1,M - C(I,J) = BETA*C(I,J) - 170 CONTINUE - END IF - DO 190 L = 1,K - TEMP = ALPHA*CONJG(B(J,L)) - DO 180 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 180 CONTINUE - 190 CONTINUE - 200 CONTINUE - ELSE -* -* Form C := alpha*A*B**T + beta*C -* - DO 250 J = 1,N - IF (BETA.EQ.ZERO) THEN - DO 210 I = 1,M - C(I,J) = ZERO - 210 CONTINUE - ELSE IF (BETA.NE.ONE) THEN - DO 220 I = 1,M - C(I,J) = BETA*C(I,J) - 220 CONTINUE - END IF - DO 240 L = 1,K - TEMP = ALPHA*B(J,L) - DO 230 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 230 CONTINUE - 240 CONTINUE - 250 CONTINUE - END IF - ELSE IF (CONJA) THEN - IF (CONJB) THEN -* -* Form C := alpha*A**H*B**H + beta*C. -* - DO 280 J = 1,N - DO 270 I = 1,M - TEMP = ZERO - DO 260 L = 1,K - TEMP = TEMP + CONJG(A(L,I))*CONJG(B(J,L)) - 260 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 270 CONTINUE - 280 CONTINUE - ELSE -* -* Form C := alpha*A**H*B**T + beta*C -* - DO 310 J = 1,N - DO 300 I = 1,M - TEMP = ZERO - DO 290 L = 1,K - TEMP = TEMP + CONJG(A(L,I))*B(J,L) - 290 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 300 CONTINUE - 310 CONTINUE - END IF - ELSE - IF (CONJB) THEN -* -* Form C := alpha*A**T*B**H + beta*C -* - DO 340 J = 1,N - DO 330 I = 1,M - TEMP = ZERO - DO 320 L = 1,K - TEMP = TEMP + A(L,I)*CONJG(B(J,L)) - 320 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 330 CONTINUE - 340 CONTINUE - ELSE -* -* Form C := alpha*A**T*B**T + beta*C -* - DO 370 J = 1,N - DO 360 I = 1,M - TEMP = ZERO - DO 350 L = 1,K - TEMP = TEMP + A(L,I)*B(J,L) - 350 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 360 CONTINUE - 370 CONTINUE - END IF - END IF -* - RETURN -* -* End of CGEMM -* - END +*/ +void cgemm(char transa, char transb, int m, int n, int k, complex float alpha, + complex float a[], int lda, complex float b[], int ldb, complex float beta, + complex float c[], int ldc) { + +// -- Reference BLAS level3 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + +// Local Scalars .. + int i, j, l, nrowa, nrowb, info, ncola; + complex float temp; + bool conja, conjb, nota, notb; + + complex float one = (1.0e0, 0.0e0); + complex float zero = (0.0e0, 0.0e0); + +// Set NOTA and NOTB as true if A and B respectively are not +// conjugated or transposed, set CONJA and CONJB as true if A and +// B respectively are to be transposed but not conjugated and set +// NROWA, NCOLA and NROWB as the number of rows and columns of A +// and the number of rows of B respectively. + + nota = ( ( transa == 'N' ) || ( transa == 'n' ) ); + notb = ( ( transb == 'N' ) || ( transb == 'n' ) ); + conja = ( ( transa == 'C' ) || ( transa == 'c' ) ); + conjb = ( ( transb == 'C' ) || ( transb == 'c' ) ); + + if ( nota ) { + nrowa = m; + ncola = k; + } + else { + nrowa = k; + ncola = m; + } + + if (notb) { + nrowb = k; + } + else { + nrowb = n; + } + +// Test the input parameters. + + info = 0; + if ( ! ( transa == 'N' || transa == 'n' || + transa == 'C' || transa == 'c' || + transa == 'T' || transa == 't' ) ) { + info = 1; + } + + if ( ! ( transb == 'N' || transb == 'n' || + transb == 'C' || transb == 'c' || + transb == 'T' || transb == 't' ) ) { + info = 2; + } + + if ( m < 0 ) { + info = 3; + } + else if ( n < 0 ) { + info = 4; + } + else if ( k < 0 ) { + info = 5; + } + else if ( lda < i4_max( 1, nrowa ) ) { + info = 8; + } + else if ( ldb < i4_max( 1, nrowb ) ) { + info = 10; + } + else if ( ldc < i4_max( 1, m ) ) { + info = 13; + } + + if (info != 0) { + xerbla('cgemm ', info); + return; + } + +// Quick return if possible. + + if (( m == 0 ) || (n == 0) || (((alpha == zero) || ( k == 0 )) && ( beta == one ) ) ) { + return; + } + +// And when alpha.eq.zero. + + if ( alpha == zero ) { + if ( beta == zero ) { + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = zero; + } + } + } + else { + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = beta * c[i+j*ldc]; + } + } + } + return; + } + +// Start the operations. + + if ( notb ) { + if ( nota ) { + +// Form C := alpha*A*B + beta*C. + + for ( j = 0; j < n; j++ ) { + if ( beta == zero ) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = zero; + } + } + else if ( beta != one ) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = beta*c[i+j*ldc]; + } + } + + for ( l = 0; l < k; l++ ) { + if ( b[l+j*ldb] != zero ) { + temp = alpha * b[l+j*ldb]; + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = c[i+j*ldc] + temp * a[i+l*lda]; + } + } + } + } + } + else if (conja) { + +// Form C := alpha*A**H*B + beta*C. + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + (~a[l+i*lda])* b[l+j*ldb]; + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + else { + +// Form C := alpha*A**T*B + beta*C + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + a[l+i*lda] * b[l+j*ldb]; + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + } + else if (nota) { + if (conjb) { + +// Form C := alpha*A*B**H + beta*C. + + for ( j = 0; j < n; j++ ) { + if (beta == zero) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = zero; + } + } + else if (beta != one) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = beta*c[i+j*ldc]; + } + } + for ( l = 0; l < k; l++ ) { + if (b[j+l*ldb] != zero) { + temp = alpha*(~b[j+l*ldb]); + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = c[i+j*ldc] + temp*a[i+l*lda]; + } + } + } + } + } + else { + +// Form C := alpha*A*B**T + beta*C + + for ( j = 0; j < n; j++ ) { + if (beta == zero) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = zero; + } + } + else if (beta != one) { + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = beta*c[i+j*ldc]; + } + } + for ( l = 0; l < k; l++ ) { + if (b[j+l*ldb] != zero) { + temp = alpha*b[j+l*ldb]; + for ( i = 0; i < m; i++ ) { + c[i+j*ldc] = c[i+j*ldc] + temp*a[i+l*lda]; + } + } + } + } + } + } + else if (conja) { + if (conjb) { + +// Form C := alpha*A**H*B**H + beta*C. + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + (~a[l+i*lda])*(~b[j+l*ldb]); + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + else { + +// Form C := alpha*A**H*B**T + beta*C + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + (~a[l+i*lda])*b[j+l*ldb]; + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + } + else { + if (conjb) { + +// Form C := alpha*A**T*B**H + beta*C + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + a[l+i*lda]*(~b[j+l*ldb]); + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + else { + +// Form C := alpha*A**T*B**T + beta*C + + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = zero; + for ( l = 0; l < k; l++ ) { + temp = temp + a[l+i*lda]*b[j+l*ldb]; + } + if (beta == zero) { + c[i+j*ldc] = alpha*temp; + } + else { + c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc]; + } + } + } + } + } + + return; + + // End of CGEMM . + +} \ No newline at end of file diff --git a/src/cgemv.c b/src/cgemv.c index 574be07..5dc0370 100644 --- a/src/cgemv.c +++ b/src/cgemv.c @@ -1,3 +1,9 @@ +#include +#include + +#include "blas.h" + +/* *> \brief \b CGEMV * * =========== DOCUMENTATION =========== @@ -154,194 +160,185 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX ALPHA,BETA - INTEGER INCX,INCY,LDA,M,N - CHARACTER TRANS -* .. -* .. Array Arguments .. - COMPLEX A(LDA,*),X(*),Y(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ONE - PARAMETER (ONE= (1.0E+0,0.0E+0)) - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY - LOGICAL NOCONJ -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 1 - ELSE IF (M.LT.0) THEN - INFO = 2 - ELSE IF (N.LT.0) THEN - INFO = 3 - ELSE IF (LDA.LT.MAX(1,M)) THEN - INFO = 6 - ELSE IF (INCX.EQ.0) THEN - INFO = 8 - ELSE IF (INCY.EQ.0) THEN - INFO = 11 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CGEMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. - + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN -* - NOCONJ = LSAME(TRANS,'T') -* -* Set LENX and LENY, the lengths of the vectors x and y, and set -* up the start points in X and Y. -* - IF (LSAME(TRANS,'N')) THEN - LENX = N - LENY = M - ELSE - LENX = M - LENY = N - END IF - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (LENX-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (LENY-1)*INCY - END IF -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* -* First form y := beta*y. -* - IF (BETA.NE.ONE) THEN - IF (INCY.EQ.1) THEN - IF (BETA.EQ.ZERO) THEN - DO 10 I = 1,LENY - Y(I) = ZERO - 10 CONTINUE - ELSE - DO 20 I = 1,LENY - Y(I) = BETA*Y(I) - 20 CONTINUE - END IF - ELSE - IY = KY - IF (BETA.EQ.ZERO) THEN - DO 30 I = 1,LENY - Y(IY) = ZERO - IY = IY + INCY - 30 CONTINUE - ELSE - DO 40 I = 1,LENY - Y(IY) = BETA*Y(IY) - IY = IY + INCY - 40 CONTINUE - END IF - END IF - END IF - IF (ALPHA.EQ.ZERO) RETURN - IF (LSAME(TRANS,'N')) THEN -* -* Form y := alpha*A*x + y. -* - JX = KX - IF (INCY.EQ.1) THEN - DO 60 J = 1,N - TEMP = ALPHA*X(JX) - DO 50 I = 1,M - Y(I) = Y(I) + TEMP*A(I,J) - 50 CONTINUE - JX = JX + INCX - 60 CONTINUE - ELSE - DO 80 J = 1,N - TEMP = ALPHA*X(JX) - IY = KY - DO 70 I = 1,M - Y(IY) = Y(IY) + TEMP*A(I,J) - IY = IY + INCY - 70 CONTINUE - JX = JX + INCX - 80 CONTINUE - END IF - ELSE -* -* Form y := alpha*A**T*x + y or y := alpha*A**H*x + y. -* - JY = KY - IF (INCX.EQ.1) THEN - DO 110 J = 1,N - TEMP = ZERO - IF (NOCONJ) THEN - DO 90 I = 1,M - TEMP = TEMP + A(I,J)*X(I) - 90 CONTINUE - ELSE - DO 100 I = 1,M - TEMP = TEMP + CONJG(A(I,J))*X(I) - 100 CONTINUE - END IF - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 110 CONTINUE - ELSE - DO 140 J = 1,N - TEMP = ZERO - IX = KX - IF (NOCONJ) THEN - DO 120 I = 1,M - TEMP = TEMP + A(I,J)*X(IX) - IX = IX + INCX - 120 CONTINUE - ELSE - DO 130 I = 1,M - TEMP = TEMP + CONJG(A(I,J))*X(IX) - IX = IX + INCX - 130 CONTINUE - END IF - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 140 CONTINUE - END IF - END IF -* - RETURN -* -* End of CGEMV -* - END +*/ +void cgemv(char trans, int m, int n, complex float alpha, complex float a[], int lda, complex float x[], int incx, complex float beta, complex float y[], int incy) { + +// -- Reference BLAS level2 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + + +// .. + +// ===================================================================== + +// .. Parameters .. + complex float one = (1.0E+0,0.0E+0 * I); + complex float zero = (0.0E+0,0.0E+0 * I); +// .. +// .. Local Scalars .. + complex float temp; + int i,info,ix,iy,j,jx,jy,kx,ky,lenx,leny; + bool noconj; + +// .. + +// Test the input parameters. + + info = 0; + if (!lsame(trans,'N') && !lsame(trans,'T') && + + !lsame(trans,'C')) { + info = 1; + } else if (m < 0) { + info = 2; + } else if (n < 0) { + info = 3; + } else if (lda < MAX(1,m)) { + info = 6; + } else if (incx == 0) { + info = 8; + } else if (incy == 0) { + info = 11; + } + if (info != 0) { + xerbla('CGEMV ',info); + return; + } + +// Quick return if possible. + + if ((m == 0) || (n == 0) || + ((alpha == zero) && (beta == one))) { + return; + } + + noconj = lsame(trans,'T'); + +// Set LENX and LENY, the lengths of the vectors x and y, and set +// up the start points in X and Y. + + if (lsame(trans,'N')) { + lenx = n; + leny = m; + } else { + lenx = m; + leny = n; + } + if (incx > 0) { + kx = 0; + } else { + kx = 0 - (lenx-1)*incx; + } + if (incy > 0) { + ky = 0; + } else { + ky = 0 - (leny-1)*incy; + } + +// Start the operations. In this version the elements of A are +// accessed sequentially with one pass through A. + +// First form y := beta*y. + + if (beta != one) { + if (incy == 1) { + if (beta == zero) { + for ( i = 0; i < leny; i++) { + y[i] = zero; + } + } else { + for ( i = 0; i < leny; i++) { + y[i] = beta*y[i]; + } + } + } else { + iy = ky; + if (beta == zero) { + for ( i = 0; i < leny; i++) { + y[iy] = zero; + iy = iy + incy; + } + } else { + for ( i = 0; i < leny; i++) { + y[iy] = beta*y[iy]; + iy = iy + incy; + } + } + } + } + if (alpha == zero) {return;} + if (lsame(trans,'N')) { + +// Form y := alpha*A*x + y. + + jx = kx; + if (incy == 1) { + for ( j = 0; j < n; j++) { + if (x[jx] != zero) { + temp = alpha*x[jx]; + for ( i = 0; i < m; i++) { + y[i] = y[i] + temp*a[i+j*lda]; + } + } + jx = jx + incx; + } + } else { + for ( j = 0; j < n; j++) { + if (x[jx] != zero) { + temp = alpha*x[jx]; + iy = ky; + for ( i = 0; i < m; i++) { + y[iy] = y[iy] + temp*a[i+j*lda]; + iy = iy + incy; + } + } + jx = jx + incx; + } + } + } else { + +// Form y := alpha*A**T*x + y or y := alpha*A**H*x + y. + + jy = ky; + if (incx == 1) { + for ( j = 0; j < n; j++) { + temp = zero; + if (noconj) { + for ( i = 0; i < m; i++) { + temp = temp + a[i+j*lda]*x[i]; + } + } else { + for ( i = 0; i < m; i++) { + temp = temp + (~a[i+j*lda])*x[i]; + } + } + y[jy] = y[jy] + alpha*temp; + jy = jy + incy; + } + } else { + for ( j = 0; j < n; j++) { + temp = zero; + ix = kx; + if (noconj) { + for ( i = 0; i < m; i++) { + temp = temp + a[i+j*lda]*x[ix]; + ix = ix + incx; + } + } else { + for ( i = 0; i < m; i++) { + temp = temp + (~a[i+j*lda])*x[ix]; + ix = ix + incx; + } + } + y[jy] = y[jy] + alpha*temp; + jy = jy + incy; + } + } + } + + return; + + // End of CGEMV . + +} \ No newline at end of file diff --git a/src/cgerc.c b/src/cgerc.c index 716628d..e17a133 100644 --- a/src/cgerc.c +++ b/src/cgerc.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CGERC * * =========== DOCUMENTATION =========== @@ -126,99 +131,89 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CGERC(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX ALPHA - INTEGER INCX,INCY,LDA,M,N -* .. -* .. Array Arguments .. - COMPLEX A(LDA,*),X(*),Y(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,J,JY,KX -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (M.LT.0) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - ELSE IF (LDA.LT.MAX(1,M)) THEN - INFO = 9 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CGERC ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* - IF (INCY.GT.0) THEN - JY = 1 - ELSE - JY = 1 - (N-1)*INCY - END IF - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*CONJG(Y(JY)) - DO 10 I = 1,M - A(I,J) = A(I,J) + X(I)*TEMP - 10 CONTINUE - END IF - JY = JY + INCY - 20 CONTINUE - ELSE - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (M-1)*INCX - END IF - DO 40 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*CONJG(Y(JY)) - IX = KX - DO 30 I = 1,M - A(I,J) = A(I,J) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - END IF - JY = JY + INCY - 40 CONTINUE - END IF -* - RETURN -* -* End of CGERC -* - END +*/ +void cgerc(int m, int n, complex float alpha, complex float x[], int incx, + complex float y[], int incy, complex float a[], int lda) { + +// -- Reference BLAS level2 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + +// ===================================================================== + +// .. Parameters .. + complex float zero = (0.0E+0, 0.0E+0 * I); + +// .. Local Scalars .. + complex float temp; + int i,info,ix,j,jy,kx; + +// Test the input parameters. + + info = 0; + if (m < 0) { + info = 1; + } else if (n < 0) { + info = 2; + } else if (incx == 0) { + info = 5; + } else if (incy == 0) { + info = 7; + } else if (lda < max(1,m)) { + info = 9; + } + if (info != 0) { + xerbla('CGERC ',info); + return; + } + +// Quick return if possible. + + if ((m == 0) || (n == 0) || (alpha == zero)) return; + +// Start the operations. In this version the elements of A are +// accessed sequentially with one pass through A. + + if (incy > 0) { + jy = 0; + } + else { + jy = 0 - (n-1)*incy; + } + if (incx == 1) { + for ( j = 0; j < n; j++ ) { + if (y[jy] != zero) { + temp = alpha*(~y[jy]); + for ( i = 0; i < m; i++ ) { + a[i+j*lda] = a[i+j*lda] + x[i]*temp; + } + } + jy = jy + incy; + } + } + else { + if (incx > 0) { + kx = 0; + } + else { + kx = 0 - (m-1)*incx; + } + for ( j = 0; j < n; j++ ) { + if (y[jy] != zero) { + temp = alpha*(~y[jy]); + ix = kx; + for ( i = 0; i < m; i++ ) { + a[i+j*lda] = a[i+j*lda] + x[ix]*temp; + ix = ix + incx; + } + } + jy = jy + incy; + } + } + + return; + + // End of CGERC . + +} \ No newline at end of file diff --git a/src/cscal.c b/src/cscal.c index 4110988..85a4d8b 100644 --- a/src/cscal.c +++ b/src/cscal.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CSCAL * * =========== DOCUMENTATION =========== @@ -74,48 +79,41 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CSCAL(N,CA,CX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX CA - INTEGER INCX,N -* .. -* .. Array Arguments .. - COMPLEX CX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I,NINCX -* .. -* .. Parameters .. - COMPLEX ONE - PARAMETER (ONE= (1.0E+0,0.0E+0)) -* .. - IF (N.LE.0 .OR. INCX.LE.0 .OR. CA.EQ.ONE) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - DO I = 1,N - CX(I) = CA*CX(I) - END DO - ELSE -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO I = 1,NINCX,INCX - CX(I) = CA*CX(I) - END DO - END IF - RETURN -* -* End of CSCAL -* - END +*/ +void cscal(int n, float complex ca, float complex *cx, int incx) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 +// +// ===================================================================== + +// .. Local Scalars .. + int i; + + if ( n <= 0 || incx <= 0 ) { + return; + } + + if (incx == 1) { +// +// code for increment equal to 1 +// + for (i = 0; i < n; i++) { + cx[i] = ca * cx[i]; + } + } else { +// +// code for increment not equal to 1 +// + for (i = 0; i < n; i++) { + cx[i*incx] = ca * cx[i*incx]; + } + } + return; + + // End of CSCAL + +} + diff --git a/src/csrot.c b/src/csrot.c index 2b4b92b..2e354ef 100644 --- a/src/csrot.c +++ b/src/csrot.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CSROT * * =========== DOCUMENTATION =========== @@ -94,60 +99,61 @@ *> \ingroup complex_blas_level1 * * ===================================================================== - SUBROUTINE CSROT( N, CX, INCX, CY, INCY, C, S ) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX, INCY, N - REAL C, S -* .. -* .. Array Arguments .. - COMPLEX CX( * ), CY( * ) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I, IX, IY - COMPLEX CTEMP -* .. -* .. Executable Statements .. -* - IF( N.LE.0 ) - $ RETURN - IF( INCX.EQ.1 .AND. INCY.EQ.1 ) THEN -* -* code for both increments equal to 1 -* - DO I = 1, N - CTEMP = C*CX( I ) + S*CY( I ) - CY( I ) = C*CY( I ) - S*CX( I ) - CX( I ) = CTEMP - END DO - ELSE -* -* code for unequal increments or equal increments not equal -* to 1 -* - IX = 1 - IY = 1 - IF( INCX.LT.0 ) - $ IX = ( -N+1 )*INCX + 1 - IF( INCY.LT.0 ) - $ IY = ( -N+1 )*INCY + 1 - DO I = 1, N - CTEMP = C*CX( IX ) + S*CY( IY ) - CY( IY ) = C*CY( IY ) - S*CX( IX ) - CX( IX ) = CTEMP - IX = IX + INCX - IY = IY + INCY - END DO - END IF - RETURN -* -* End of CSROT -* - END +*/ +void csrot ( int n, float complex cx[], int incx, float complex cy[], + int incy, float c, float s ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + float complex ctemp; + int i, ix, iy; + + if ( n <= 0 ) { + return; + } + + if ( incx == 1 && incy == 1 ) { + + // code for both increments equal to 1 + + for ( i = 0; i < n; i++ ) { + ctemp = c * cx[i] + s * cy[i]; + cy[i] = c * cy[i] - s * cx[i]; + cx[i] = ctemp; + } + } else { + + // code for unequal increments or equal increments not equal + // to 1 + + if ( 0 <= incx ) { + ix = 0; + } else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } else { + iy = ( -n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + ctemp = c * cx[ix] + s * cy[iy]; + cy[iy] = c * cy[iy] - s * cx[ix]; + cx[ix] = ctemp; + ix = ix + incx; + iy = iy + incy; + } + } + + return; + + // End of CSROT + +} diff --git a/src/csscal.c b/src/csscal.c index b9a8ca5..7467b1e 100644 --- a/src/csscal.c +++ b/src/csscal.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CSSCAL * * =========== DOCUMENTATION =========== @@ -74,51 +79,40 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CSSCAL(N,SA,CX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL SA - INTEGER INCX,N -* .. -* .. Array Arguments .. - COMPLEX CX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I,NINCX -* .. -* .. Parameters .. - REAL ONE - PARAMETER (ONE=1.0E+0) -* .. -* .. Intrinsic Functions .. - INTRINSIC AIMAG,CMPLX,REAL -* .. - IF (N.LE.0 .OR. INCX.LE.0 .OR. SA.EQ.ONE) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - DO I = 1,N - CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) - END DO - ELSE -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO I = 1,NINCX,INCX - CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) - END DO - END IF - RETURN -* -* End of CSSCAL -* - END +*/ +void csscal ( int n, float sa, float complex cx[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + int i; + + if ( n <= 0 || incx <= 0 ) { + return; + } + + if ( incx == 1 ) { + + // code for increment equal to 1 + + for ( i = 0; i < n; i++ ) { + cx[i] = sa * cx[i]; + } + } else { + + // code for increment not equal to 1 + + for ( i = 0; i < n; i++ ) { + cx[i*incx] = sa * cx[i*incx]; + } + } + + return; + + // End of CSSCAL + +} diff --git a/src/cswap.c b/src/cswap.c index 310bf18..64d6b9d 100644 --- a/src/cswap.c +++ b/src/cswap.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b CSWAP * * =========== DOCUMENTATION =========== @@ -77,53 +82,56 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CSWAP(N,CX,INCX,CY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - COMPLEX CX(*),CY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - COMPLEX CTEMP - INTEGER I,IX,IY -* .. - IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 - DO I = 1,N - CTEMP = CX(I) - CX(I) = CY(I) - CY(I) = CTEMP - END DO - ELSE -* -* code for unequal increments or equal increments not equal -* to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - CTEMP = CX(IX) - CX(IX) = CY(IY) - CY(IY) = CTEMP - IX = IX + INCX - IY = IY + INCY - END DO - END IF - RETURN -* -* End of CSWAP -* - END +*/ +void cswap ( int n, float complex cx[], int incx, float complex cy[], + int incy ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// ===================================================================== + +// .. Local Scalars .. + float complex ctemp; + int i, ix, iy; + + if ( n <= 0 ) return; + + if ( incx == 1 && incy == 1 ) { + // code for both increments equal to 1 + for ( i = 0; i < n; i++ ) { + ctemp = cx[i]; + cx[i] = cy[i]; + cy[i] = ctemp; + } + } else { + + // code for unequal increments or equal increments not equal + // to 1 + + if ( 0 <= incx ) { + ix = 0; + } else { + ix = ( -n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } else { + iy = ( -n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + ctemp = cx[ix]; + cx[ix] = cy[iy]; + cy[iy] = ctemp; + ix = ix + incx; + iy = iy + incy; + } + } + + return; + + // End of CSWAP + +} diff --git a/src/ctrmv.c b/src/ctrmv.c index 2404b3d..bc67b4a 100644 --- a/src/ctrmv.c +++ b/src/ctrmv.c @@ -1,3 +1,9 @@ +#include +#include + +#include "blas.h" + +/* *> \brief \b CTRMV * * =========== DOCUMENTATION =========== @@ -143,228 +149,231 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE CTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,LDA,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - COMPLEX A(LDA,*),X(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,KX - LOGICAL NOCONJ,NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (LDA.LT.MAX(1,N)) THEN - INFO = 6 - ELSE IF (INCX.EQ.0) THEN - INFO = 8 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CTRMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOCONJ = LSAME(TRANS,'T') - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := A*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*A(I,J) - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*A(J,J) - END IF - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 I = 1,J - 1 - X(IX) = X(IX) + TEMP*A(I,J) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*A(J,J) - END IF - JX = JX + INCX - 40 CONTINUE - END IF - ELSE - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*A(I,J) - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*A(J,J) - END IF - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 I = N,J + 1,-1 - X(IX) = X(IX) + TEMP*A(I,J) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*A(J,J) - END IF - JX = JX - INCX - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A**T*x or x := A**H*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 110 J = N,1,-1 - TEMP = X(J) - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 90 I = J - 1,1,-1 - TEMP = TEMP + A(I,J)*X(I) - 90 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J)) - DO 100 I = J - 1,1,-1 - TEMP = TEMP + CONJG(A(I,J))*X(I) - 100 CONTINUE - END IF - X(J) = TEMP - 110 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 140 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 120 I = J - 1,1,-1 - IX = IX - INCX - TEMP = TEMP + A(I,J)*X(IX) - 120 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J)) - DO 130 I = J - 1,1,-1 - IX = IX - INCX - TEMP = TEMP + CONJG(A(I,J))*X(IX) - 130 CONTINUE - END IF - X(JX) = TEMP - JX = JX - INCX - 140 CONTINUE - END IF - ELSE - IF (INCX.EQ.1) THEN - DO 170 J = 1,N - TEMP = X(J) - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 150 I = J + 1,N - TEMP = TEMP + A(I,J)*X(I) - 150 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J)) - DO 160 I = J + 1,N - TEMP = TEMP + CONJG(A(I,J))*X(I) - 160 CONTINUE - END IF - X(J) = TEMP - 170 CONTINUE - ELSE - JX = KX - DO 200 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 180 I = J + 1,N - IX = IX + INCX - TEMP = TEMP + A(I,J)*X(IX) - 180 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J)) - DO 190 I = J + 1,N - IX = IX + INCX - TEMP = TEMP + CONJG(A(I,J))*X(IX) - 190 CONTINUE - END IF - X(JX) = TEMP - JX = JX + INCX - 200 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of CTRMV -* - END +*/ +void ctrmv(char uplo, char trans, char diag, int n, complex float a[], int lda, complex float x[], int incx) { + +// -- Reference BLAS level2 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Parameters .. + complex float zero = (0.0E+0,0.0E+0*I); +// .. +// .. Local Scalars .. + complex float temp; + int i,info,ix,j,jx,kx; + bool noconj,nounit; + +// Test the input parameters. +// + info = 0; + if ( !lsame( uplo, 'U' ) && !lsame( uplo, 'L' ) ) { + info = 1; + } else if ( !lsame( trans, 'N' ) && !lsame( trans, 'T' ) && + !lsame( trans, 'C' ) ) { + info = 2; + } else if ( !lsame( diag, 'U' ) && !lsame( diag, 'N' ) ) { + info = 3; + } else if (n < 0) { + info = 4; + } else if ( lda < i4_max( 1, n ) ) { + info = 6; + } else if ( incx == 0 ) { + info = 8; + } + if ( info != 0 ) { + xerbla('CTRMV', info); + return; + } +// +// Quick return if possible. +// + if ( n == 0 ) { + return; + } +// + noconj = lsame( trans, 'T' ); + nounit = lsame( diag, 'N' ); + + // Set up the start point in X if the increment is not unity. This + // will be ( N - 1 )*INCX too small for descending loops. + + if ( incx <= 0 ) { + kx = 0 - ( n-1 ) * incx; + } else if ( incx != 1 ) { + kx = 0; + } + + // Start the operations. In this version the elements of A are + // accessed sequentially with one pass through A. + + if ( lsame( trans, 'N' ) ) { + + // Form x := A*x. + + if ( lsame(uplo,'U') ) { + if (incx == 1) { + for ( j = 0; j < n; j++ ) { + if ( x[j] != zero ) { + temp = x[j]; + for ( i = 0; i < j - 1; i++ ) { + x[i] = x[i] + temp * a[i+j*lda]; + } + if ( nounit ) { + x[j] = x[j] * a[j+j*lda]; + } + } + } + } else { + jx = kx; + for ( j = 0; j < n; j++ ) { + if (x[jx] != zero) { + temp = x[jx]; + ix = kx; + for ( i = 0; i < j; i++ ) { + x[ix] = x[ix] + temp * a[i+j*lda]; + ix = ix + incx; + } + if ( nounit ) { + x[jx] = x[jx] * a[j+j*lda]; + } + } + jx = jx + incx; + } + } + } else { + if ( incx == 1 ) { + for ( j = n - 1; 0 <= j; j-- ) { + if ( x[j] != zero) { + temp = x[j]; + for ( i = n - 1; j < i; i-- ) { + x[i] = x[i] + temp * a[i+j*lda]; + } + if ( nounit ) { + x[j] = x[j] * a[j+j*lda]; + } + } + } + } else { + kx = kx + ( n - 1 ) * incx; + jx = kx; + for ( j = n - 1; 0 <= j; j-- ) { + if (x[jx] != zero) { + temp = x[jx]; + ix = kx; + for ( i = n - 1; j < i; i-- ) { + x[ix] = x[ix] + temp * a[i+j*lda]; + ix = ix - incx; + } + if ( nounit ) { + x[jx] = x[jx] * a[j+j*lda]; + } + } + jx = jx - incx; + } + } + } + } else { + + // Form x := A**T*x or x := A**H*x. + + if ( lsame(uplo,'U') ) { + if ( incx == 1 ) { + for ( j = n - 1; 0 <= j; j-- ) { + temp = x[j]; + if ( noconj ) { + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j - 1; 0 <= i; i-- ) { + temp = temp + a[i+j*lda] * x[i]; + } + } else { + if (nounit) { + temp = temp * (~a[j+j*lda]); + } + for ( i = j-1; 0 <= i; i-- ) { + temp = temp + (~a[i+j*lda])*x[i]; + } + } + x[j] = temp; + } + } else { + jx = kx + ( n-1 ) * incx; + for ( j = n - 1; 0 <= j; j-- ) { + temp = x[jx]; + ix = jx; + if ( noconj ) { + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j - 1; 0 <= i; i-- ) { + ix = ix - incx; + temp = temp + a[i+j*lda] * x[ix]; + } + } else { + if (nounit) { + temp = temp * (~a[j+j*lda]); + } + for ( i = j-1; 0 <= i; i-- ) { + ix = ix - incx; + temp = temp + (~a[i+j*lda]) * x[ix]; + } + } + x[jx] = temp; + jx = jx - incx; + } + } + } else { + if ( incx == 1 ) { + for ( j = 0; j < n; j++ ) { + temp = x[j]; + if ( noconj ) { + if (nounit) { + temp = temp*a[j+j*lda]; + } + for ( i = j + 1; i < n; i++ ) { + temp = temp + a[i+j*lda] * x[i]; + } + } else { + if (nounit) temp = temp * (~a[j+j*lda]); + for ( i = j + 1; i < n; i++ ) { + temp = temp + (~a[i+j*lda]) * x[i]; + } + } + x[j] = temp; + } + } else { + jx = kx; + for ( j = 0; j < n; j++ ) { + temp = x[jx]; + ix = jx; + if ( noconj ) { + if ( nounit ) { + temp = temp*a[j+j*lda]; + } + for ( i = j + 1; i < n; i++ ) { + ix = ix + incx; + temp = temp + a[i+j*lda] * x[ix]; + } + } else { + if (nounit) temp = temp * (~a[j+j*lda]); + for ( i = j + 1; i < n; i++ ) { + ix = ix + incx; + temp = temp + (~a[i+j*lda]) * x[ix]; + } + } + x[jx] = temp; + jx = jx + incx; + } + } + } + } + + return; + + // End of CTRMV . + +} diff --git a/src/ddot.c b/src/ddot.c index 683a04b..e94f417 100644 --- a/src/ddot.c +++ b/src/ddot.c @@ -1,3 +1,4 @@ +/* *> \brief \b DDOT * * =========== DOCUMENTATION =========== @@ -78,71 +79,63 @@ *> \endverbatim *> * ===================================================================== - DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - DOUBLE PRECISION DX(*),DY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - DOUBLE PRECISION DTEMP - INTEGER I,IX,IY,M,MP1 -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. - DDOT = 0.0d0 - DTEMP = 0.0d0 - IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* -* -* clean-up loop -* - M = MOD(N,5) - IF (M.NE.0) THEN - DO I = 1,M - DTEMP = DTEMP + DX(I)*DY(I) - END DO - IF (N.LT.5) THEN - DDOT=DTEMP - RETURN - END IF - END IF - MP1 = M + 1 - DO I = MP1,N,5 - DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + - $ DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - DTEMP = DTEMP + DX(IX)*DY(IY) - IX = IX + INCX - IY = IY + INCY - END DO - END IF - DDOT = DTEMP - RETURN -* -* End of DDOT -* - END +*/ +double ddot(int n, double dx[], int incx, double dy[], int incy) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + +// ===================================================================== +// +// .. Local Scalars .. + double dtemp; + int i, ix, iy, m; + +// .. + dtemp = 0; + if (n <= 0) return dtemp; + if (incx == 1 && incy == 1) { +// +// code for both increments equal to 1 +// +// +// clean-up loop +// + m = n % 5; + if (m != 0) { + for (i = 0; i < m; i++) { + dtemp = dtemp + dx[i]*dy[i]; + } + if (n < 5) return dtemp; + } + // mp1 = m + 1; + for (i = m; i < n; i += 5) { + dtemp += dx[i ]*dy[i ] + + dx[i+1]*dy[i+1] + + dx[i+2]*dy[i+2] + + dx[i+3]*dy[i+3] + + dx[i+4]*dy[i+4]; + } + } else { +// +// code for unequal increments or equal increments +// not equal to 1 +// + ix = 0; + iy = 0; + if (incx < 0) ix = (-n+1)*incx; + if (incy < 0) iy = (-n+1)*incy; + for (i = 0; i < n; i++) { + dtemp += dx[ix]*dy[iy]; + ix += incx; + iy += incy; + } + } + + return dtemp; + + // End DDOT + +} diff --git a/src/icamax.c b/src/icamax.c index c103cc8..db8a7ee 100644 --- a/src/icamax.c +++ b/src/icamax.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b ICAMAX * * =========== DOCUMENTATION =========== @@ -67,61 +72,60 @@ *> \endverbatim *> * ===================================================================== - INTEGER FUNCTION ICAMAX(N,CX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - COMPLEX CX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - REAL SMAX - INTEGER I,IX -* .. -* .. External Functions .. - REAL SCABS1 - EXTERNAL SCABS1 -* .. - ICAMAX = 0 - IF (N.LT.1 .OR. INCX.LE.0) RETURN - ICAMAX = 1 - IF (N.EQ.1) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - SMAX = SCABS1(CX(1)) - DO I = 2,N - IF (SCABS1(CX(I)).GT.SMAX) THEN - ICAMAX = I - SMAX = SCABS1(CX(I)) - END IF - END DO - ELSE -* -* code for increment not equal to 1 -* - IX = 1 - SMAX = SCABS1(CX(1)) - IX = IX + INCX - DO I = 2,N - IF (SCABS1(CX(IX)).GT.SMAX) THEN - ICAMAX = I - SMAX = SCABS1(CX(IX)) - END IF - IX = IX + INCX - END DO - END IF - RETURN -* -* End of ICAMAX -* - END +*/ +int icamax ( int n, float complex x[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + + // .. Local Scalars .. + float smax; + int i, ix, value; + + value = 0; + + if ( n < 1 || incx <= 0 ) { + return value; + } + + value = 1; + + if ( n == 1 ) { + return value; + } + + if ( incx != 1 ) { + ix = 0; + smax = cabs1 ( x[0] ); + ix = ix + incx; + + for ( i = 1; i < n; i++ ) + { + if ( smax < cabs1 ( x[ix] ) ) + { + value = i + 1; + smax = cabs1 ( x[ix] ); + } + ix = ix + incx; + } + } + else { + smax = cabs1 ( x[0] ); + for ( i = 1; i < n; i++ ) + { + if ( smax < cabs1 ( x[i] ) ) + { + value = i + 1; + smax = cabs1 ( x[i] ); + } + } + } + + return value; + + // End of ICAMAX + +} diff --git a/src/isamax.c b/src/isamax.c index 8ad44ad..80e699c 100644 --- a/src/isamax.c +++ b/src/isamax.c @@ -1,3 +1,6 @@ +#include "blas.h" + +/* *> \brief \b ISAMAX * * =========== DOCUMENTATION =========== @@ -67,60 +70,58 @@ *> \endverbatim *> * ===================================================================== - INTEGER FUNCTION ISAMAX(N,SX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - REAL SX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - REAL SMAX - INTEGER I,IX -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS -* .. - ISAMAX = 0 - IF (N.LT.1 .OR. INCX.LE.0) RETURN - ISAMAX = 1 - IF (N.EQ.1) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - SMAX = ABS(SX(1)) - DO I = 2,N - IF (ABS(SX(I)).GT.SMAX) THEN - ISAMAX = I - SMAX = ABS(SX(I)) - END IF - END DO - ELSE -* -* code for increment not equal to 1 -* - IX = 1 - SMAX = ABS(SX(1)) - IX = IX + INCX - DO I = 2,N - IF (ABS(SX(IX)).GT.SMAX) THEN - ISAMAX = I - SMAX = ABS(SX(IX)) - END IF - IX = IX + INCX - END DO - END IF - RETURN -* -* End of ISAMAX -* - END +*/ +int isamax ( int n, float dx[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + + +// ===================================================================== + +// .. Local Scalars .. + float dmax; + int i, ix, value; + + value = 0; + + if ( n < 1 || incx <= 0 ) { + return value; + } + + value = 1; + + if ( n == 1 ) { + return value; + } + + if ( incx == 1 ) { + dmax =fabs ( dx[0] ); + + for ( i = 1; i < n; i++ ) { + if ( dmax < fabs ( dx[i] ) ) { + value = i + 1; + dmax = fabs ( dx[i] ); + } + } + } + else { + ix = 0; + dmax = fabs ( dx[0] ); + ix = ix + incx; + + for ( i = 1; i < n; i++ ) { + if ( dmax < fabs ( dx[ix] ) ) { + value = i + 1; + dmax = fabs ( dx[ix] ); + } + ix = ix + incx; + } + } + + return value; + + // End of ISAMAX + +} \ No newline at end of file diff --git a/src/izamax.c b/src/izamax.c index e4779a1..ef027cb 100644 --- a/src/izamax.c +++ b/src/izamax.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b IZAMAX * * =========== DOCUMENTATION =========== @@ -67,61 +72,50 @@ *> \endverbatim *> * ===================================================================== - INTEGER FUNCTION IZAMAX(N,ZX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - COMPLEX*16 ZX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - DOUBLE PRECISION DMAX - INTEGER I,IX -* .. -* .. External Functions .. - DOUBLE PRECISION DCABS1 - EXTERNAL DCABS1 -* .. - IZAMAX = 0 - IF (N.LT.1 .OR. INCX.LE.0) RETURN - IZAMAX = 1 - IF (N.EQ.1) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - DMAX = DCABS1(ZX(1)) - DO I = 2,N - IF (DCABS1(ZX(I)).GT.DMAX) THEN - IZAMAX = I - DMAX = DCABS1(ZX(I)) - END IF - END DO - ELSE -* -* code for increment not equal to 1 -* - IX = 1 - DMAX = DCABS1(ZX(1)) - IX = IX + INCX - DO I = 2,N - IF (DCABS1(ZX(IX)).GT.DMAX) THEN - IZAMAX = I - DMAX = DCABS1(ZX(IX)) - END IF - IX = IX + INCX - END DO - END IF - RETURN -* -* End of IZAMAX -* - END +*/ +int izamax(int n, complex double zx[], int incx) { + +// -- Reference BLAS level1 routine (version 3.4.0) -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +// November 2011 + +// .. Local Scalars .. + double dmax; + int i, ix, ret_val; + + ret_val = 0; + if (n < 1 || incx <= 0) return ret_val; + ret_val = 1; + if (n == 1) return ret_val; + if (incx == 1) { + +// code for increment equal to 1 + dmax = dcabs1(zx[1]); + for (i = 1; i < n; i++) { + if (dcabs1(zx[i]) > dmax) { + ret_val = i; + dmax = dcabs1(zx[i]); + } + } + } else { + +// code for increment not equal to 1 + + ix = 1; + dmax = dcabs1(zx[1]); + ix += incx; + for (i = 1; i < n; i++) { + if (dcabs1(zx[ix]) > dmax) { + ret_val = i; + dmax = dcabs1(zx[ix]); + } + ix += incx; + } + } + + return ret_val; + + // End of IZAMAX + +} diff --git a/src/lsame.c b/src/lsame.c index 6aa4007..98a3b4e 100644 --- a/src/lsame.c +++ b/src/lsame.c @@ -1,3 +1,6 @@ +#include + +/* *> \brief \b LSAME * * =========== DOCUMENTATION =========== @@ -49,74 +52,33 @@ *> \ingroup aux_blas * * ===================================================================== - LOGICAL FUNCTION LSAME(CA,CB) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - CHARACTER CA,CB -* .. -* -* ===================================================================== -* -* .. Intrinsic Functions .. - INTRINSIC ICHAR -* .. -* .. Local Scalars .. - INTEGER INTA,INTB,ZCODE -* .. -* -* Test if the characters are equal -* - LSAME = CA .EQ. CB - IF (LSAME) RETURN -* -* Now test for equivalence if both characters are alphabetic. -* - ZCODE = ICHAR('Z') -* -* Use 'Z' rather than 'A' so that ASCII can be detected on Prime -* machines, on which ICHAR returns a value with bit 8 set. -* ICHAR('A') on Prime machines returns 193 which is the same as -* ICHAR('A') on an EBCDIC machine. -* - INTA = ICHAR(CA) - INTB = ICHAR(CB) -* - IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN -* -* ASCII is assumed - ZCODE is the ASCII code of either lower or -* upper case 'Z'. -* - IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32 - IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32 -* - ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN -* -* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or -* upper case 'Z'. -* - IF (INTA.GE.129 .AND. INTA.LE.137 .OR. - + INTA.GE.145 .AND. INTA.LE.153 .OR. - + INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64 - IF (INTB.GE.129 .AND. INTB.LE.137 .OR. - + INTB.GE.145 .AND. INTB.LE.153 .OR. - + INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64 -* - ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN -* -* ASCII is assumed, on Prime machines - ZCODE is the ASCII code -* plus 128 of either lower or upper case 'Z'. -* - IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32 - IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32 - END IF - LSAME = INTA .EQ. INTB -* -* RETURN -* -* End of LSAME -* - END +*/ +bool lsame ( char ca, char cb ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// Test if the characters are equal + + if ( ca == cb ) { + return true; + } + + if ( 'A' <= ca && ca <= 'Z' ) { + if ( ca - 'A' == cb - 'a' ) { + return true; + } + } else if ( 'a' <= ca && ca <= 'z' ) { + if ( ca - 'a' == cb - 'A' ) { + return true; + } + } + + return false; + + // End of LSAME + +} \ No newline at end of file diff --git a/src/sasum.c b/src/sasum.c index 6c2e0c9..7043d57 100644 --- a/src/sasum.c +++ b/src/sasum.c @@ -1,3 +1,6 @@ +#include "blas.h" + +/* *> \brief \b SASUM * * =========== DOCUMENTATION =========== @@ -68,65 +71,26 @@ *> \endverbatim *> * ===================================================================== - REAL FUNCTION SASUM(N,SX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - REAL SX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - REAL STEMP - INTEGER I,M,MP1,NINCX -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS,MOD -* .. - SASUM = 0.0e0 - STEMP = 0.0e0 - IF (N.LE.0 .OR. INCX.LE.0) RETURN - IF (INCX.EQ.1) THEN -* code for increment equal to 1 -* -* -* clean-up loop -* - M = MOD(N,6) - IF (M.NE.0) THEN - DO I = 1,M - STEMP = STEMP + ABS(SX(I)) - END DO - IF (N.LT.6) THEN - SASUM = STEMP - RETURN - END IF - END IF - MP1 = M + 1 - DO I = MP1,N,6 - STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) + - $ ABS(SX(I+2)) + ABS(SX(I+3)) + - $ ABS(SX(I+4)) + ABS(SX(I+5)) - END DO - ELSE -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO I = 1,NINCX,INCX - STEMP = STEMP + ABS(SX(I)) - END DO - END IF - SASUM = STEMP - RETURN -* -* End of SASUM -* - END +*/ +float sasum ( int n, float x[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + int i, j; + float value; + + value = 0.0; + j = 0; + + for ( i = 0; i < n; i++ ) { + value = value + fabs ( x[j] ); + j = j + incx; + } + + return value; +} \ No newline at end of file diff --git a/src/saxpy.c b/src/saxpy.c index ded238a..771bedd 100644 --- a/src/saxpy.c +++ b/src/saxpy.c @@ -1,3 +1,6 @@ +#include "blas.h" + +/* *> \brief \b SAXPY * * =========== DOCUMENTATION =========== @@ -85,68 +88,65 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL SA - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - REAL SX(*),SY(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I,IX,IY,M,MP1 -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. - IF (N.LE.0) RETURN - IF (SA.EQ.0.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN -* -* code for both increments equal to 1 -* -* -* clean-up loop -* - M = MOD(N,4) - IF (M.NE.0) THEN - DO I = 1,M - SY(I) = SY(I) + SA*SX(I) - END DO - END IF - IF (N.LT.4) RETURN - MP1 = M + 1 - DO I = MP1,N,4 - SY(I) = SY(I) + SA*SX(I) - SY(I+1) = SY(I+1) + SA*SX(I+1) - SY(I+2) = SY(I+2) + SA*SX(I+2) - SY(I+3) = SY(I+3) + SA*SX(I+3) - END DO - ELSE -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO I = 1,N - SY(IY) = SY(IY) + SA*SX(IX) - IX = IX + INCX - IY = IY + INCY - END DO - END IF - RETURN -* -* End of SAXPY -* - END +*/ +void saxpy ( int n, float da, float dx[], int incx, float dy[], int incy ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + int i, ix, iy, m; + + if ( n <= 0 ) { + return; + } + + if ( da == 0.0 ) { + return; + } + /* + Code for unequal increments or equal increments + not equal to 1. + */ + if ( incx != 1 || incy != 1 ) { + if ( 0 <= incx ) { + ix = 0; + } + else { + ix = ( - n + 1 ) * incx; + } + + if ( 0 <= incy ) { + iy = 0; + } + else { + iy = ( - n + 1 ) * incy; + } + + for ( i = 0; i < n; i++ ) { + dy[iy] = dy[iy] + da * dx[ix]; + ix = ix + incx; + iy = iy + incy; + } + } else { + // Code for both increments equal to 1. + m = n % 4; + + for ( i = 0; i < m; i++ ) { + dy[i] = dy[i] + da * dx[i]; + } + + for ( i = m; i < n; i = i + 4 ) { + dy[i ] = dy[i ] + da * dx[i ]; + dy[i+1] = dy[i+1] + da * dx[i+1]; + dy[i+2] = dy[i+2] + da * dx[i+2]; + dy[i+3] = dy[i+3] + da * dx[i+3]; + } + + } + + return; +} \ No newline at end of file diff --git a/src/scabs1.c b/src/scabs1.c index 9bacf09..197a9be 100644 --- a/src/scabs1.c +++ b/src/scabs1.c @@ -1,3 +1,8 @@ +# include + +#include "blas.h" + +/* *> \brief \b SCABS1 * * =========== DOCUMENTATION =========== @@ -42,24 +47,12 @@ *> \ingroup single_blas_level1 * * ===================================================================== - REAL FUNCTION SCABS1(Z) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX Z -* .. -* -* ===================================================================== -* -* .. Intrinsic Functions .. - INTRINSIC ABS,AIMAG,REAL -* .. - SCABS1 = ABS(REAL(Z)) + ABS(AIMAG(Z)) - RETURN -* -* End of SCABS1 -* - END +*/ +float cabs1 ( float complex z ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + + return fabs ( creal ( z ) ) + fabs ( cimag ( z ) ); +} diff --git a/src/scasum.c b/src/scasum.c index 8fd67e9..3ec5fd8 100644 --- a/src/scasum.c +++ b/src/scasum.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b SCASUM * * =========== DOCUMENTATION =========== @@ -68,50 +73,41 @@ *> \endverbatim *> * ===================================================================== - REAL FUNCTION SCASUM(N,CX,INCX) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - COMPLEX CX(*) -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - REAL STEMP - INTEGER I,NINCX -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS,AIMAG,REAL -* .. - SCASUM = 0.0e0 - STEMP = 0.0e0 - IF (N.LE.0 .OR. INCX.LE.0) RETURN - IF (INCX.EQ.1) THEN -* -* code for increment equal to 1 -* - DO I = 1,N - STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) - END DO - ELSE -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO I = 1,NINCX,INCX - STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) - END DO - END IF - SCASUM = STEMP - RETURN -* -* End of SCASUM -* - END +*/ +float scasum ( int n, float complex x[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + int i, ix; + float value; + + value = 0.0; + + if ( n <= 0 || incx <= 0 ) { + return value; + } + + if ( incx == 1 ) { + for ( i = 0; i < n; i++ ) + { + value = value + fabs ( creal ( x[i] ) ) + + fabs ( cimag ( x[i] ) ); + } + } + else { + ix = 0; + for ( i = 0; i < n; i++ ) + { + value = value + fabs ( creal ( x[ix] ) ) + + fabs ( cimag ( x[ix] ) ); + ix = ix + incx; + } + } + + return value; +} \ No newline at end of file diff --git a/src/scnrm2.c b/src/scnrm2.c index 1c2b2d4..d481ed0 100644 --- a/src/scnrm2.c +++ b/src/scnrm2.c @@ -1,209 +1,144 @@ -!> \brief \b SCNRM2 -! -! =========== DOCUMENTATION =========== -! -! Online html documentation available at -! http://www.netlib.org/lapack/explore-html/ -! -! Definition: -! =========== -! -! REAL FUNCTION SCNRM2(N,X,INCX) -! -! .. Scalar Arguments .. -! INTEGER INCX,N -! .. -! .. Array Arguments .. -! COMPLEX X(*) -! .. -! -! -!> \par Purpose: -! ============= -!> -!> \verbatim -!> -!> SCNRM2 returns the euclidean norm of a vector via the function -!> name, so that -!> -!> SCNRM2 := sqrt( x**H*x ) -!> \endverbatim -! -! Arguments: -! ========== -! -!> \param[in] N -!> \verbatim -!> N is INTEGER -!> number of elements in input vector(s) -!> \endverbatim -!> -!> \param[in] X -!> \verbatim -!> X is COMPLEX array, dimension (N) -!> complex vector with N elements -!> \endverbatim -!> -!> \param[in] INCX -!> \verbatim -!> INCX is INTEGER, storage spacing between elements of X -!> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n -!> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n -!> If INCX = 0, x isn't a vector so there is no need to call -!> this subroutine. If you call it anyway, it will count x(1) -!> in the vector norm N times. -!> \endverbatim -! -! Authors: -! ======== -! -!> \author Edward Anderson, Lockheed Martin -! -!> \date August 2016 -! -!> \ingroup single_blas_level1 -! -!> \par Contributors: -! ================== -!> -!> Weslley Pereira, University of Colorado Denver, USA -! -!> \par Further Details: -! ===================== -!> -!> \verbatim -!> -!> Anderson E. (2017) -!> Algorithm 978: Safe Scaling in the Level 1 BLAS -!> ACM Trans Math Softw 44:1--28 -!> https://doi.org/10.1145/3061665 -!> -!> Blue, James L. (1978) -!> A Portable Fortran Program to Find the Euclidean Norm of a Vector -!> ACM Trans Math Softw 4:15--23 -!> https://doi.org/10.1145/355769.355771 -!> -!> \endverbatim -!> -! ===================================================================== -function SCNRM2( n, x, incx ) - integer, parameter :: wp = kind(1.e0) - real(wp) :: SCNRM2 -! -! -- Reference BLAS level1 routine (version 3.9.1) -- -! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -! March 2021 -! -! .. Constants .. - real(wp), parameter :: zero = 0.0_wp - real(wp), parameter :: one = 1.0_wp - real(wp), parameter :: maxN = huge(0.0_wp) -! .. -! .. Blue's scaling constants .. - real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( & - (minexponent(0._wp) - 1) * 0.5_wp) - real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( & - (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp) - real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( & - (minexponent(0._wp) - digits(0._wp)) * 0.5_wp)) - real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( & - (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp)) -! .. -! .. Scalar Arguments .. - integer :: incx, n -! .. -! .. Array Arguments .. - complex(wp) :: x(*) -! .. -! .. Local Scalars .. - integer :: i, ix - logical :: notbig - real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin -! -! Quick return if possible -! - SCNRM2 = zero - if( n <= 0 ) return -! - scl = one - sumsq = zero -! -! Compute the sum of squares in 3 accumulators: -! abig -- sums of squares scaled down to avoid overflow -! asml -- sums of squares scaled up to avoid underflow -! amed -- sums of squares that do not require scaling -! The thresholds and multipliers are -! tbig -- values bigger than this are scaled down by sbig -! tsml -- values smaller than this are scaled up by ssml -! - notbig = .true. - asml = zero - amed = zero - abig = zero - ix = 1 - if( incx < 0 ) ix = 1 - (n-1)*incx - do i = 1, n - ax = abs(real(x(ix))) - if (ax > tbig) then - abig = abig + (ax*sbig)**2 - notbig = .false. - else if (ax < tsml) then - if (notbig) asml = asml + (ax*ssml)**2 - else - amed = amed + ax**2 - end if - ax = abs(aimag(x(ix))) - if (ax > tbig) then - abig = abig + (ax*sbig)**2 - notbig = .false. - else if (ax < tsml) then - if (notbig) asml = asml + (ax*ssml)**2 - else - amed = amed + ax**2 - end if - ix = ix + incx - end do -! -! Combine abig and amed or amed and asml if more than one -! accumulator was used. -! - if (abig > zero) then -! -! Combine abig and amed if abig > 0. -! - if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then - abig = abig + (amed*sbig)*sbig - end if - scl = one / sbig - sumsq = abig - else if (asml > zero) then -! -! Combine amed and asml if asml > 0. -! - if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then - amed = sqrt(amed) - asml = sqrt(asml) / ssml - if (asml > amed) then - ymin = amed - ymax = asml - else - ymin = asml - ymax = amed - end if - scl = one - sumsq = ymax**2*( one + (ymin/ymax)**2 ) - else - scl = one / ssml - sumsq = asml - end if - else -! -! Otherwise all values are mid-range -! - scl = one - sumsq = amed - end if - SCNRM2 = scl*sqrt( sumsq ) - return -end function +#include + +#include "blas.h" + +/* +*> \brief \b SCNRM2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* REAL FUNCTION SCNRM2(N,X,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* COMPLEX X(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SCNRM2 returns the euclidean norm of a vector via the function +*> name, so that +*> +*> SCNRM2 := sqrt( x**H*x ) +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> number of elements in input vector(s) +*> \endverbatim +*> +*> \param[in] X +*> \verbatim +*> X is COMPLEX array, dimension (N) +*> complex vector with N elements +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER, storage spacing between elements of X +*> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n +*> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n +*> If INCX = 0, x isn't a vector so there is no need to call +*> this subroutine. If you call it anyway, it will count x(1) +*> in the vector norm N times. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Edward Anderson, Lockheed Martin +* +*> \date August 2016 +* +*> \ingroup single_blas_level1 +* +*> \par Contributors: +* ================== +*> +*> Weslley Pereira, University of Colorado Denver, USA +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Anderson E. (2017) +*> Algorithm 978: Safe Scaling in the Level 1 BLAS +*> ACM Trans Math Softw 44:1--28 +*> https://doi.org/10.1145/3061665 +*> +*> Blue, James L. (1978) +*> A Portable Fortran Program to Find the Euclidean Norm of a Vector +*> ACM Trans Math Softw 4:15--23 +*> https://doi.org/10.1145/355769.355771 +*> +*> \endverbatim +*> +* ===================================================================== +*/ +float scnrm2 ( int n, float complex x[], int incx ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + int i, ix; + float scale, ssq, temp, value; + + if ( n < 1 || incx < 1 ) { + value = 0.0; + } else { + scale = 0.0; + ssq = 1.0; + ix = 0; + + for ( i = 0; i < n; i++ ) { + if ( creal ( x[ix] ) != 0.0 ) { + temp = fabs ( creal ( x[ix] ) ); + if ( scale < temp ) { + ssq = 1.0 + ssq * pow ( scale / temp, 2 ); + scale = temp; + } else { + ssq = ssq + pow ( temp / scale, 2 ); + } + } + + if ( cimag ( x[ix] ) != 0.0 ) { + temp = fabs ( cimag ( x[ix] ) ); + if ( scale < temp ) + { + ssq = 1.0 + ssq * pow ( scale / temp, 2 ); + scale = temp; + } + else + { + ssq = ssq + pow ( temp / scale, 2 ); + } + } + ix = ix + incx; + } + value = scale * sqrt ( ssq ); + } + + return value; + +} \ No newline at end of file diff --git a/src/sgemv.c b/src/sgemv.c index 0517b12..3d82e44 100644 --- a/src/sgemv.c +++ b/src/sgemv.c @@ -1,3 +1,6 @@ +# include "blas0.h" + +/* *> \brief \b SGEMV * * =========== DOCUMENTATION =========== @@ -152,176 +155,164 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL ALPHA,BETA - INTEGER INCX,INCY,LDA,M,N - CHARACTER TRANS -* .. -* .. Array Arguments .. - REAL A(LDA,*),X(*),Y(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ONE,ZERO - PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 1 - ELSE IF (M.LT.0) THEN - INFO = 2 - ELSE IF (N.LT.0) THEN - INFO = 3 - ELSE IF (LDA.LT.MAX(1,M)) THEN - INFO = 6 - ELSE IF (INCX.EQ.0) THEN - INFO = 8 - ELSE IF (INCY.EQ.0) THEN - INFO = 11 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('SGEMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. - + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN -* -* Set LENX and LENY, the lengths of the vectors x and y, and set -* up the start points in X and Y. -* - IF (LSAME(TRANS,'N')) THEN - LENX = N - LENY = M - ELSE - LENX = M - LENY = N - END IF - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (LENX-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (LENY-1)*INCY - END IF -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* -* First form y := beta*y. -* - IF (BETA.NE.ONE) THEN - IF (INCY.EQ.1) THEN - IF (BETA.EQ.ZERO) THEN - DO 10 I = 1,LENY - Y(I) = ZERO - 10 CONTINUE - ELSE - DO 20 I = 1,LENY - Y(I) = BETA*Y(I) - 20 CONTINUE - END IF - ELSE - IY = KY - IF (BETA.EQ.ZERO) THEN - DO 30 I = 1,LENY - Y(IY) = ZERO - IY = IY + INCY - 30 CONTINUE - ELSE - DO 40 I = 1,LENY - Y(IY) = BETA*Y(IY) - IY = IY + INCY - 40 CONTINUE - END IF - END IF - END IF - IF (ALPHA.EQ.ZERO) RETURN - IF (LSAME(TRANS,'N')) THEN -* -* Form y := alpha*A*x + y. -* - JX = KX - IF (INCY.EQ.1) THEN - DO 60 J = 1,N - TEMP = ALPHA*X(JX) - DO 50 I = 1,M - Y(I) = Y(I) + TEMP*A(I,J) - 50 CONTINUE - JX = JX + INCX - 60 CONTINUE - ELSE - DO 80 J = 1,N - TEMP = ALPHA*X(JX) - IY = KY - DO 70 I = 1,M - Y(IY) = Y(IY) + TEMP*A(I,J) - IY = IY + INCY - 70 CONTINUE - JX = JX + INCX - 80 CONTINUE - END IF - ELSE -* -* Form y := alpha*A**T*x + y. -* - JY = KY - IF (INCX.EQ.1) THEN - DO 100 J = 1,N - TEMP = ZERO - DO 90 I = 1,M - TEMP = TEMP + A(I,J)*X(I) - 90 CONTINUE - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 100 CONTINUE - ELSE - DO 120 J = 1,N - TEMP = ZERO - IX = KX - DO 110 I = 1,M - TEMP = TEMP + A(I,J)*X(IX) - IX = IX + INCX - 110 CONTINUE - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 120 CONTINUE - END IF - END IF -* - RETURN -* -* End of SGEMV -* - END +*/ +void sgemv ( char trans, int m, int n, float alpha, float a[], int lda, + float x[], int incx, float beta, float y[], int incy ) { + +// -- Reference BLAS level2 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + float temp; + int i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny; + + // Test the input parameters. + + info = 0; + if ( ! lsame ( trans, 'N' ) && + ! lsame ( trans, 'T' ) && + ! lsame ( trans, 'C' ) ) { + info = 1; + } else if ( m < 0 ) { + info = 2; + } else if ( n < 0 ) { + info = 3; + } else if ( lda < i4_max ( 1, m ) ) { + info = 6; + } else if ( incx == 0 ) { + info = 8; + } else if ( incy == 0 ) { + info = 11; + } + + if ( info != 0 ) { + xerbla ( "SGEMV", info ); + return; + } + + // Quick return if possible. + + if ( ( m == 0 ) || + ( n == 0 ) || + ( ( alpha == 0.0 ) && ( beta == 1.0 ) ) ) { + return; + } + /* + Set LENX and LENY, the lengths of the vectors x and y, and set + up the start points in X and Y. + */ + if ( lsame ( trans, 'N' ) ) { + lenx = n; + leny = m; + } else { + lenx = m; + leny = n; + } + + if ( 0 < incx ) { + kx = 0; + } else { + kx = 0 - ( lenx - 1 ) * incx; + } + + if ( 0 < incy ) { + ky = 0; + } else { + ky = 0 - ( leny - 1 ) * incy; + } + /* + Start the operations. In this version the elements of A are + accessed sequentially with one pass through A. + + First form y := beta*y. + */ + if ( beta != 1.0 ) { + if ( incy == 1 ) { + if ( beta == 0.0 ) { + for ( i = 0; i < leny; i++ ) { + y[i] = 0.0; + } + } else { + for ( i = 0; i < leny; i++ ) { + y[i] = beta * y[i]; + } + } + } else { + iy = ky; + if ( beta == 0.0 ) { + for ( i = 0; i < leny; i++ ) { + y[iy] = 0.0; + iy = iy + incy; + } + } else { + for ( i = 0; i < leny; i++ ) { + y[iy] = beta * y[iy]; + iy = iy + incy; + } + } + } + } + + if ( alpha == 0.0 ) { + return; + } + + // Form y := alpha*A*x + y. + + if ( lsame ( trans, 'N' ) ) { + jx = kx; + if ( incy == 1 ) { + for ( j = 0; j < n; j++ ) { + if ( x[jx] != 0.0 ) { + temp = alpha * x[jx]; + for ( i = 0; i < m; i++ ) { + y[i] = y[i] + temp * a[i+j*lda]; + } + } + jx = jx + incx; + } + } else { + for ( j = 0; j < n; j++ ) { + if ( x[jx] != 0.0 ) { + temp = alpha * x[jx]; + iy = ky; + for ( i = 0; i < m; i++ ) + { + y[iy] = y[iy] + temp * a[i+j*lda]; + iy = iy + incy; + } + } + jx = jx + incx; + } + } + } else { + // Form y := alpha*A'*x + y. + jy = ky; + if ( incx == 1 ) { + for ( j = 0; j < n; j++ ) { + temp = 0.0; + for ( i = 0; i < m; i++ ) { + temp = temp + a[i+j*lda] * x[i]; + } + y[jy] = y[jy] + alpha * temp; + jy = jy + incy; + } + } else { + for ( j = 0; j < n; j++ ) { + temp = 0.0; + ix = kx; + for ( i = 0; i < m; i++ ) { + temp = temp + a[i+j*lda] * x[ix]; + ix = ix + incx; + } + y[jy] = y[jy] + alpha * temp; + jy = jy + incy; + } + } + } + + return; +} \ No newline at end of file diff --git a/src/sger.c b/src/sger.c index 9dfe4a2..c1d8a06 100644 --- a/src/sger.c +++ b/src/sger.c @@ -1,3 +1,6 @@ +#include "blas.h" + +/* *> \brief \b SGER * * =========== DOCUMENTATION =========== @@ -126,99 +129,85 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL ALPHA - INTEGER INCX,INCY,LDA,M,N -* .. -* .. Array Arguments .. - REAL A(LDA,*),X(*),Y(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,J,JY,KX -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (M.LT.0) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - ELSE IF (LDA.LT.MAX(1,M)) THEN - INFO = 9 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('SGER ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* - IF (INCY.GT.0) THEN - JY = 1 - ELSE - JY = 1 - (N-1)*INCY - END IF - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*Y(JY) - DO 10 I = 1,M - A(I,J) = A(I,J) + X(I)*TEMP - 10 CONTINUE - END IF - JY = JY + INCY - 20 CONTINUE - ELSE - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (M-1)*INCX - END IF - DO 40 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*Y(JY) - IX = KX - DO 30 I = 1,M - A(I,J) = A(I,J) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - END IF - JY = JY + INCY - 40 CONTINUE - END IF -* - RETURN -* -* End of SGER -* - END +*/ +void sger ( int m, int n, float alpha, float x[], int incx, float y[], + int incy, float a[], int lda ) { + +// -- Reference BLAS level2 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + float temp; + int i, info, ix, j, jy, kx; + + // Test the input parameters. + + info = 0; + if ( m < 0 ) { + info = 1; + } else if ( n < 0 ) { + info = 2; + } else if ( incx == 0 ) { + info = 5; + } else if ( incy == 0 ) { + info = 7; + } else if ( lda < i4_max ( 1, m ) ) { + info = 9; + } + + if ( info != 0 ) { + xerbla ( "SGER", info ); + return; + } + + // Quick return if possible. + + if ( m == 0 || n == 0 || alpha == 0.0 ) { + return; + } + /* + Start the operations. In this version the elements of A are + accessed sequentially with one pass through A. + */ + if ( 0 < incy ) { + jy = 0; + } else { + jy = 0 - ( n - 1 ) * incy; + } + + if ( incx == 1 ) { + for ( j = 0; j < n; j++ ) { + if ( y[jy] != 0.0 ) { + temp = alpha * y[jy]; + for ( i = 0; i < m; i++ ) { + a[i+j*lda] = a[i+j*lda] + x[i] * temp; + } + } + jy = jy + incy; + } + } else { + if ( 0 < incx ) { + kx = 0; + } + else { + kx = 0 - ( m - 1 ) * incx; + } + for ( j = 0; j < n; j++ ) { + if ( y[jy] != 0.0 ) { + temp = alpha * y[jy]; + ix = kx; + for ( i = 0; i < m; i++ ) { + a[i+j*lda] = a[i+j*lda] + x[ix] * temp; + ix = ix + incx; + } + } + jy = jy + incy; + } + } + + return; + +} \ No newline at end of file diff --git a/src/strmv.c b/src/strmv.c index 6f5a3e8..06802bb 100644 --- a/src/strmv.c +++ b/src/strmv.c @@ -1,3 +1,6 @@ +#include "blas.h" + +/* *> \brief \b STRMV * * =========== DOCUMENTATION =========== @@ -143,197 +146,190 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE STRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) -* -* -- Reference BLAS level2 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX,LDA,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - REAL A(LDA,*),X(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,J,JX,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (LDA.LT.MAX(1,N)) THEN - INFO = 6 - ELSE IF (INCX.EQ.0) THEN - INFO = 8 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('STRMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := A*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*A(I,J) - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*A(J,J) - END IF - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 I = 1,J - 1 - X(IX) = X(IX) + TEMP*A(I,J) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*A(J,J) - END IF - JX = JX + INCX - 40 CONTINUE - END IF - ELSE - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*A(I,J) - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*A(J,J) - END IF - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 I = N,J + 1,-1 - X(IX) = X(IX) + TEMP*A(I,J) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*A(J,J) - END IF - JX = JX - INCX - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A**T*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 100 J = N,1,-1 - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 90 I = J - 1,1,-1 - TEMP = TEMP + A(I,J)*X(I) - 90 CONTINUE - X(J) = TEMP - 100 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 120 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 110 I = J - 1,1,-1 - IX = IX - INCX - TEMP = TEMP + A(I,J)*X(IX) - 110 CONTINUE - X(JX) = TEMP - JX = JX - INCX - 120 CONTINUE - END IF - ELSE - IF (INCX.EQ.1) THEN - DO 140 J = 1,N - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 130 I = J + 1,N - TEMP = TEMP + A(I,J)*X(I) - 130 CONTINUE - X(J) = TEMP - 140 CONTINUE - ELSE - JX = KX - DO 160 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*A(J,J) - DO 150 I = J + 1,N - IX = IX + INCX - TEMP = TEMP + A(I,J)*X(IX) - 150 CONTINUE - X(JX) = TEMP - JX = JX + INCX - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of STRMV -* - END +*/ +void strmv(char uplo, char trans, char diag, int n, float a[], int lda, + float x[], int incx) { +// -- Reference BLAS level2 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + + int i, info, ix, j, jx, kx, nounit; + float temp; + + // Test the input parameters. + + info = 0; + if ( ! lsame ( uplo, 'U' ) && ! lsame ( uplo, 'L' ) ) { + info = 1; + } else if ( ! lsame ( trans, 'N' ) && ! lsame ( trans, 'T' ) && + ! lsame ( trans, 'C' ) ) { + info = 2; + } else if ( ! lsame ( diag, 'U' ) && ! lsame ( diag, 'N' ) ) { + info = 3; + } else if ( n < 0 ) { + info = 4; + } else if ( lda < i4_max ( 1, n ) ) { + info = 6; + } else if ( incx == 0 ) { + info = 8; + } + + if ( info != 0 ) { + xerbla ( "STRMV", info ); + return; + } + + // Quick return if possible. + + if ( n == 0 ) { + return; + } + + nounit = lsame ( diag, 'N' ); + /* + Set up the start point in X if the increment is not unity. This + will be ( N - 1 ) * INCX too small for descending loops. + */ + if ( incx <= 0 ) { + kx = 0 - ( n - 1 ) * incx; + } else if ( incx != 1 ) { + kx = 0; + } + /* + Start the operations. In this version the elements of A are + accessed sequentially with one pass through A. + */ + if ( lsame ( trans, 'N' ) ) { + + // Form x := A*x. + + if ( lsame ( uplo, 'U' ) ) { + if ( incx == 1 ) { + for ( j = 0; j < n; j++ ) { + if ( x[j] != 0.0 ) { + temp = x[j]; + for ( i = 0; i < j; i++ ) { + x[i] = x[i] + temp * a[i+j*lda]; + } + if ( nounit ) { + x[j] = x[j] * a[j+j*lda]; + } + } + } + } else { + jx = kx; + for ( j = 0; j < n; j++ ) { + if ( x[jx] != 0.0 ) { + temp = x[jx]; + ix = kx; + for ( i = 0; i < j; i++ ) { + x[ix] = x[ix] + temp * a[i+j*lda]; + ix = ix + incx; + } + if ( nounit ) { + x[jx] = x[jx] * a[j+j*lda]; + } + } + jx = jx + incx; + } + } + } else { + if ( incx == 1 ) { + for ( j = n - 1; 0 <= j; j-- ) { + if ( x[j] != 0.0 ) { + temp = x[j]; + for ( i = n - 1; j < i; i-- ) { + x[i] = x[i] + temp * a[i+j*lda]; + } + if ( nounit ) { + x[j] = x[j] * a[j+j*lda]; + } + } + } + } else { + kx = kx + ( n - 1 ) * incx; + jx = kx; + for ( j = n - 1; 0 <= j; j-- ) { + if ( x[jx] != 0.0 ) { + temp = x[jx]; + ix = kx; + for ( i = n - 1; j < i; i-- ) { + x[ix] = x[ix] + temp * a[i+j*lda]; + ix = ix - incx; + } + if ( nounit ) { + x[jx] = x[jx] * a[j+j*lda]; + } + } + jx = jx - incx; + } + } + } + } else { + // Form x := A'*x. + if ( lsame ( uplo, 'U' ) ) { + if ( incx == 1 ) { + for ( j = n - 1; 0 <= j; j-- ) { + temp = x[j]; + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j - 1; 0 <= i; i-- ) { + temp = temp + a[i+j*lda] * x[i]; + } + x[j] = temp; + } + } else { + jx = kx + ( n - 1 ) * incx; + for ( j = n - 1; 0 <= j; j-- ) { + temp = x[jx]; + ix = jx; + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j - 1; 0 <= i; i-- ) { + ix = ix - incx; + temp = temp + a[i+j*lda] * x[ix]; + } + x[jx] = temp; + jx = jx - incx; + } + } + } else { + if ( incx == 1 ) { + for ( j = 0; j < n; j++ ) { + temp = x[j]; + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j + 1; i < n; i++ ) { + temp = temp + a[i+j*lda] * x[i]; + } + x[j] = temp; + } + } else { + jx = kx; + for ( j = 0; j < n; j++ ) { + temp = x[jx]; + ix = jx; + if ( nounit ) { + temp = temp * a[j+j*lda]; + } + for ( i = j + 1; i < n; i++ ) { + ix = ix + incx; + temp = temp + a[i+j*lda] * x[ix]; + } + x[jx] = temp; + jx = jx + incx; + } + } + } + } + + return; + + // End of STRMV + +} \ No newline at end of file diff --git a/src/strsm.c b/src/strsm.c index 6e97c52..34e8644 100644 --- a/src/strsm.c +++ b/src/strsm.c @@ -1,3 +1,8 @@ +#include + +#include "blas.h" + +/* *> \brief \b STRSM * * =========== DOCUMENTATION =========== @@ -177,264 +182,263 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE STRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) -* -* -- Reference BLAS level3 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL ALPHA - INTEGER LDA,LDB,M,N - CHARACTER DIAG,SIDE,TRANSA,UPLO -* .. -* .. Array Arguments .. - REAL A(LDA,*),B(LDB,*) -* .. -* -* ===================================================================== -* -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,J,K,NROWA - LOGICAL LSIDE,NOUNIT,UPPER -* .. -* .. Parameters .. - REAL ONE,ZERO - PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) -* .. -* -* Test the input parameters. -* - LSIDE = LSAME(SIDE,'L') - IF (LSIDE) THEN - NROWA = M - ELSE - NROWA = N - END IF - NOUNIT = LSAME(DIAG,'N') - UPPER = LSAME(UPLO,'U') -* - INFO = 0 - IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN - INFO = 1 - ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN - INFO = 2 - ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. - + (.NOT.LSAME(TRANSA,'T')) .AND. - + (.NOT.LSAME(TRANSA,'C'))) THEN - INFO = 3 - ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN - INFO = 4 - ELSE IF (M.LT.0) THEN - INFO = 5 - ELSE IF (N.LT.0) THEN - INFO = 6 - ELSE IF (LDA.LT.MAX(1,NROWA)) THEN - INFO = 9 - ELSE IF (LDB.LT.MAX(1,M)) THEN - INFO = 11 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('STRSM ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (M.EQ.0 .OR. N.EQ.0) RETURN -* -* And when alpha.eq.zero. -* - IF (ALPHA.EQ.ZERO) THEN - DO 20 J = 1,N - DO 10 I = 1,M - B(I,J) = ZERO - 10 CONTINUE - 20 CONTINUE - RETURN - END IF -* -* Start the operations. -* - IF (LSIDE) THEN - IF (LSAME(TRANSA,'N')) THEN -* -* Form B := alpha*inv( A )*B. -* - IF (UPPER) THEN - DO 60 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 30 I = 1,M - B(I,J) = ALPHA*B(I,J) - 30 CONTINUE - END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF - 50 CONTINUE - 60 CONTINUE - ELSE - DO 100 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 70 I = 1,M - B(I,J) = ALPHA*B(I,J) - 70 CONTINUE - END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF - 90 CONTINUE - 100 CONTINUE - END IF - ELSE -* -* Form B := alpha*inv( A**T )*B. -* - IF (UPPER) THEN - DO 130 J = 1,N - DO 120 I = 1,M - TEMP = ALPHA*B(I,J) - DO 110 K = 1,I - 1 - TEMP = TEMP - A(K,I)*B(K,J) - 110 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(I,I) - B(I,J) = TEMP - 120 CONTINUE - 130 CONTINUE - ELSE - DO 160 J = 1,N - DO 150 I = M,1,-1 - TEMP = ALPHA*B(I,J) - DO 140 K = I + 1,M - TEMP = TEMP - A(K,I)*B(K,J) - 140 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(I,I) - B(I,J) = TEMP - 150 CONTINUE - 160 CONTINUE - END IF - END IF - ELSE - IF (LSAME(TRANSA,'N')) THEN -* -* Form B := alpha*B*inv( A ). -* - IF (UPPER) THEN - DO 210 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 170 I = 1,M - B(I,J) = ALPHA*B(I,J) - 170 CONTINUE - END IF - DO 190 K = 1,J - 1 - IF (A(K,J).NE.ZERO) THEN - DO 180 I = 1,M - B(I,J) = B(I,J) - A(K,J)*B(I,K) - 180 CONTINUE - END IF - 190 CONTINUE - IF (NOUNIT) THEN - TEMP = ONE/A(J,J) - DO 200 I = 1,M - B(I,J) = TEMP*B(I,J) - 200 CONTINUE - END IF - 210 CONTINUE - ELSE - DO 260 J = N,1,-1 - IF (ALPHA.NE.ONE) THEN - DO 220 I = 1,M - B(I,J) = ALPHA*B(I,J) - 220 CONTINUE - END IF - DO 240 K = J + 1,N - IF (A(K,J).NE.ZERO) THEN - DO 230 I = 1,M - B(I,J) = B(I,J) - A(K,J)*B(I,K) - 230 CONTINUE - END IF - 240 CONTINUE - IF (NOUNIT) THEN - TEMP = ONE/A(J,J) - DO 250 I = 1,M - B(I,J) = TEMP*B(I,J) - 250 CONTINUE - END IF - 260 CONTINUE - END IF - ELSE -* -* Form B := alpha*B*inv( A**T ). -* - IF (UPPER) THEN - DO 310 K = N,1,-1 - IF (NOUNIT) THEN - TEMP = ONE/A(K,K) - DO 270 I = 1,M - B(I,K) = TEMP*B(I,K) - 270 CONTINUE - END IF - DO 290 J = 1,K - 1 - IF (A(J,K).NE.ZERO) THEN - TEMP = A(J,K) - DO 280 I = 1,M - B(I,J) = B(I,J) - TEMP*B(I,K) - 280 CONTINUE - END IF - 290 CONTINUE - IF (ALPHA.NE.ONE) THEN - DO 300 I = 1,M - B(I,K) = ALPHA*B(I,K) - 300 CONTINUE - END IF - 310 CONTINUE - ELSE - DO 360 K = 1,N - IF (NOUNIT) THEN - TEMP = ONE/A(K,K) - DO 320 I = 1,M - B(I,K) = TEMP*B(I,K) - 320 CONTINUE - END IF - DO 340 J = K + 1,N - IF (A(J,K).NE.ZERO) THEN - TEMP = A(J,K) - DO 330 I = 1,M - B(I,J) = B(I,J) - TEMP*B(I,K) - 330 CONTINUE - END IF - 340 CONTINUE - IF (ALPHA.NE.ONE) THEN - DO 350 I = 1,M - B(I,K) = ALPHA*B(I,K) - 350 CONTINUE - END IF - 360 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of STRSM -* - END +*/ +void strsm ( char side, char uplo, char transa, char diag, int m, int n, + float alpha, float a[], int lda, float b[], int ldb ) { + +// -- Reference BLAS level3 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + +// ===================================================================== + +// .. Local Scalars .. + float temp; + int i, info, j, k, nrowa; + bool lside, nounit, upper; + + // Test the input parameters. + + lside = lsame ( side, 'L' ); + + if ( lside ) { + nrowa = m; + } + else { + nrowa = n; + } + + nounit = lsame ( diag, 'N' ); + upper = lsame ( uplo, 'U' ); + + info = 0; + + if ( ( ! lside ) && ( ! lsame ( side, 'R' ) ) ) { + info = 1; + } + else if ( ( ! upper ) && ( ! lsame ( uplo, 'L' ) ) ) { + info = 2; + } + else if ( ( ! lsame ( transa, 'N' ) ) && + ( ! lsame ( transa, 'T' ) ) && + ( ! lsame ( transa, 'C' ) ) ) { + info = 3; + } + else if ( ( ! lsame ( diag, 'U' ) ) && ( ! lsame ( diag, 'N' ) ) ) { + info = 4; + } + else if ( m < 0 ) { + info = 5; + } + else if ( n < 0 ) { + info = 6; + } + else if ( lda < i4_max ( 1, nrowa ) ) { + info = 9; + } + else if ( ldb < i4_max ( 1, m ) ) { + info = 11; + } + + if ( info != 0 ) { + xerbla ( "STRSM", info ); + return; + } + + // Quick return if possible. + + if ( n == 0 ) { + return; + } + + // and when alpha is 0.0. + + if ( alpha == 0.0 ) { + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = 0.0; + } + } + return; + } + + // Start the operations. + + if ( lside ) { + + // Form B := alpha*inv( a )*B. + + if ( lsame ( transa, 'N' ) ) { + if ( upper ) { + for ( j = 0; j < n; j++ ) { + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = alpha * b[i+j*ldb]; + } + } + for ( k = m - 1; 0 <= k; k-- ) { + if ( b[k+j*ldb] != 0.0 ) { + if ( nounit ) { + b[k+j*ldb] = b[k+j*ldb] / a[k+k*lda]; + } + for ( i = 0; i < k; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - b[k+j*ldb] * a[i+k*lda]; + } + } + } + } + } else { + for ( j = 0; j < n; j++ ) { + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = alpha * b[i+j*ldb]; + } + } + for ( k = 0; k < m; k++ ) { + if ( b[k+j*ldb] != 0.0 ) { + if ( nounit ) { + b[k+j*ldb] = b[k+j*ldb] / a[k+k*lda]; + } + for ( i = k + 1; i < m; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - b[k+j*ldb] * a[i+k*lda]; + } + } + } + } + } + } else { + // Form B := alpha*inv( A' )*B. + if ( upper ) { + for ( j = 0; j < n; j++ ) { + for ( i = 0; i < m; i++ ) { + temp = alpha * b[i+j*ldb]; + for ( k = 0; k < i; k++ ) { + temp = temp - a[k+i*lda] * b[k+j*ldb]; + } + if ( nounit ) { + temp = temp / a[i+i*lda]; + } + b[i+j*ldb] = temp; + } + } + } else { + for ( j = 0; j < n; j++ ) { + for ( i = m - 1; 0 <= i; i-- ) { + temp = alpha * b[i+j*ldb]; + for ( k = i + 1; k < m; k++ ) { + temp = temp - a[k+i*lda] * b[k+j*ldb]; + } + if ( nounit ) { + temp = temp / a[i+i*lda]; + } + b[i+j*ldb] = temp; + } + } + } + } + } else { + // Form B := alpha*B*inv( A ). + if ( lsame ( transa, 'N' ) ) { + if ( upper ) { + for ( j = 0; j < n; j++ ) { + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = alpha * b[i+j*ldb]; + } + } + for ( k = 0; k < j; k++ ) { + if ( a[k+j*lda] != 0.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - a[k+j*lda] * b[i+k*ldb]; + } + } + } + if ( nounit ) { + temp = 1.0 / a[j+j*lda]; + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = temp * b[i+j*ldb]; + } + } + } + } + else { + for ( j = n - 1; 0 <= j; j-- ) { + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = alpha * b[i+j*ldb]; + } + } + for ( k = j + 1; k < n; k++ ) { + if ( a[k+j*lda] != 0.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - a[k+j*lda] * b[i+k*ldb]; + } + } + } + if ( nounit ) { + temp = 1.0 / a[j+j*lda]; + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = temp * b[i+j*ldb]; + } + } + } + } + } + + // Form B := alpha*B*inv( A' ). + + else { + if ( upper ) { + for ( k = n - 1; 0 <= k; k-- ) { + if ( nounit ) { + temp = 1.0 / a[k+k*lda]; + for ( i = 0; i < m; i++ ) { + b[i+k*ldb] = temp * b[i+k*ldb]; + } + } + for ( j = 0; j < k; j++ ) { + if ( a[j+k*lda] != 0.0 ) { + temp = a[j+k*lda]; + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - temp * b[i+k*ldb]; + } + } + } + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+k*ldb] = alpha * b[i+k*ldb]; + } + } + } + } + else { + for ( k = 0; k < n; k++ ) { + if ( nounit ) { + temp = 1.0 / a[k+k*lda]; + for ( i = 0; i < m; i++ ) { + b[i+k*ldb] = temp * b[i+k*ldb]; + } + } + for ( j = k + 1; j < n; j++ ) { + if ( a[j+k*lda] != 0.0 ) { + temp = a[j+k*lda]; + for ( i = 0; i < m; i++ ) { + b[i+j*ldb] = b[i+j*ldb] - temp * b[i+k*ldb]; + } + } + } + if ( alpha != 1.0 ) { + for ( i = 0; i < m; i++ ) { + b[i+k*ldb] = alpha * b[i+k*ldb]; + } + } + } + } + } + } + + return; + + // End of STRSM + +} \ No newline at end of file diff --git a/src/xerbla.c b/src/xerbla.c index 0d73705..1cdd597 100644 --- a/src/xerbla.c +++ b/src/xerbla.c @@ -1,3 +1,8 @@ +#include + +#include "errquit.h" + +/* *> \brief \b XERBLA * * =========== DOCUMENTATION =========== @@ -56,31 +61,17 @@ *> \ingroup aux_blas * * ===================================================================== - SUBROUTINE XERBLA( SRNAME, INFO ) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - CHARACTER*(*) SRNAME - INTEGER INFO -* .. -* -* ===================================================================== -* -* .. Intrinsic Functions .. - INTRINSIC LEN_TRIM -* .. -* .. Executable Statements .. -* - WRITE( *, FMT = 9999 )SRNAME( 1:LEN_TRIM( SRNAME ) ), INFO -* - STOP -* - 9999 FORMAT( ' ** On entry to ', A, ' parameter number ', I2, ' had ', - $ 'an illegal value' ) -* -* End of XERBLA -* - END +*/ +void xerbla( char *srname, int info ) { + +// -- Reference BLAS level1 routine -- +// -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- + + printf("xerbla:double: lapack error\n"); + printf(" ** On entry to %s parameter number %d had an illegal value", srname, info); + exit(1); + + // End of XERBLA + +} \ No newline at end of file