diff --git a/internal/testdata/dlaqr5test/dgemm.f b/internal/testdata/dlaqr5test/dgemm.f
deleted file mode 100644
index 4bae243a..00000000
--- a/internal/testdata/dlaqr5test/dgemm.f
+++ /dev/null
@@ -1,384 +0,0 @@
-*> \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/dlaqr5test/dlacpy.f b/internal/testdata/dlaqr5test/dlacpy.f
deleted file mode 100644
index a9a23c94..00000000
--- a/internal/testdata/dlaqr5test/dlacpy.f
+++ /dev/null
@@ -1,156 +0,0 @@
-*> \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/dlaqr5test/dlamch.f b/internal/testdata/dlaqr5test/dlamch.f
deleted file mode 100644
index 22a16218..00000000
--- a/internal/testdata/dlaqr5test/dlamch.f
+++ /dev/null
@@ -1,189 +0,0 @@
-*> \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/dlaqr5test/dlapy2.f b/internal/testdata/dlaqr5test/dlapy2.f
deleted file mode 100644
index d43b0d5d..00000000
--- a/internal/testdata/dlaqr5test/dlapy2.f
+++ /dev/null
@@ -1,104 +0,0 @@
-*> \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/dlaqr5test/dlaqr5data.zip b/internal/testdata/dlaqr5test/dlaqr5data.zip
deleted file mode 100644
index 215bfb61..00000000
Binary files a/internal/testdata/dlaqr5test/dlaqr5data.zip and /dev/null differ
diff --git a/internal/testdata/dlaqr5test/dlarfg.f b/internal/testdata/dlaqr5test/dlarfg.f
deleted file mode 100644
index ce91d33c..00000000
--- a/internal/testdata/dlaqr5test/dlarfg.f
+++ /dev/null
@@ -1,196 +0,0 @@
-*> \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/dlaqr5test/dnrm2.f b/internal/testdata/dlaqr5test/dnrm2.f
deleted file mode 100644
index 5ea257a2..00000000
--- a/internal/testdata/dlaqr5test/dnrm2.f
+++ /dev/null
@@ -1,112 +0,0 @@
-*> \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/dlaqr5test/dscal.f b/internal/testdata/dlaqr5test/dscal.f
deleted file mode 100644
index 3337de8e..00000000
--- a/internal/testdata/dlaqr5test/dscal.f
+++ /dev/null
@@ -1,110 +0,0 @@
-*> \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/dlaqr5test/dtrmm.f b/internal/testdata/dlaqr5test/dtrmm.f
deleted file mode 100644
index cbd5ce70..00000000
--- a/internal/testdata/dlaqr5test/dtrmm.f
+++ /dev/null
@@ -1,415 +0,0 @@
-*> \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/dlaqr5test/lsame.f b/internal/testdata/dlaqr5test/lsame.f
deleted file mode 100644
index 315304c3..00000000
--- a/internal/testdata/dlaqr5test/lsame.f
+++ /dev/null
@@ -1,125 +0,0 @@
-*> \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/dlaqr5test/main.go b/internal/testdata/dlaqr5test/main.go
index a34dbbe2..1ef02faf 100644
--- a/internal/testdata/dlaqr5test/main.go
+++ b/internal/testdata/dlaqr5test/main.go
@@ -2,81 +2,49 @@
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
-// This program generates test data for Dlaqr5. Test cases are stored as text
-// files inside dlaqr5data.zip archive which is then read by testlapack/dlaqr5.go.
+// This program generates test data for Dlaqr5. Test cases are stored in
+// gzip-compressed JSON file testlapack/testdata/dlaqr5data.json.gz which is
+// read during testing by testlapack/dlaqr5.go.
//
// This program uses cgo to call Fortran version of DLAQR5. Therefore, matrices
-// generated by hessrand and eye are in column-major format but are written into
-// test case files in row-major format. See writeCase and writeCaseWant for
-// details.
+// passed to the Fortran routine are in column-major format but are written into
+// the output file in row-major format.
package main
-// void dlaqr5_(int* wantt, int* wantz, int* kacc22, int* n, int* ktop, int* kbot, int* nshfts,
-// double* sr, double* si, double* h, int* ldh, int* iloz, int* ihiz,
-// double* z, int* ldz, double* v, int* ldv, double* u, int* ldu,
-// int* nv, double* wv, int* ldwv, int* nh, double* wh, int* ldwh);
-import "C"
-
import (
- "archive/zip"
- "fmt"
- "io"
+ "compress/gzip"
+ "encoding/json"
"log"
"math/rand"
"os"
+ "path/filepath"
+
+ "github.com/gonum/lapack/internal/testdata/netlib"
)
-func fortranDlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot int, nshfts int, sr, si []float64, h []float64,
- ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int,
- u []float64, ldu int, nh int, wh []float64, ldwh int, nv int, wv []float64, ldwv int) {
- func() {
- wt := C.int(0)
- if wantt {
- wt = 1
- }
- wz := C.int(0)
- if wantz {
- wz = 1
- }
- kacc22 := C.int(kacc22)
- n := C.int(n)
- ktop := C.int(ktop + 1)
- kbot := C.int(kbot + 1)
- nshfts := C.int(nshfts)
- ldh := C.int(ldh)
- iloz := C.int(iloz + 1)
- ihiz := C.int(ihiz + 1)
- ldz := C.int(ldz)
- ldv := C.int(ldv)
- ldu := C.int(ldu)
- nh := C.int(nh)
- ldwh := C.int(ldwh)
- nv := C.int(nv)
- ldwv := C.int(ldwv)
- C.dlaqr5_((*C.int)(&wt), (*C.int)(&wz), (*C.int)(&kacc22),
- (*C.int)(&n), (*C.int)(&ktop), (*C.int)(&kbot),
- (*C.int)(&nshfts), (*C.double)(&sr[0]), (*C.double)(&si[0]),
- (*C.double)(&h[0]), (*C.int)(&ldh),
- (*C.int)(&iloz), (*C.int)(&ihiz), (*C.double)(&z[0]), (*C.int)(&ldz),
- (*C.double)(&v[0]), (*C.int)(&ldv),
- (*C.double)(&u[0]), (*C.int)(&ldu),
- (*C.int)(&nh), (*C.double)(&wh[0]), (*C.int)(&ldwh),
- (*C.int)(&nv), (*C.double)(&wv[0]), (*C.int)(&ldwv))
- }()
+type Dlaqr5Test struct {
+ WantT bool
+ N int
+ NShifts int
+ KTop, KBot int
+ ShiftR, ShiftI []float64
+ H []float64
+
+ HWant []float64
+ ZWant []float64
}
func main() {
- const tmpl = "wantt%v_n%v_nshfts%v_ktop%v_kbot%v.txt"
-
- file, err := os.Create("dlaqr5data.zip")
+ file, err := os.Create(filepath.FromSlash("../../../testlapack/testdata/dlaqr5data.json.gz"))
if err != nil {
log.Fatal(err)
}
defer file.Close()
- zipfile := zip.NewWriter(file)
+ w := gzip.NewWriter(file)
rnd := rand.New(rand.NewSource(1))
+ var tests []Dlaqr5Test
for _, wantt := range []bool{true, false} {
for _, n := range []int{2, 3, 4, 5, 6, 7, 11} {
for k := 0; k <= min(5, n); k++ {
@@ -89,14 +57,13 @@ func main() {
sr, si := shiftpairs(npairs, rnd)
nshfts := len(sr)
- v := make([]float64, nshfts/2*3)
- u := make([]float64, (3*nshfts-3)*(3*nshfts-3))
+ v := genrand(nshfts/2, 3, rnd)
+ u := genrand(3*nshfts-3, 3*nshfts-3, rnd)
+ wh := genrand(3*nshfts-3, n, rnd)
nh := n
- wh := make([]float64, (3*nshfts-3)*n)
+ wv := genrand(n, 3*nshfts-3, rnd)
nv := n
- wv := make([]float64, n*(3*nshfts-3))
- z := eye(n)
h := hessrand(n, rnd)
if ktop > 0 {
h[ktop+(ktop-1)*n] = 0
@@ -104,34 +71,54 @@ func main() {
if kbot < n-1 {
h[kbot+1+kbot*n] = 0
}
+ hin := make([]float64, len(h))
+ copy(hin, h)
+ z := eye(n)
- w, err := zipfile.Create(fmt.Sprintf(tmpl, wantt, n, nshfts, ktop, kbot))
- if err != nil {
- log.Fatal(err)
- }
- writeCase(w, wantt, n, nshfts, ktop, kbot, sr, si, h)
- fortranDlaqr5(wantt, true, 2,
- n, ktop, kbot,
+ netlib.Dlaqr5(wantt, true, 2,
+ n, ktop+1, kbot+1,
nshfts, sr, si,
h, n,
- 0, n-1, z, n,
+ 1, n, z, n,
v, 3,
u, 3*nshfts-3,
nh, wh, nh,
nv, wv, 3*nshfts-3)
- writeCaseWant(w, n, h, z)
+
+ tests = append(tests, Dlaqr5Test{
+ WantT: wantt,
+ N: n,
+ NShifts: nshfts,
+ KTop: ktop,
+ KBot: kbot,
+ ShiftR: sr,
+ ShiftI: si,
+ H: rowMajor(n, n, hin),
+ HWant: rowMajor(n, n, h),
+ ZWant: rowMajor(n, n, z),
+ })
}
}
}
}
}
+ json.NewEncoder(w).Encode(tests)
- err = zipfile.Close()
+ err = w.Close()
if err != nil {
log.Fatal(err)
}
}
+// genrand returns a general r×c matrix with random entries.
+func genrand(r, c int, rnd *rand.Rand) []float64 {
+ m := make([]float64, r*c)
+ for i := range m {
+ m[i] = rnd.NormFloat64()
+ }
+ return m
+}
+
// eye returns an identity matrix of order n.
func eye(n int) []float64 {
m := make([]float64, n*n)
@@ -176,36 +163,18 @@ func shiftpairs(k int, rnd *rand.Rand) (sr, si []float64) {
return sr, si
}
-// writeCase writes into w given data with one value per line. h is assumed in
-// column-major order and is written in row-major.
-func writeCase(w io.Writer, wantt bool, n, nshfts, ktop, kbot int, sr, si, h []float64) {
- fmt.Fprintln(w, wantt, n, nshfts, ktop, kbot)
- for _, v := range sr {
- fmt.Fprintln(w, v)
+// rowMajor returns the given r×c column-major matrix a in row-major format.
+func rowMajor(r, c int, a []float64) []float64 {
+ if len(a) != r*c {
+ panic("testdata: slice length mismatch")
}
- for _, v := range si {
- fmt.Fprintln(w, v)
- }
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- fmt.Fprintln(w, h[i+j*n])
- }
- }
-}
-
-// writeCaseWant writes into w given data with one value per line. h and z are
-// assumed in column-major order and are written in row-major.
-func writeCaseWant(w io.Writer, n int, h, z []float64) {
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- fmt.Fprintln(w, h[i+j*n])
- }
- }
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- fmt.Fprintln(w, z[i+j*n])
+ m := make([]float64, len(a))
+ for i := 0; i < r; i++ {
+ for j := 0; j < c; j++ {
+ m[i*c+j] = a[i+j*r]
}
}
+ return m
}
func min(a, b int) int {
diff --git a/internal/testdata/dlaqr5test/run.bash b/internal/testdata/dlaqr5test/run.bash
deleted file mode 100755
index 6cf6ed74..00000000
--- a/internal/testdata/dlaqr5test/run.bash
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/usr/bin/env bash
-
-go build
-./dlaqr5test
-rm -f dlaqr5test
diff --git a/internal/testdata/dlaqr5test/xerbla.f b/internal/testdata/dlaqr5test/xerbla.f
deleted file mode 100644
index eb1c037d..00000000
--- a/internal/testdata/dlaqr5test/xerbla.f
+++ /dev/null
@@ -1,89 +0,0 @@
-*> \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
diff --git a/internal/testdata/dlaqr5test/dlabad.f b/internal/testdata/netlib/dlabad.f
similarity index 100%
rename from internal/testdata/dlaqr5test/dlabad.f
rename to internal/testdata/netlib/dlabad.f
diff --git a/internal/testdata/dlaqr5test/dlaqr1.f b/internal/testdata/netlib/dlaqr1.f
similarity index 100%
rename from internal/testdata/dlaqr5test/dlaqr1.f
rename to internal/testdata/netlib/dlaqr1.f
diff --git a/internal/testdata/dlaqr5test/dlaqr5.f b/internal/testdata/netlib/dlaqr5.f
similarity index 100%
rename from internal/testdata/dlaqr5test/dlaqr5.f
rename to internal/testdata/netlib/dlaqr5.f
diff --git a/internal/testdata/dlaqr5test/dlaset.f b/internal/testdata/netlib/dlaset.f
similarity index 100%
rename from internal/testdata/dlaqr5test/dlaset.f
rename to internal/testdata/netlib/dlaset.f
diff --git a/internal/testdata/netlib/netlib.go b/internal/testdata/netlib/netlib.go
index 01173873..182e074a 100644
--- a/internal/testdata/netlib/netlib.go
+++ b/internal/testdata/netlib/netlib.go
@@ -5,6 +5,11 @@
package netlib
// void dlahr2_(int* n, int* k, int* nb, double* a, int* lda, double* tau, double* t, int* ldt, double* y, int* ldy);
+//
+// void dlaqr5_(int* wantt, int* wantz, int* kacc22, int* n, int* ktop, int* kbot, int* nshfts,
+// double* sr, double* si, double* h, int* ldh, int* iloz, int* ihiz,
+// double* z, int* ldz, double* v, int* ldv, double* u, int* ldu,
+// int* nv, double* wv, int* ldwv, int* nh, double* wh, int* ldwh);
import "C"
func Dlahr2(n, k, nb int, a []float64, lda int, tau, t []float64, ldt int, y []float64, ldy int) {
@@ -22,3 +27,42 @@ func Dlahr2(n, k, nb int, a []float64, lda int, tau, t []float64, ldt int, y []f
(*C.double)(&y[0]), (*C.int)(&ldy))
}()
}
+
+func Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot int, nshfts int, sr, si []float64, h []float64,
+ ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int,
+ u []float64, ldu int, nh int, wh []float64, ldwh int, nv int, wv []float64, ldwv int) {
+ func() {
+ wt := C.int(0)
+ if wantt {
+ wt = 1
+ }
+ wz := C.int(0)
+ if wantz {
+ wz = 1
+ }
+ kacc22 := C.int(kacc22)
+ n := C.int(n)
+ ktop := C.int(ktop)
+ kbot := C.int(kbot)
+ nshfts := C.int(nshfts)
+ ldh := C.int(ldh)
+ iloz := C.int(iloz)
+ ihiz := C.int(ihiz)
+ ldz := C.int(ldz)
+ ldv := C.int(ldv)
+ ldu := C.int(ldu)
+ nh := C.int(nh)
+ ldwh := C.int(ldwh)
+ nv := C.int(nv)
+ ldwv := C.int(ldwv)
+ C.dlaqr5_((*C.int)(&wt), (*C.int)(&wz), (*C.int)(&kacc22),
+ (*C.int)(&n), (*C.int)(&ktop), (*C.int)(&kbot),
+ (*C.int)(&nshfts), (*C.double)(&sr[0]), (*C.double)(&si[0]),
+ (*C.double)(&h[0]), (*C.int)(&ldh),
+ (*C.int)(&iloz), (*C.int)(&ihiz), (*C.double)(&z[0]), (*C.int)(&ldz),
+ (*C.double)(&v[0]), (*C.int)(&ldv),
+ (*C.double)(&u[0]), (*C.int)(&ldu),
+ (*C.int)(&nh), (*C.double)(&wh[0]), (*C.int)(&ldwh),
+ (*C.int)(&nv), (*C.double)(&wv[0]), (*C.int)(&ldwv))
+ }()
+}
diff --git a/testlapack/dlaqr5.go b/testlapack/dlaqr5.go
index 176738e1..8365a1ef 100644
--- a/testlapack/dlaqr5.go
+++ b/testlapack/dlaqr5.go
@@ -5,115 +5,210 @@
package testlapack
import (
- "archive/zip"
+ "compress/gzip"
+ "encoding/json"
"fmt"
- "io"
"log"
+ "math"
+ "math/rand"
+ "os"
+ "path/filepath"
"testing"
+
+ "github.com/gonum/blas"
+ "github.com/gonum/blas/blas64"
)
type Dlaqr5er interface {
Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, nshfts int, sr, si []float64, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int, u []float64, ldu int, nh int, wh []float64, ldwh int, nv int, wv []float64, ldwv int)
}
+type Dlaqr5test struct {
+ WantT bool
+ N int
+ NShifts int
+ KTop, KBot int
+ ShiftR, ShiftI []float64
+ H []float64
+
+ HWant []float64
+ ZWant []float64
+}
+
func Dlaqr5Test(t *testing.T, impl Dlaqr5er) {
- r, err := zip.OpenReader("../internal/testdata/dlaqr5test/dlaqr5data.zip")
+ // Test without using reference data.
+ rnd := rand.New(rand.NewSource(1))
+ for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 30} {
+ for _, extra := range []int{0, 1, 20} {
+ for _, kacc22 := range []int{0, 1, 2} {
+ for cas := 0; cas < 100; cas++ {
+ testDlaqr5(t, impl, n, extra, kacc22, rnd)
+ }
+ }
+ }
+ }
+
+ // Test using reference data computed by the reference netlib
+ // implementation.
+ file, err := os.Open(filepath.FromSlash("../testlapack/testdata/dlaqr5data.json.gz"))
+ if err != nil {
+ log.Fatal(err)
+ }
+ defer file.Close()
+ r, err := gzip.NewReader(file)
if err != nil {
log.Fatal(err)
}
defer r.Close()
- for _, f := range r.File {
- tc, err := f.Open()
- if err != nil {
- log.Fatal(err)
- }
- wantt, n, nshfts, ktop, kbot, sr, si, hOrig, hwant, zwant := readDlaqr5Case(tc)
- tc.Close()
+ var tests []Dlaqr5test
+ json.NewDecoder(r).Decode(&tests)
+ for _, test := range tests {
+ wantt := test.WantT
+ n := test.N
+ nshfts := test.NShifts
+ ktop := test.KTop
+ kbot := test.KBot
+ sr := test.ShiftR
+ si := test.ShiftI
- v := make([]float64, nshfts/2*3)
- u := make([]float64, (3*nshfts-3)*(3*nshfts-3))
- nh := n
- wh := make([]float64, (3*nshfts-3)*n)
- nv := n
- wv := make([]float64, n*(3*nshfts-3))
+ for _, extra := range []int{0, 1, 10} {
+ v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
+ u := randomGeneral(3*nshfts-3, 3*nshfts-3, 3*nshfts-3+extra, rnd)
+ nh := n
+ wh := randomGeneral(3*nshfts-3, n, n+extra, rnd)
+ nv := n
+ wv := randomGeneral(n, 3*nshfts-3, 3*nshfts-3+extra, rnd)
+
+ h := nanGeneral(n, n, n+extra)
- for _, ldh := range []int{n, n + 1, n + 10} {
- h := make([]float64, (n-1)*ldh+n)
for _, kacc22 := range []int{0, 1, 2} {
- copyMatrix(n, n, h, ldh, hOrig)
- z := eye(n, ldh)
+ copyMatrix(n, n, h.Data, h.Stride, test.H)
+ z := eye(n, n+extra)
impl.Dlaqr5(wantt, true, kacc22,
n, ktop, kbot,
nshfts, sr, si,
- h, ldh,
- 0, n-1, z, ldh,
- v, 3,
- u, 3*nshfts-3,
- nh, wh, nh,
- nv, wv, 3*nshfts-3)
+ h.Data, h.Stride,
+ 0, n-1, z.Data, z.Stride,
+ v.Data, v.Stride,
+ u.Data, u.Stride,
+ nh, wh.Data, wh.Stride,
+ nv, wv.Data, wv.Stride)
- if !equalApprox(n, n, h, ldh, hwant, 1e-13) {
- t.Errorf("Case %v, kacc22=%v: unexpected matrix H\nh =%v\nhwant=%v", f.Name, kacc22, h, hwant)
+ prefix := fmt.Sprintf("wantt=%v, n=%v, nshfts=%v, ktop=%v, kbot=%v, extra=%v, kacc22=%v",
+ wantt, n, nshfts, ktop, kbot, extra, kacc22)
+ if !equalApprox(n, n, h.Data, h.Stride, test.HWant, 1e-13) {
+ t.Errorf("Case %v: unexpected matrix H\nh =%v\nhwant=%v", prefix, h.Data, test.HWant)
}
- if !equalApprox(n, n, z, ldh, zwant, 1e-13) {
- t.Errorf("Case %v, kacc22=%v: unexpected matrix Z\nz =%v\nzwant=%v", f.Name, kacc22, z, zwant)
+ if !equalApprox(n, n, z.Data, z.Stride, test.ZWant, 1e-13) {
+ t.Errorf("Case %v: unexpected matrix Z\nz =%v\nzwant=%v", prefix, z.Data, test.ZWant)
}
}
}
}
}
-// readDlaqr5Case reads and returns test data written by internal/testdata/dlaqr5test/main.go.
-func readDlaqr5Case(r io.Reader) (wantt bool, n, nshfts, ktop, kbot int, sr, si []float64, h, hwant, zwant []float64) {
- _, err := fmt.Fscanln(r, &wantt, &n, &nshfts, &ktop, &kbot)
- if err != nil {
- log.Fatal(err)
- }
-
- sr = make([]float64, nshfts)
- si = make([]float64, nshfts)
- h = make([]float64, n*n)
- hwant = make([]float64, n*n)
- zwant = make([]float64, n*n)
-
- for i := range sr {
- _, err = fmt.Fscanln(r, &sr[i])
- if err != nil {
- log.Fatal(err)
- }
- }
- for i := range si {
- _, err = fmt.Fscanln(r, &si[i])
- if err != nil {
- log.Fatal(err)
- }
- }
+func testDlaqr5(t *testing.T, impl Dlaqr5er, n, extra, kacc22 int, rnd *rand.Rand) {
+ wantt := true
+ wantz := true
+ nshfts := 2 * n
+ sr := make([]float64, nshfts)
+ si := make([]float64, nshfts)
for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- _, err = fmt.Fscanln(r, &h[i*n+j])
- if err != nil {
- log.Fatal(err)
+ re := rnd.NormFloat64()
+ im := rnd.NormFloat64()
+ sr[2*i], sr[2*i+1] = re, re
+ si[2*i], si[2*i+1] = im, -im
+ }
+ ktop := rnd.Intn(n)
+ kbot := rnd.Intn(n)
+ if kbot < ktop {
+ ktop, kbot = kbot, ktop
+ }
+
+ v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
+ u := randomGeneral(3*nshfts-3, 3*nshfts-3, 3*nshfts-3+extra, rnd)
+ nh := n
+ wh := randomGeneral(3*nshfts-3, n, n+extra, rnd)
+ nv := n
+ wv := randomGeneral(n, 3*nshfts-3, 3*nshfts-3+extra, rnd)
+
+ h := randomHessenberg(n, n+extra, rnd)
+ if ktop > 0 {
+ h.Data[ktop*h.Stride+ktop-1] = 0
+ }
+ if kbot < n-1 {
+ h.Data[(kbot+1)*h.Stride+kbot] = 0
+ }
+ hCopy := h
+ hCopy.Data = make([]float64, len(h.Data))
+ copy(hCopy.Data, h.Data)
+
+ z := eye(n, n+extra)
+
+ impl.Dlaqr5(wantt, wantz, kacc22,
+ n, ktop, kbot,
+ nshfts, sr, si,
+ h.Data, h.Stride,
+ 0, n-1, z.Data, z.Stride,
+ v.Data, v.Stride,
+ u.Data, u.Stride,
+ nh, wh.Data, wh.Stride,
+ nv, wv.Data, wv.Stride)
+
+ prefix := fmt.Sprintf("Case n=%v, extra=%v, kacc22=%v", n, extra, kacc22)
+
+ if !generalOutsideAllNaN(h) {
+ t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
+ }
+ if !generalOutsideAllNaN(z) {
+ t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
+ }
+ if !generalOutsideAllNaN(u) {
+ t.Errorf("%v: out-of-range write to U\n%v", prefix, u.Data)
+ }
+ if !generalOutsideAllNaN(v) {
+ t.Errorf("%v: out-of-range write to V\n%v", prefix, v.Data)
+ }
+ if !generalOutsideAllNaN(wh) {
+ t.Errorf("%v: out-of-range write to WH\n%v", prefix, wh.Data)
+ }
+ if !generalOutsideAllNaN(wv) {
+ t.Errorf("%v: out-of-range write to WV\n%v", prefix, wv.Data)
+ }
+
+ for i := 0; i < n; i++ {
+ for j := 0; j < i-1; j++ {
+ if h.Data[i*h.Stride+j] != 0 {
+ t.Errorf("%v: H is not Hessenberg, H[%v,%v]!=0", prefix, i, j)
}
}
}
+ if !isOrthonormal(z) {
+ t.Errorf("%v: Z is not orthogonal", prefix)
+ }
+ // Construct Z^T * HOrig * Z and check that it is equal to H from Dlaqr5.
+ hz := blas64.General{
+ Rows: n,
+ Cols: n,
+ Stride: n,
+ Data: make([]float64, n*n),
+ }
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hCopy, z, 0, hz)
+ zhz := blas64.General{
+ Rows: n,
+ Cols: n,
+ Stride: n,
+ Data: make([]float64, n*n),
+ }
+ blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hz, 0, zhz)
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
- _, err = fmt.Fscanln(r, &hwant[i*n+j])
- if err != nil {
- log.Fatal(err)
+ diff := zhz.Data[i*zhz.Stride+j] - h.Data[i*h.Stride+j]
+ if math.Abs(diff) > 1e-13 {
+ t.Errorf("%v: Z^T*HOrig*Z and H are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
}
}
}
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- _, err = fmt.Fscanln(r, &zwant[i*n+j])
- if err != nil {
- log.Fatal(err)
- }
- }
- }
-
- return wantt, n, nshfts, ktop, kbot, sr, si, h, hwant, zwant
}
diff --git a/testlapack/general.go b/testlapack/general.go
index 5d4dbcc8..e5639ee6 100644
--- a/testlapack/general.go
+++ b/testlapack/general.go
@@ -60,6 +60,22 @@ func randomGeneral(r, c, stride int, rnd *rand.Rand) blas64.General {
return ans
}
+// randomHessenberg allocates a new n×n Hessenberg matrix filled with zeros
+// under the first subdiagonal and with random numbers elsewhere. Out-of-range
+// elements are filled with NaN values.
+func randomHessenberg(n, stride int, rnd *rand.Rand) blas64.General {
+ ans := nanGeneral(n, n, stride)
+ for i := 0; i < n; i++ {
+ for j := 0; j < i-1; j++ {
+ ans.Data[i*ans.Stride+j] = 0
+ }
+ for j := max(0, i-1); j < n; j++ {
+ ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
+ }
+ }
+ return ans
+}
+
// nanTriangular allocates a new r×c triangular matrix filled with NaN values.
func nanTriangular(uplo blas.Uplo, n, stride int) blas64.Triangular {
return blas64.Triangular{
@@ -751,11 +767,14 @@ func equalApproxTriangular(upper bool, n int, a []float64, lda int, b []float64,
return true
}
-// eye returns an identity matrix of order n and stride ld.
-func eye(n, ld int) []float64 {
- m := make([]float64, (n-1)*ld+n)
+// eye returns an identity matrix of given order and stride.
+func eye(n, stride int) blas64.General {
+ ans := nanGeneral(n, n, stride)
for i := 0; i < n; i++ {
- m[i*ld+i] = 1
+ for j := 0; j < n; j++ {
+ ans.Data[i*ans.Stride+j] = 0
+ }
+ ans.Data[i*ans.Stride+i] = 1
}
- return m
+ return ans
}
diff --git a/testlapack/testdata/dlaqr5data.json.gz b/testlapack/testdata/dlaqr5data.json.gz
new file mode 100644
index 00000000..9e629480
Binary files /dev/null and b/testlapack/testdata/dlaqr5data.json.gz differ