diff --git a/internal/testdata/netlib/daxpy.f b/internal/testdata/netlib/daxpy.f new file mode 100644 index 00000000..64a02d68 --- /dev/null +++ b/internal/testdata/netlib/daxpy.f @@ -0,0 +1,115 @@ +*> \brief \b DAXPY +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION DA +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DAXPY constant times a vector plus a vector. +*> uses unrolled loops for increments equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DAXPY(N,DA,DX,INCX,DY,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 +* +* .. Scalar Arguments .. + DOUBLE PRECISION DA + INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. + DOUBLE PRECISION DX(*),DY(*) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I,IX,IY,M,MP1 +* .. +* .. Intrinsic Functions .. + INTRINSIC MOD +* .. + IF (N.LE.0) RETURN + IF (DA.EQ.0.0d0) 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 + DY(I) = DY(I) + DA*DX(I) + END DO + END IF + IF (N.LT.4) RETURN + MP1 = M + 1 + DO I = MP1,N,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) + 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 + DY(IY) = DY(IY) + DA*DX(IX) + IX = IX + INCX + IY = IY + INCY + END DO + END IF + RETURN + END diff --git a/internal/testdata/netlib/dcopy.f b/internal/testdata/netlib/dcopy.f new file mode 100644 index 00000000..d9d5ac7a --- /dev/null +++ b/internal/testdata/netlib/dcopy.f @@ -0,0 +1,115 @@ +*> \brief \b DCOPY +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DCOPY copies a vector, x, to a vector, y. +*> uses unrolled loops for increments equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DCOPY(N,DX,INCX,DY,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 +* +* .. Scalar Arguments .. + INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. + DOUBLE PRECISION DX(*),DY(*) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I,IX,IY,M,MP1 +* .. +* .. Intrinsic Functions .. + INTRINSIC MOD +* .. + 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,7) + IF (M.NE.0) THEN + DO I = 1,M + DY(I) = DX(I) + END DO + IF (N.LT.7) RETURN + END IF + MP1 = M + 1 + DO I = MP1,N,7 + DY(I) = DX(I) + DY(I+1) = DX(I+1) + DY(I+2) = DX(I+2) + DY(I+3) = DX(I+3) + DY(I+4) = DX(I+4) + DY(I+5) = DX(I+5) + DY(I+6) = DX(I+6) + 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 + DY(IY) = DX(IX) + IX = IX + INCX + IY = IY + INCY + END DO + END IF + RETURN + END diff --git a/internal/testdata/netlib/dgemm.f b/internal/testdata/netlib/dgemm.f new file mode 100644 index 00000000..4bae243a --- /dev/null +++ b/internal/testdata/netlib/dgemm.f @@ -0,0 +1,384 @@ +*> \brief \b DGEMM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA,BETA +* INTEGER K,LDA,LDB,LDC,M,N +* CHARACTER TRANSA,TRANSB +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGEMM performs one of the matrix-matrix operations +*> +*> C := alpha*op( A )*op( B ) + beta*C, +*> +*> where op( X ) is one of +*> +*> op( X ) = X or op( X ) = X**T, +*> +*> alpha and beta are scalars, and A, B and C are matrices, with op( A ) +*> an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is CHARACTER*1 +*> On entry, TRANSA specifies the form of op( A ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSA = 'N' or 'n', op( A ) = A. +*> +*> TRANSA = 'T' or 't', op( A ) = A**T. +*> +*> TRANSA = 'C' or 'c', op( A ) = A**T. +*> \endverbatim +*> +*> \param[in] TRANSB +*> \verbatim +*> TRANSB is CHARACTER*1 +*> On entry, TRANSB specifies the form of op( B ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSB = 'N' or 'n', op( B ) = B. +*> +*> TRANSB = 'T' or 't', op( B ) = B**T. +*> +*> TRANSB = 'C' or 'c', op( B ) = B**T. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of the matrix +*> op( A ) and of the matrix C. M must be at least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of the matrix +*> op( B ) and the number of columns of the matrix C. N must be +*> at least zero. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> On entry, K specifies the number of columns of the matrix +*> op( A ) and the number of rows of the matrix op( B ). K must +*> be at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is +*> k when TRANSA = 'N' or 'n', and is m otherwise. +*> Before entry with TRANSA = 'N' or 'n', the leading m by k +*> part of the array A must contain the matrix A, otherwise +*> the leading k by m part of the array A must contain the +*> matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. When TRANSA = 'N' or 'n' then +*> LDA must be at least max( 1, m ), otherwise LDA must be at +*> least max( 1, k ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is +*> n when TRANSB = 'N' or 'n', and is k otherwise. +*> Before entry with TRANSB = 'N' or 'n', the leading k by n +*> part of the array B must contain the matrix B, otherwise +*> the leading n by k part of the array B must contain the +*> matrix B. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> On entry, LDB specifies the first dimension of B as declared +*> in the calling (sub) program. When TRANSB = 'N' or 'n' then +*> LDB must be at least max( 1, k ), otherwise LDB must be at +*> least max( 1, n ). +*> \endverbatim +*> +*> \param[in] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION. +*> On entry, BETA specifies the scalar beta. When BETA is +*> supplied as zero then C need not be set on input. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is DOUBLE PRECISION array of DIMENSION ( LDC, n ). +*> Before entry, the leading m by n part of the array C must +*> contain the matrix C, except when beta is zero, in which +*> case C need not be set on entry. +*> On exit, the array C is overwritten by the m by n matrix +*> ( alpha*op( A )*op( B ) + beta*C ). +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> On entry, LDC specifies the first dimension of C as declared +*> in the calling (sub) program. LDC must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup double_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 3 Blas routine. +*> +*> -- Written on 8-February-1989. +*> Jack Dongarra, Argonne National Laboratory. +*> Iain Duff, AERE Harwell. +*> Jeremy Du Croz, Numerical Algorithms Group Ltd. +*> Sven Hammarling, Numerical Algorithms Group Ltd. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* +* -- Reference BLAS level3 routine (version 3.6.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + DOUBLE PRECISION ALPHA,BETA + INTEGER K,LDA,LDB,LDC,M,N + CHARACTER TRANSA,TRANSB +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* ===================================================================== +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Local Scalars .. + DOUBLE PRECISION TEMP + INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB + LOGICAL NOTA,NOTB +* .. +* .. Parameters .. + DOUBLE PRECISION ONE,ZERO + PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) +* .. +* +* Set NOTA and NOTB as true if A and B respectively are not +* transposed and set NROWA, NCOLA and NROWB as the number of rows +* and columns of A and the number of rows of B respectively. +* + NOTA = LSAME(TRANSA,'N') + NOTB = LSAME(TRANSB,'N') + IF (NOTA) THEN + NROWA = M + NCOLA = K + ELSE + NROWA = K + NCOLA = M + END IF + IF (NOTB) THEN + NROWB = K + ELSE + NROWB = N + END IF +* +* Test the input parameters. +* + INFO = 0 + IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND. + + (.NOT.LSAME(TRANSA,'T'))) THEN + INFO = 1 + ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .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('DGEMM ',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 if 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 +* +* Form C := alpha*A**T*B + beta*C +* + DO 120 J = 1,N + DO 110 I = 1,M + TEMP = ZERO + DO 100 L = 1,K + TEMP = TEMP + 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 + END IF + ELSE + IF (NOTA) THEN +* +* Form C := alpha*A*B**T + beta*C +* + DO 170 J = 1,N + IF (BETA.EQ.ZERO) THEN + DO 130 I = 1,M + C(I,J) = ZERO + 130 CONTINUE + ELSE IF (BETA.NE.ONE) THEN + DO 140 I = 1,M + C(I,J) = BETA*C(I,J) + 140 CONTINUE + END IF + DO 160 L = 1,K + TEMP = ALPHA*B(J,L) + DO 150 I = 1,M + C(I,J) = C(I,J) + TEMP*A(I,L) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + ELSE +* +* Form C := alpha*A**T*B**T + beta*C +* + DO 200 J = 1,N + DO 190 I = 1,M + TEMP = ZERO + DO 180 L = 1,K + TEMP = TEMP + A(L,I)*B(J,L) + 180 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = ALPHA*TEMP + ELSE + C(I,J) = ALPHA*TEMP + BETA*C(I,J) + END IF + 190 CONTINUE + 200 CONTINUE + END IF + END IF +* + RETURN +* +* End of DGEMM . +* + END diff --git a/internal/testdata/netlib/dgemv.f b/internal/testdata/netlib/dgemv.f new file mode 100644 index 00000000..e04cc07c --- /dev/null +++ b/internal/testdata/netlib/dgemv.f @@ -0,0 +1,330 @@ +*> \brief \b DGEMV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA,BETA +* INTEGER INCX,INCY,LDA,M,N +* CHARACTER TRANS +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*),Y(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGEMV performs one of the matrix-vector operations +*> +*> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, +*> +*> where alpha and beta are scalars, x and y are vectors and A is an +*> m by n matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> On entry, TRANS specifies the operation to be performed as +*> follows: +*> +*> TRANS = 'N' or 'n' y := alpha*A*x + beta*y. +*> +*> TRANS = 'T' or 't' y := alpha*A**T*x + beta*y. +*> +*> TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of the matrix A. +*> M must be at least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry, the leading m by n part of the array A must +*> contain the matrix of coefficients. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, m ). +*> \endverbatim +*> +*> \param[in] X +*> \verbatim +*> X is DOUBLE PRECISION array of DIMENSION at least +*> ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' +*> and at least +*> ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. +*> Before entry, the incremented array X must contain the +*> vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> \endverbatim +*> +*> \param[in] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION. +*> On entry, BETA specifies the scalar beta. When BETA is +*> supplied as zero then Y need not be set on input. +*> \endverbatim +*> +*> \param[in,out] Y +*> \verbatim +*> Y is DOUBLE PRECISION array of DIMENSION at least +*> ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' +*> and at least +*> ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. +*> Before entry with BETA non-zero, the incremented array Y +*> must contain the vector y. On exit, Y is overwritten by the +*> updated vector y. +*> \endverbatim +*> +*> \param[in] INCY +*> \verbatim +*> INCY is INTEGER +*> On entry, INCY specifies the increment for the elements of +*> Y. INCY must not be zero. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup double_blas_level2 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 2 Blas routine. +*> The vector and matrix arguments are not referenced when N = 0, or M = 0 +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) +* +* -- Reference BLAS level2 routine (version 3.6.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + DOUBLE PRECISION ALPHA,BETA + INTEGER INCX,INCY,LDA,M,N + CHARACTER TRANS +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(LDA,*),X(*),Y(*) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE,ZERO + PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) +* .. +* .. Local Scalars .. + DOUBLE PRECISION 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('DGEMV ',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 DGEMV . +* + END diff --git a/internal/testdata/netlib/dlacpy.f b/internal/testdata/netlib/dlacpy.f new file mode 100644 index 00000000..a9a23c94 --- /dev/null +++ b/internal/testdata/netlib/dlacpy.f @@ -0,0 +1,156 @@ +*> \brief \b DLACPY copies all or part of one two-dimensional array to another. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLACPY + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDB, M, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ), B( LDB, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLACPY copies all or part of a two-dimensional matrix A to another +*> matrix B. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies the part of the matrix A to be copied to B. +*> = 'U': Upper triangular part +*> = 'L': Lower triangular part +*> Otherwise: All of the matrix A +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> The m by n matrix A. If UPLO = 'U', only the upper triangle +*> or trapezoid is accessed; if UPLO = 'L', only the lower +*> triangle or trapezoid is accessed. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (LDB,N) +*> On exit, B = A in the locations specified by UPLO. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,M). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) +* +* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* September 2012 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDB, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I, J +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. Executable Statements .. +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO 20 J = 1, N + DO 10 I = 1, MIN( J, M ) + B( I, J ) = A( I, J ) + 10 CONTINUE + 20 CONTINUE + ELSE IF( LSAME( UPLO, 'L' ) ) THEN + DO 40 J = 1, N + DO 30 I = J, M + B( I, J ) = A( I, J ) + 30 CONTINUE + 40 CONTINUE + ELSE + DO 60 J = 1, N + DO 50 I = 1, M + B( I, J ) = A( I, J ) + 50 CONTINUE + 60 CONTINUE + END IF + RETURN +* +* End of DLACPY +* + END diff --git a/internal/testdata/netlib/dlahr2.f b/internal/testdata/netlib/dlahr2.f new file mode 100644 index 00000000..9d15979c --- /dev/null +++ b/internal/testdata/netlib/dlahr2.f @@ -0,0 +1,326 @@ +*> \brief \b DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLAHR2 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) +* +* .. Scalar Arguments .. +* INTEGER K, LDA, LDT, LDY, N, NB +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), +* $ Y( LDY, NB ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1) +*> matrix A so that elements below the k-th subdiagonal are zero. The +*> reduction is performed by an orthogonal similarity transformation +*> Q**T * A * Q. The routine returns the matrices V and T which determine +*> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T. +*> +*> This is an auxiliary routine called by DGEHRD. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The offset for the reduction. Elements below the k-th +*> subdiagonal in the first NB columns are reduced to zero. +*> K < N. +*> \endverbatim +*> +*> \param[in] NB +*> \verbatim +*> NB is INTEGER +*> The number of columns to be reduced. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N-K+1) +*> On entry, the n-by-(n-k+1) general matrix A. +*> On exit, the elements on and above the k-th subdiagonal in +*> the first NB columns are overwritten with the corresponding +*> elements of the reduced matrix; the elements below the k-th +*> subdiagonal, with the array TAU, represent the matrix Q as a +*> product of elementary reflectors. The other columns of A are +*> unchanged. See Further Details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION array, dimension (NB) +*> The scalar factors of the elementary reflectors. See Further +*> Details. +*> \endverbatim +*> +*> \param[out] T +*> \verbatim +*> T is DOUBLE PRECISION array, dimension (LDT,NB) +*> The upper triangular matrix T. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= NB. +*> \endverbatim +*> +*> \param[out] Y +*> \verbatim +*> Y is DOUBLE PRECISION array, dimension (LDY,NB) +*> The n-by-nb matrix Y. +*> \endverbatim +*> +*> \param[in] LDY +*> \verbatim +*> LDY is INTEGER +*> The leading dimension of the array Y. LDY >= N. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The matrix Q is represented as a product of nb elementary reflectors +*> +*> Q = H(1) H(2) . . . H(nb). +*> +*> Each H(i) has the form +*> +*> H(i) = I - tau * v * v**T +*> +*> where tau is a real scalar, and v is a real vector with +*> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in +*> A(i+k+1:n,i), and tau in TAU(i). +*> +*> The elements of the vectors v together form the (n-k+1)-by-nb matrix +*> V which is needed, with T and Y, to apply the transformation to the +*> unreduced part of the matrix, using an update of the form: +*> A := (I - V*T*V**T) * (A - Y*V**T). +*> +*> The contents of A on exit are illustrated by the following example +*> with n = 7, k = 3 and nb = 2: +*> +*> ( a a a a a ) +*> ( a a a a a ) +*> ( a a a a a ) +*> ( h h a a a ) +*> ( v1 h a a a ) +*> ( v1 v2 a a a ) +*> ( v1 v2 a a a ) +*> +*> where a denotes an element of the original matrix A, h denotes a +*> modified element of the upper Hessenberg matrix H, and vi denotes an +*> element of the vector defining H(i). +*> +*> This subroutine is a slight modification of LAPACK-3.0's DLAHRD +*> incorporating improvements proposed by Quintana-Orti and Van de +*> Gejin. Note that the entries of A(1:K,2:NB) differ from those +*> returned by the original LAPACK-3.0's DLAHRD routine. (This +*> subroutine is not backward compatible with LAPACK-3.0's DLAHRD.) +*> \endverbatim +* +*> \par References: +* ================ +*> +*> Gregorio Quintana-Orti and Robert van de Geijn, "Improving the +*> performance of reduction to Hessenberg form," ACM Transactions on +*> Mathematical Software, 32(2):180-194, June 2006. +*> +* ===================================================================== + SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) +* +* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* September 2012 +* +* .. Scalar Arguments .. + INTEGER K, LDA, LDT, LDY, N, NB +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), + $ Y( LDY, NB ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, + $ ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I + DOUBLE PRECISION EI +* .. +* .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DLACPY, + $ DLARFG, DSCAL, DTRMM, DTRMV +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.1 ) + $ RETURN +* + DO 10 I = 1, NB + IF( I.GT.1 ) THEN +* +* Update A(K+1:N,I) +* +* Update I-th column of A - Y * V**T +* + CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, + $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) +* +* Apply I - V * T**T * V**T to this column (call it b) from the +* left, using the last column of T as workspace +* +* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) +* ( V2 ) ( b2 ) +* +* where V1 is unit lower triangular +* +* w := V1**T * b1 +* + CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) + CALL DTRMV( 'Lower', 'Transpose', 'UNIT', + $ I-1, A( K+1, 1 ), + $ LDA, T( 1, NB ), 1 ) +* +* w := w + V2**T * b2 +* + CALL DGEMV( 'Transpose', N-K-I+1, I-1, + $ ONE, A( K+I, 1 ), + $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) +* +* w := T**T * w +* + CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT', + $ I-1, T, LDT, + $ T( 1, NB ), 1 ) +* +* b2 := b2 - V2*w +* + CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, + $ A( K+I, 1 ), + $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) +* +* b1 := b1 - V1*w +* + CALL DTRMV( 'Lower', 'NO TRANSPOSE', + $ 'UNIT', I-1, + $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) + CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) +* + A( K+I-1, I-1 ) = EI + END IF +* +* Generate the elementary reflector H(I) to annihilate +* A(K+I+1:N,I) +* + CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, + $ TAU( I ) ) + EI = A( K+I, I ) + A( K+I, I ) = ONE +* +* Compute Y(K+1:N,I) +* + CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, + $ ONE, A( K+1, I+1 ), + $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) + CALL DGEMV( 'Transpose', N-K-I+1, I-1, + $ ONE, A( K+I, 1 ), LDA, + $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) + CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, + $ Y( K+1, 1 ), LDY, + $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) + CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) +* +* Compute T(1:I,I) +* + CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) + CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT', + $ I-1, T, LDT, + $ T( 1, I ), 1 ) + T( I, I ) = TAU( I ) +* + 10 CONTINUE + A( K+NB, NB ) = EI +* +* Compute Y(1:K,1:NB) +* + CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) + CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', + $ 'UNIT', K, NB, + $ ONE, A( K+1, 1 ), LDA, Y, LDY ) + IF( N.GT.K+NB ) + $ CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, + $ NB, N-K-NB, ONE, + $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, + $ LDY ) + CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', + $ 'NON-UNIT', K, NB, + $ ONE, T, LDT, Y, LDY ) +* + RETURN +* +* End of DLAHR2 +* + END diff --git a/internal/testdata/netlib/dlamch.f b/internal/testdata/netlib/dlamch.f new file mode 100644 index 00000000..22a16218 --- /dev/null +++ b/internal/testdata/netlib/dlamch.f @@ -0,0 +1,189 @@ +*> \brief \b DLAMCH +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAMCH determines double precision machine parameters. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] CMACH +*> \verbatim +*> Specifies the value to be returned by DLAMCH: +*> = 'E' or 'e', DLAMCH := eps +*> = 'S' or 's , DLAMCH := sfmin +*> = 'B' or 'b', DLAMCH := base +*> = 'P' or 'p', DLAMCH := eps*base +*> = 'N' or 'n', DLAMCH := t +*> = 'R' or 'r', DLAMCH := rnd +*> = 'M' or 'm', DLAMCH := emin +*> = 'U' or 'u', DLAMCH := rmin +*> = 'L' or 'l', DLAMCH := emax +*> = 'O' or 'o', DLAMCH := rmax +*> where +*> eps = relative machine precision +*> sfmin = safe minimum, such that 1/sfmin does not overflow +*> base = base of the machine +*> prec = eps*base +*> t = number of (base) digits in the mantissa +*> rnd = 1.0 when rounding occurs in addition, 0.0 otherwise +*> emin = minimum exponent before (gradual) underflow +*> rmin = underflow threshold - base**(emin-1) +*> emax = largest exponent before overflow +*> rmax = overflow threshold - (base**emax)*(1-eps) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) +* +* -- LAPACK auxiliary routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + CHARACTER CMACH +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + DOUBLE PRECISION RND, EPS, SFMIN, SMALL, RMACH +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions .. + INTRINSIC DIGITS, EPSILON, HUGE, MAXEXPONENT, + $ MINEXPONENT, RADIX, TINY +* .. +* .. Executable Statements .. +* +* +* Assume rounding, not chopping. Always. +* + RND = ONE +* + IF( ONE.EQ.RND ) THEN + EPS = EPSILON(ZERO) * 0.5 + ELSE + EPS = EPSILON(ZERO) + END IF +* + IF( LSAME( CMACH, 'E' ) ) THEN + RMACH = EPS + ELSE IF( LSAME( CMACH, 'S' ) ) THEN + SFMIN = TINY(ZERO) + SMALL = ONE / HUGE(ZERO) + IF( SMALL.GE.SFMIN ) THEN +* +* Use SMALL plus a bit, to avoid the possibility of rounding +* causing overflow when computing 1/sfmin. +* + SFMIN = SMALL*( ONE+EPS ) + END IF + RMACH = SFMIN + ELSE IF( LSAME( CMACH, 'B' ) ) THEN + RMACH = RADIX(ZERO) + ELSE IF( LSAME( CMACH, 'P' ) ) THEN + RMACH = EPS * RADIX(ZERO) + ELSE IF( LSAME( CMACH, 'N' ) ) THEN + RMACH = DIGITS(ZERO) + ELSE IF( LSAME( CMACH, 'R' ) ) THEN + RMACH = RND + ELSE IF( LSAME( CMACH, 'M' ) ) THEN + RMACH = MINEXPONENT(ZERO) + ELSE IF( LSAME( CMACH, 'U' ) ) THEN + RMACH = tiny(zero) + ELSE IF( LSAME( CMACH, 'L' ) ) THEN + RMACH = MAXEXPONENT(ZERO) + ELSE IF( LSAME( CMACH, 'O' ) ) THEN + RMACH = HUGE(ZERO) + ELSE + RMACH = ZERO + END IF +* + DLAMCH = RMACH + RETURN +* +* End of DLAMCH +* + END +************************************************************************ +*> \brief \b DLAMC3 +*> \details +*> \b Purpose: +*> \verbatim +*> DLAMC3 is intended to force A and B to be stored prior to doing +*> the addition of A and B , for use in situations where optimizers +*> might hold one of these in a register. +*> \endverbatim +*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. +*> \date November 2015 +*> \ingroup auxOTHERauxiliary +*> +*> \param[in] A +*> \verbatim +*> A is a DOUBLE PRECISION +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is a DOUBLE PRECISION +*> The values A and B. +*> \endverbatim +*> + DOUBLE PRECISION FUNCTION DLAMC3( A, B ) +* +* -- LAPACK auxiliary routine (version 3.6.0) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2010 +* +* .. Scalar Arguments .. + DOUBLE PRECISION A, B +* .. +* ===================================================================== +* +* .. Executable Statements .. +* + DLAMC3 = A + B +* + RETURN +* +* End of DLAMC3 +* + END +* +************************************************************************ diff --git a/internal/testdata/netlib/dlapy2.f b/internal/testdata/netlib/dlapy2.f new file mode 100644 index 00000000..d43b0d5d --- /dev/null +++ b/internal/testdata/netlib/dlapy2.f @@ -0,0 +1,104 @@ +*> \brief \b DLAPY2 returns sqrt(x2+y2). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLAPY2 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION X, Y +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary +*> overflow. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] X +*> \verbatim +*> X is DOUBLE PRECISION +*> \endverbatim +*> +*> \param[in] Y +*> \verbatim +*> Y is DOUBLE PRECISION +*> X and Y specify the values x and y. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) +* +* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* September 2012 +* +* .. Scalar Arguments .. + DOUBLE PRECISION X, Y +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D0 ) +* .. +* .. Local Scalars .. + DOUBLE PRECISION W, XABS, YABS, Z +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +* .. +* .. Executable Statements .. +* + XABS = ABS( X ) + YABS = ABS( Y ) + W = MAX( XABS, YABS ) + Z = MIN( XABS, YABS ) + IF( Z.EQ.ZERO ) THEN + DLAPY2 = W + ELSE + DLAPY2 = W*SQRT( ONE+( Z / W )**2 ) + END IF + RETURN +* +* End of DLAPY2 +* + END diff --git a/internal/testdata/netlib/dlarfg.f b/internal/testdata/netlib/dlarfg.f new file mode 100644 index 00000000..ce91d33c --- /dev/null +++ b/internal/testdata/netlib/dlarfg.f @@ -0,0 +1,196 @@ +*> \brief \b DLARFG generates an elementary reflector (Householder matrix). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARFG + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* DOUBLE PRECISION ALPHA, TAU +* .. +* .. Array Arguments .. +* DOUBLE PRECISION X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLARFG generates a real elementary reflector H of order n, such +*> that +*> +*> H * ( alpha ) = ( beta ), H**T * H = I. +*> ( x ) ( 0 ) +*> +*> where alpha and beta are scalars, and x is an (n-1)-element real +*> vector. H is represented in the form +*> +*> H = I - tau * ( 1 ) * ( 1 v**T ) , +*> ( v ) +*> +*> where tau is a real scalar and v is a real (n-1)-element +*> vector. +*> +*> If the elements of x are all zero, then tau = 0 and H is taken to be +*> the unit matrix. +*> +*> Otherwise 1 <= tau <= 2. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the elementary reflector. +*> \endverbatim +*> +*> \param[in,out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION +*> On entry, the value alpha. +*> On exit, it is overwritten with the value beta. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension +*> (1+(N-2)*abs(INCX)) +*> On entry, the vector x. +*> On exit, it is overwritten with the vector v. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between elements of X. INCX > 0. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION +*> The value tau. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +* ===================================================================== + SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.4.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* September 2012 +* +* .. Scalar Arguments .. + INTEGER INCX, N + DOUBLE PRECISION ALPHA, TAU +* .. +* .. Array Arguments .. + DOUBLE PRECISION X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 + EXTERNAL DLAMCH, DLAPY2, DNRM2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, SIGN +* .. +* .. External Subroutines .. + EXTERNAL DSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.1 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = DNRM2( N-1, X, INCX ) +* + IF( XNORM.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + RSAFMN = ONE / SAFMIN + 10 CONTINUE + KNT = KNT + 1 + CALL DSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHA = ALPHA*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = DNRM2( N-1, X, INCX ) + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + END IF + TAU = ( BETA-ALPHA ) / BETA + CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of DLARFG +* + END diff --git a/internal/testdata/netlib/dnrm2.f b/internal/testdata/netlib/dnrm2.f new file mode 100644 index 00000000..5ea257a2 --- /dev/null +++ b/internal/testdata/netlib/dnrm2.f @@ -0,0 +1,112 @@ +*> \brief \b DNRM2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DNRM2(N,X,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION X(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DNRM2 returns the euclidean norm of a vector via the function +*> name, so that +*> +*> DNRM2 := sqrt( x'*x ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> -- This version written on 25-October-1982. +*> Modified on 14-October-1993 to inline the call to DLASSQ. +*> Sven Hammarling, Nag Ltd. +*> \endverbatim +*> +* ===================================================================== + DOUBLE PRECISION FUNCTION DNRM2(N,X,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 +* +* .. Scalar Arguments .. + INTEGER INCX,N +* .. +* .. Array Arguments .. + DOUBLE PRECISION X(*) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE,ZERO + PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) +* .. +* .. Local Scalars .. + DOUBLE PRECISION ABSXI,NORM,SCALE,SSQ + INTEGER IX +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS,SQRT +* .. + IF (N.LT.1 .OR. INCX.LT.1) THEN + NORM = ZERO + ELSE IF (N.EQ.1) THEN + NORM = ABS(X(1)) + ELSE + SCALE = ZERO + SSQ = ONE +* The following loop is equivalent to this call to the LAPACK +* auxiliary routine: +* CALL DLASSQ( N, X, INCX, SCALE, SSQ ) +* + DO 10 IX = 1,1 + (N-1)*INCX,INCX + IF (X(IX).NE.ZERO) THEN + ABSXI = ABS(X(IX)) + IF (SCALE.LT.ABSXI) THEN + SSQ = ONE + SSQ* (SCALE/ABSXI)**2 + SCALE = ABSXI + ELSE + SSQ = SSQ + (ABSXI/SCALE)**2 + END IF + END IF + 10 CONTINUE + NORM = SCALE*SQRT(SSQ) + END IF +* + DNRM2 = NORM + RETURN +* +* End of DNRM2. +* + END diff --git a/internal/testdata/netlib/dscal.f b/internal/testdata/netlib/dscal.f new file mode 100644 index 00000000..3337de8e --- /dev/null +++ b/internal/testdata/netlib/dscal.f @@ -0,0 +1,110 @@ +*> \brief \b DSCAL +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DSCAL(N,DA,DX,INCX) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION DA +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSCAL scales a vector by a constant. +*> uses unrolled loops for increment equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 3/93 to return if incx .le. 0. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DSCAL(N,DA,DX,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 +* +* .. Scalar Arguments .. + DOUBLE PRECISION DA + INTEGER INCX,N +* .. +* .. Array Arguments .. + DOUBLE PRECISION DX(*) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I,M,MP1,NINCX +* .. +* .. Intrinsic Functions .. + INTRINSIC MOD +* .. + 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,5) + IF (M.NE.0) THEN + DO I = 1,M + DX(I) = DA*DX(I) + END DO + IF (N.LT.5) RETURN + END IF + MP1 = M + 1 + DO I = MP1,N,5 + DX(I) = DA*DX(I) + DX(I+1) = DA*DX(I+1) + DX(I+2) = DA*DX(I+2) + DX(I+3) = DA*DX(I+3) + DX(I+4) = DA*DX(I+4) + END DO + ELSE +* +* code for increment not equal to 1 +* + NINCX = N*INCX + DO I = 1,NINCX,INCX + DX(I) = DA*DX(I) + END DO + END IF + RETURN + END diff --git a/internal/testdata/netlib/dtrmm.f b/internal/testdata/netlib/dtrmm.f new file mode 100644 index 00000000..cbd5ce70 --- /dev/null +++ b/internal/testdata/netlib/dtrmm.f @@ -0,0 +1,415 @@ +*> \brief \b DTRMM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA +* INTEGER LDA,LDB,M,N +* CHARACTER DIAG,SIDE,TRANSA,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),B(LDB,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRMM performs one of the matrix-matrix operations +*> +*> B := alpha*op( A )*B, or B := alpha*B*op( A ), +*> +*> where alpha is a scalar, B is an m by n matrix, A is a unit, or +*> non-unit, upper or lower triangular matrix and op( A ) is one of +*> +*> op( A ) = A or op( A ) = A**T. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> On entry, SIDE specifies whether op( A ) multiplies B from +*> the left or right as follows: +*> +*> SIDE = 'L' or 'l' B := alpha*op( A )*B. +*> +*> SIDE = 'R' or 'r' B := alpha*B*op( A ). +*> \endverbatim +*> +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix A is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is CHARACTER*1 +*> On entry, TRANSA specifies the form of op( A ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSA = 'N' or 'n' op( A ) = A. +*> +*> TRANSA = 'T' or 't' op( A ) = A**T. +*> +*> TRANSA = 'C' or 'c' op( A ) = A**T. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit triangular +*> as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of B. M must be at +*> least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of B. N must be +*> at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. When alpha is +*> zero then A is not referenced and B need not be set before +*> entry. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m +*> when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. +*> Before entry with UPLO = 'U' or 'u', the leading k by k +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading k by k +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. When SIDE = 'L' or 'l' then +*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +*> then LDA must be at least max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array of DIMENSION ( LDB, n ). +*> Before entry, the leading m by n part of the array B must +*> contain the matrix B, and on exit is overwritten by the +*> transformed matrix. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> On entry, LDB specifies the first dimension of B as declared +*> in the calling (sub) program. LDB must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 3 Blas routine. +*> +*> -- Written on 8-February-1989. +*> Jack Dongarra, Argonne National Laboratory. +*> Iain Duff, AERE Harwell. +*> Jeremy Du Croz, Numerical Algorithms Group Ltd. +*> Sven Hammarling, Numerical Algorithms Group Ltd. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* -- 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 +* +* .. Scalar Arguments .. + DOUBLE PRECISION ALPHA + INTEGER LDA,LDB,M,N + CHARACTER DIAG,SIDE,TRANSA,UPLO +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(LDA,*),B(LDB,*) +* .. +* +* ===================================================================== +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Local Scalars .. + DOUBLE PRECISION TEMP + INTEGER I,INFO,J,K,NROWA + LOGICAL LSIDE,NOUNIT,UPPER +* .. +* .. Parameters .. + DOUBLE PRECISION ONE,ZERO + PARAMETER (ONE=1.0D+0,ZERO=0.0D+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('DTRMM ',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*A*B. +* + IF (UPPER) THEN + DO 50 J = 1,N + DO 40 K = 1,M + IF (B(K,J).NE.ZERO) THEN + TEMP = ALPHA*B(K,J) + DO 30 I = 1,K - 1 + B(I,J) = B(I,J) + TEMP*A(I,K) + 30 CONTINUE + IF (NOUNIT) TEMP = TEMP*A(K,K) + B(K,J) = TEMP + END IF + 40 CONTINUE + 50 CONTINUE + ELSE + DO 80 J = 1,N + DO 70 K = M,1,-1 + IF (B(K,J).NE.ZERO) THEN + TEMP = ALPHA*B(K,J) + B(K,J) = TEMP + IF (NOUNIT) B(K,J) = B(K,J)*A(K,K) + DO 60 I = K + 1,M + B(I,J) = B(I,J) + TEMP*A(I,K) + 60 CONTINUE + END IF + 70 CONTINUE + 80 CONTINUE + END IF + ELSE +* +* Form B := alpha*A**T*B. +* + IF (UPPER) THEN + DO 110 J = 1,N + DO 100 I = M,1,-1 + TEMP = B(I,J) + IF (NOUNIT) TEMP = TEMP*A(I,I) + DO 90 K = 1,I - 1 + TEMP = TEMP + A(K,I)*B(K,J) + 90 CONTINUE + B(I,J) = ALPHA*TEMP + 100 CONTINUE + 110 CONTINUE + ELSE + DO 140 J = 1,N + DO 130 I = 1,M + TEMP = B(I,J) + IF (NOUNIT) TEMP = TEMP*A(I,I) + DO 120 K = I + 1,M + TEMP = TEMP + A(K,I)*B(K,J) + 120 CONTINUE + B(I,J) = ALPHA*TEMP + 130 CONTINUE + 140 CONTINUE + END IF + END IF + ELSE + IF (LSAME(TRANSA,'N')) THEN +* +* Form B := alpha*B*A. +* + IF (UPPER) THEN + DO 180 J = N,1,-1 + TEMP = ALPHA + IF (NOUNIT) TEMP = TEMP*A(J,J) + DO 150 I = 1,M + B(I,J) = TEMP*B(I,J) + 150 CONTINUE + DO 170 K = 1,J - 1 + IF (A(K,J).NE.ZERO) THEN + TEMP = ALPHA*A(K,J) + DO 160 I = 1,M + B(I,J) = B(I,J) + TEMP*B(I,K) + 160 CONTINUE + END IF + 170 CONTINUE + 180 CONTINUE + ELSE + DO 220 J = 1,N + TEMP = ALPHA + IF (NOUNIT) TEMP = TEMP*A(J,J) + DO 190 I = 1,M + B(I,J) = TEMP*B(I,J) + 190 CONTINUE + DO 210 K = J + 1,N + IF (A(K,J).NE.ZERO) THEN + TEMP = ALPHA*A(K,J) + DO 200 I = 1,M + B(I,J) = B(I,J) + TEMP*B(I,K) + 200 CONTINUE + END IF + 210 CONTINUE + 220 CONTINUE + END IF + ELSE +* +* Form B := alpha*B*A**T. +* + IF (UPPER) THEN + DO 260 K = 1,N + DO 240 J = 1,K - 1 + IF (A(J,K).NE.ZERO) THEN + TEMP = ALPHA*A(J,K) + DO 230 I = 1,M + B(I,J) = B(I,J) + TEMP*B(I,K) + 230 CONTINUE + END IF + 240 CONTINUE + TEMP = ALPHA + IF (NOUNIT) TEMP = TEMP*A(K,K) + IF (TEMP.NE.ONE) THEN + DO 250 I = 1,M + B(I,K) = TEMP*B(I,K) + 250 CONTINUE + END IF + 260 CONTINUE + ELSE + DO 300 K = N,1,-1 + DO 280 J = K + 1,N + IF (A(J,K).NE.ZERO) THEN + TEMP = ALPHA*A(J,K) + DO 270 I = 1,M + B(I,J) = B(I,J) + TEMP*B(I,K) + 270 CONTINUE + END IF + 280 CONTINUE + TEMP = ALPHA + IF (NOUNIT) TEMP = TEMP*A(K,K) + IF (TEMP.NE.ONE) THEN + DO 290 I = 1,M + B(I,K) = TEMP*B(I,K) + 290 CONTINUE + END IF + 300 CONTINUE + END IF + END IF + END IF +* + RETURN +* +* End of DTRMM . +* + END diff --git a/internal/testdata/netlib/dtrmv.f b/internal/testdata/netlib/dtrmv.f new file mode 100644 index 00000000..71459fe7 --- /dev/null +++ b/internal/testdata/netlib/dtrmv.f @@ -0,0 +1,342 @@ +*> \brief \b DTRMV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,LDA,N +* CHARACTER DIAG,TRANS,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRMV performs one of the matrix-vector operations +*> +*> x := A*x, or x := A**T*x, +*> +*> where x is an n element vector and A is an n by n unit, or non-unit, +*> upper or lower triangular matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> On entry, TRANS specifies the operation to be performed as +*> follows: +*> +*> TRANS = 'N' or 'n' x := A*x. +*> +*> TRANS = 'T' or 't' x := A**T*x. +*> +*> TRANS = 'C' or 'c' x := A**T*x. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit +*> triangular as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the order of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry with UPLO = 'U' or 'u', the leading n by n +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading n by n +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array of dimension at least +*> ( 1 + ( n - 1 )*abs( INCX ) ). +*> Before entry, the incremented array X must contain the n +*> element vector x. On exit, X is overwritten with the +*> tranformed vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level2 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 2 Blas routine. +*> The vector and matrix arguments are not referenced when N = 0, or M = 0 +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* -- 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 +* +* .. Scalar Arguments .. + INTEGER INCX,LDA,N + CHARACTER DIAG,TRANS,UPLO +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(LDA,*),X(*) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER (ZERO=0.0D+0) +* .. +* .. Local Scalars .. + DOUBLE PRECISION 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('DTRMV ',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 DTRMV . +* + END diff --git a/internal/testdata/netlib/lsame.f b/internal/testdata/netlib/lsame.f new file mode 100644 index 00000000..315304c3 --- /dev/null +++ b/internal/testdata/netlib/lsame.f @@ -0,0 +1,125 @@ +*> \brief \b LSAME +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* LOGICAL FUNCTION LSAME( CA, CB ) +* +* .. Scalar Arguments .. +* CHARACTER CA, CB +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> LSAME returns .TRUE. if CA is the same letter as CB regardless of +*> case. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] CA +*> \verbatim +*> \endverbatim +*> +*> \param[in] CB +*> \verbatim +*> CA and CB specify the single characters to be compared. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + LOGICAL FUNCTION LSAME( CA, CB ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER CA, CB +* .. +* +* ===================================================================== +* +* .. Intrinsic Functions .. + INTRINSIC ICHAR +* .. +* .. Local Scalars .. + INTEGER INTA, INTB, ZCODE +* .. +* .. Executable Statements .. +* +* 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 diff --git a/internal/testdata/netlib/netlib.go b/internal/testdata/netlib/netlib.go new file mode 100644 index 00000000..01173873 --- /dev/null +++ b/internal/testdata/netlib/netlib.go @@ -0,0 +1,24 @@ +// Copyright ©2016 The gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package netlib + +// void dlahr2_(int* n, int* k, int* nb, double* a, int* lda, double* tau, double* t, int* ldt, double* y, int* ldy); +import "C" + +func Dlahr2(n, k, nb int, a []float64, lda int, tau, t []float64, ldt int, y []float64, ldy int) { + func() { + n := C.int(n) + k := C.int(k) + nb := C.int(nb) + lda := C.int(lda) + ldt := C.int(ldt) + ldy := C.int(ldy) + C.dlahr2_((*C.int)(&n), (*C.int)(&k), (*C.int)(&nb), + (*C.double)(&a[0]), (*C.int)(&lda), + (*C.double)(&tau[0]), + (*C.double)(&t[0]), (*C.int)(&ldt), + (*C.double)(&y[0]), (*C.int)(&ldy)) + }() +} diff --git a/internal/testdata/netlib/xerbla.f b/internal/testdata/netlib/xerbla.f new file mode 100644 index 00000000..eb1c037d --- /dev/null +++ b/internal/testdata/netlib/xerbla.f @@ -0,0 +1,89 @@ +*> \brief \b XERBLA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE XERBLA( SRNAME, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER*(*) SRNAME +* INTEGER INFO +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> XERBLA is an error handler for the LAPACK routines. +*> It is called by an LAPACK routine if an input parameter has an +*> invalid value. A message is printed and execution stops. +*> +*> Installers may consider modifying the STOP statement in order to +*> call system-specific exception-handling facilities. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SRNAME +*> \verbatim +*> SRNAME is CHARACTER*(*) +*> The name of the routine which called XERBLA. +*> \endverbatim +*> +*> \param[in] INFO +*> \verbatim +*> INFO is INTEGER +*> The position of the invalid parameter in the parameter list +*> of the calling routine. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup aux_blas +* +* ===================================================================== + SUBROUTINE XERBLA( SRNAME, INFO ) +* +* -- 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 +* +* .. 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