package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-25 08:10:28 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'dlaqr5-test'

Browse Source
This commit is contained in:
Vladimir Chalupecky
2016-06-02 15:37:06 +09:00
parent 2cc9c872fa 2a708f2c6e
commit e47d94d748
21 changed files with 301 additions and 2059 deletions
Download Patch File Download Diff File Expand all files Collapse all files

384
internal/testdata/dlaqr5test/dgemm.f vendored
View File

@@ -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

156
internal/testdata/dlaqr5test/dlacpy.f vendored
View File

@@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlacpy.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlacpy.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlacpy.f">
*> [TXT]</a>
*> \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

189
internal/testdata/dlaqr5test/dlamch.f vendored
View File

@@ -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
*
************************************************************************

104
internal/testdata/dlaqr5test/dlapy2.f vendored
View File

@@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlapy2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlapy2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlapy2.f">
*> [TXT]</a>
*> \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

BIN
internal/testdata/dlaqr5test/dlaqr5data.zip vendored
View File

Binary file not shown.

196
internal/testdata/dlaqr5test/dlarfg.f vendored
View File

@@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfg.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfg.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfg.f">
*> [TXT]</a>
*> \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

112
internal/testdata/dlaqr5test/dnrm2.f vendored
View File

@@ -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

110
internal/testdata/dlaqr5test/dscal.f vendored
View File

@@ -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

415
internal/testdata/dlaqr5test/dtrmm.f vendored
View File

@@ -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

125
internal/testdata/dlaqr5test/lsame.f vendored
View File

@@ -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

163
internal/testdata/dlaqr5test/main.go vendored
View File

@@ -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 {

5
internal/testdata/dlaqr5test/run.bash vendored
View File

@@ -1,5 +0,0 @@
#!/usr/bin/env bash
go build
./dlaqr5test
rm -f dlaqr5test

89
internal/testdata/dlaqr5test/xerbla.f vendored
View File

@@ -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

0
internal/testdata/dlaqr5test/dlabad.f → internal/testdata/netlib/dlabad.f vendored
View File

0
internal/testdata/dlaqr5test/dlaqr1.f → internal/testdata/netlib/dlaqr1.f vendored
View File

0
internal/testdata/dlaqr5test/dlaqr5.f → internal/testdata/netlib/dlaqr5.f vendored
View File

0
internal/testdata/dlaqr5test/dlaset.f → internal/testdata/netlib/dlaset.f vendored
View File

44
internal/testdata/netlib/netlib.go vendored
View File

@@ -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))
}()
}

239
testlapack/dlaqr5.go
View File

@@ -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
}

29
testlapack/general.go
View File

@@ -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
}

BIN
testlapack/testdata/dlaqr5data.json.gz vendored Normal file
View File

Binary file not shown.