Add new tests for Householder reconstruction functions from 3.9.1

This commit is contained in:
Martin Kroeker 2021-05-02 19:28:21 +02:00 committed by GitHub
parent 4c1d47098b
commit 88b70fba3e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 2033 additions and 134 deletions

View File

@ -40,7 +40,7 @@ set(SLINTST schkaa.f
sgennd.f sqrt04.f sqrt05.f schkqrt.f serrqrt.f schkqrtp.f serrqrtp.f
schklqt.f schklqtp.f schktsqr.f
serrlqt.f serrlqtp.f serrtsqr.f stsqr01.f slqt04.f slqt05.f
schkorhr_col.f serrorhr_col.f sorhr_col01.f)
schkorhr_col.f serrorhr_col.f sorhr_col01.f sorhr_col02.f)
if(USE_XBLAS)
list(APPEND SLINTST sdrvgbx.f sdrvgex.f sdrvsyx.f sdrvpox.f
@ -96,7 +96,7 @@ set(CLINTST cchkaa.f
cqrt04.f cqrt05.f cchkqrt.f cerrqrt.f cchkqrtp.f cerrqrtp.f
cchklqt.f cchklqtp.f cchktsqr.f
cerrlqt.f cerrlqtp.f cerrtsqr.f ctsqr01.f clqt04.f clqt05.f
cchkunhr_col.f cerrunhr_col.f cunhr_col01.f)
cchkunhr_col.f cerrunhr_col.f cunhr_col01.f cunhr_col02.f)
if(USE_XBLAS)
list(APPEND CLINTST cdrvgbx.f cdrvgex.f cdrvhex.f cdrvsyx.f cdrvpox.f
@ -142,7 +142,7 @@ set(DLINTST dchkaa.f
dqrt04.f dqrt05.f dchkqrt.f derrqrt.f dchkqrtp.f derrqrtp.f
dchklq.f dchklqt.f dchklqtp.f dchktsqr.f
derrlqt.f derrlqtp.f derrtsqr.f dtsqr01.f dlqt04.f dlqt05.f
dchkorhr_col.f derrorhr_col.f dorhr_col01.f)
dchkorhr_col.f derrorhr_col.f dorhr_col01.f dorhr_col02.f)
if(USE_XBLAS)
list(APPEND DLINTST ddrvgbx.f ddrvgex.f ddrvsyx.f ddrvpox.f
@ -198,7 +198,7 @@ set(ZLINTST zchkaa.f
zqrt04.f zqrt05.f zchkqrt.f zerrqrt.f zchkqrtp.f zerrqrtp.f
zchklqt.f zchklqtp.f zchktsqr.f
zerrlqt.f zerrlqtp.f zerrtsqr.f ztsqr01.f zlqt04.f zlqt05.f
zchkunhr_col.f zerrunhr_col.f zunhr_col01.f)
zchkunhr_col.f zerrunhr_col.f zunhr_col01.f zunhr_col02.f)
if(USE_XBLAS)
list(APPEND ZLINTST zdrvgbx.f zdrvgex.f zdrvhex.f zdrvsyx.f zdrvpox.f

View File

@ -74,7 +74,7 @@ SLINTST = schkaa.o \
sgennd.o sqrt04.o sqrt05.o schkqrt.o serrqrt.o schkqrtp.o serrqrtp.o \
schklqt.o schklqtp.o schktsqr.o \
serrlqt.o serrlqtp.o serrtsqr.o stsqr01.o slqt04.o slqt05.o \
schkorhr_col.o serrorhr_col.o sorhr_col01.o
schkorhr_col.o serrorhr_col.o sorhr_col01.o sorhr_col02.o
ifdef USEXBLAS
SLINTST += sdrvgbx.o sdrvgex.o sdrvsyx.o sdrvpox.o \
@ -123,7 +123,7 @@ CLINTST = cchkaa.o \
cqrt04.o cqrt05.o cchkqrt.o cerrqrt.o cchkqrtp.o cerrqrtp.o \
cchklqt.o cchklqtp.o cchktsqr.o \
cerrlqt.o cerrlqtp.o cerrtsqr.o ctsqr01.o clqt04.o clqt05.o \
cchkunhr_col.o cerrunhr_col.o cunhr_col01.o
cchkunhr_col.o cerrunhr_col.o cunhr_col01.o cunhr_col02.o
ifdef USEXBLAS
CLINTST += cdrvgbx.o cdrvgex.o cdrvhex.o cdrvsyx.o cdrvpox.o \
@ -167,7 +167,7 @@ DLINTST = dchkaa.o \
dqrt04.o dqrt05.o dchkqrt.o derrqrt.o dchkqrtp.o derrqrtp.o \
dchklq.o dchklqt.o dchklqtp.o dchktsqr.o \
derrlqt.o derrlqtp.o derrtsqr.o dtsqr01.o dlqt04.o dlqt05.o \
dchkorhr_col.o derrorhr_col.o dorhr_col01.o
dchkorhr_col.o derrorhr_col.o dorhr_col01.o dorhr_col02.o
ifdef USEXBLAS
DLINTST += ddrvgbx.o ddrvgex.o ddrvsyx.o ddrvpox.o \
@ -215,7 +215,7 @@ ZLINTST = zchkaa.o \
zqrt04.o zqrt05.o zchkqrt.o zerrqrt.o zchkqrtp.o zerrqrtp.o \
zchklqt.o zchklqtp.o zchktsqr.o \
zerrlqt.o zerrlqtp.o zerrtsqr.o ztsqr01.o zlqt04.o zlqt05.o \
zchkunhr_col.o zerrunhr_col.o zunhr_col01.o
zchkunhr_col.o zerrunhr_col.o zunhr_col01.o zunhr_col02.o
ifdef USEXBLAS
ZLINTST += zdrvgbx.o zdrvgex.o zdrvhex.o zdrvsyx.o zdrvpox.o \

View File

@ -24,9 +24,12 @@
*>
*> \verbatim
*>
*> CCHKUNHR_COL tests CUNHR_COL using CLATSQR and CGEMQRT. Therefore, CLATSQR
*> (used in CGEQR) and CGEMQRT (used in CGEMQR) have to be tested
*> before this test.
*> CCHKUNHR_COL tests:
*> 1) CUNGTSQR and CUNHR_COL using CLATSQR, CGEMQRT,
*> 2) CUNGTSQR_ROW and CUNHR_COL inside CGETSQRHRT
*> (which calls CLATSQR, CUNGTSQR_ROW and CUNHR_COL) using CGEMQRT.
*> Therefore, CLATSQR (part of CGEQR), CGEMQRT (part of CGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
@ -97,19 +100,16 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup complex_lin
*
* =====================================================================
SUBROUTINE CCHKUNHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
$ NBVAL, NOUT )
SUBROUTINE CCHKUNHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL,
$ NNB, NBVAL, NOUT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.7.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
* .. Scalar Arguments ..
LOGICAL TSTERR
@ -135,10 +135,11 @@
REAL RESULT( NTESTS )
* ..
* .. External Subroutines ..
EXTERNAL ALAHD, ALASUM, CERRUNHR_COL, CUNHR_COL01
EXTERNAL ALAHD, ALASUM, CERRUNHR_COL, CUNHR_COL01,
$ CUNHR_COL02
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
INTRINSIC MAX, MIN
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
@ -201,8 +202,8 @@
*
* Test CUNHR_COL
*
CALL CUNHR_COL01( M, N, MB1, NB1, NB2,
$ RESULT )
CALL CUNHR_COL01( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
@ -226,12 +227,78 @@
END DO
END DO
*
* Do for each value of M in MVAL.
*
DO I = 1, NM
M = MVAL( I )
*
* Do for each value of N in NVAL.
*
DO J = 1, NN
N = NVAL( J )
*
* Only for M >= N
*
IF ( MIN( M, N ).GT.0 .AND. M.GE.N ) THEN
*
* Do for each possible value of MB1
*
DO IMB1 = 1, NNB
MB1 = NBVAL( IMB1 )
*
* Only for MB1 > N
*
IF ( MB1.GT.N ) THEN
*
* Do for each possible value of NB1
*
DO INB1 = 1, NNB
NB1 = NBVAL( INB1 )
*
* Do for each possible value of NB2
*
DO INB2 = 1, NNB
NB2 = NBVAL( INB2 )
*
IF( NB1.GT.0 .AND. NB2.GT.0 ) THEN
*
* Test CUNHR_COL
*
CALL CUNHR_COL02( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
*
DO T = 1, NTESTS
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 ) M, N, MB1,
$ NB1, NB2, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + NTESTS
END IF
END DO
END DO
END IF
END DO
END IF
END DO
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 'M=', I5, ', N=', I5, ', MB1=', I5,
$ ', NB1=', I5, ', NB2=', I5,' test(', I2, ')=', G12.5 )
9999 FORMAT( 'CUNGTSQR and CUNHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
9998 FORMAT( 'CUNGTSQR_ROW and CUNHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
RETURN
*
* End of CCHKUNHR_COL

View File

@ -13,7 +13,7 @@
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* REAL RESULT(6)
* DOUBLE PRECISION RESULT(6)
*
*
*> \par Purpose:
@ -21,8 +21,8 @@
*>
*> \verbatim
*>
*> CUNHR_COL01 tests CUNHR_COL using CLATSQR, CGEMQRT and CUNGTSQR.
*> Therefore, CLATSQR (part of CGEQR), CGEMQRT (part CGEMQR), CUNGTSQR
*> CUNHR_COL01 tests CUNGTSQR and CUNHR_COL using CLATSQR, CGEMQRT.
*> Therefore, CLATSQR (part of CGEQR), CGEMQRT (part of CGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
@ -62,14 +62,46 @@
*> \verbatim
*> RESULT is REAL array, dimension (6)
*> Results of each of the six tests below.
*> ( C is a M-by-N random matrix, D is a N-by-M random matrix )
*>
*> RESULT(1) = | A - Q * R | / (eps * m * |A|)
*> RESULT(2) = | I - (Q**H) * Q | / (eps * m )
*> RESULT(3) = | Q * C - Q * C | / (eps * m * |C|)
*> RESULT(4) = | (Q**H) * C - (Q**H) * C | / (eps * m * |C|)
*> RESULT(5) = | (D * Q) - D * Q | / (eps * m * |D|)
*> RESULT(6) = | D * (Q**H) - D * (Q**H) | / (eps * m * |D|)
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m unitary Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in CGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m unitary matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using CGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using CGEMM.
*> \endverbatim
*
* Authors:
@ -80,18 +112,15 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup complex16_lin
*> \ingroup complex_lin
*
* =====================================================================
SUBROUTINE CUNHR_COL01( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.9.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
@ -102,10 +131,10 @@
*
* ..
* .. Local allocatable arrays
COMPLEX, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
COMPLEX , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
REAL, ALLOCATABLE :: RWORK(:)
REAL , ALLOCATABLE :: RWORK(:)
*
* .. Parameters ..
REAL ZERO
@ -218,7 +247,7 @@
* Copy the factor R into the array R.
*
SRNAMT = 'CLACPY'
CALL CLACPY( 'U', M, N, AF, M, R, M )
CALL CLACPY( 'U', N, N, AF, M, R, M )
*
* Reconstruct the orthogonal matrix Q.
*
@ -240,7 +269,7 @@
* matrix S.
*
SRNAMT = 'CLACPY'
CALL CLACPY( 'U', M, N, R, M, AF, M )
CALL CLACPY( 'U', N, N, R, M, AF, M )
*
DO I = 1, N
IF( DIAG( I ).EQ.-CONE ) THEN

View File

@ -0,0 +1,381 @@
*> \brief \b CUNHR_COL02
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE CUNHR_COL02( M, N, MB1, NB1, NB2, RESULT )
*
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* REAL RESULT(6)
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CUNHR_COL02 tests CUNGTSQR_ROW and CUNHR_COL inside CGETSQRHRT
*> (which calls CLATSQR, CUNGTSQR_ROW and CUNHR_COL) using CGEMQRT.
*> Therefore, CLATSQR (part of CGEQR), CGEMQRT (part of CGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> Number of rows in test matrix.
*> \endverbatim
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> Number of columns in test matrix.
*> \endverbatim
*> \param[in] MB1
*> \verbatim
*> MB1 is INTEGER
*> Number of row in row block in an input test matrix.
*> \endverbatim
*>
*> \param[in] NB1
*> \verbatim
*> NB1 is INTEGER
*> Number of columns in column block an input test matrix.
*> \endverbatim
*>
*> \param[in] NB2
*> \verbatim
*> NB2 is INTEGER
*> Number of columns in column block in an output test matrix.
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is REAL array, dimension (6)
*> Results of each of the six tests below.
*>
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m unitary Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in CGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m unitary matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using CGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using CGEMM.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex_lin
*
* =====================================================================
SUBROUTINE CUNHR_COL02( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
REAL RESULT(6)
*
* =====================================================================
*
* ..
* .. Local allocatable arrays
COMPLEX , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
REAL , ALLOCATABLE :: RWORK(:)
*
* .. Parameters ..
REAL ZERO
PARAMETER ( ZERO = 0.0E+0 )
COMPLEX CONE, CZERO
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ),
$ CZERO = ( 0.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
LOGICAL TESTZEROS
INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
REAL ANORM, EPS, RESID, CNORM, DNORM
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 )
COMPLEX WORKQUERY( 1 )
* ..
* .. External Functions ..
REAL SLAMCH, CLANGE, CLANSY
EXTERNAL SLAMCH, CLANGE, CLANSY
* ..
* .. External Subroutines ..
EXTERNAL CLACPY, CLARNV, CLASET, CGETSQRHRT,
$ CSCAL, CGEMM, CGEMQRT, CHERK
* ..
* .. Intrinsic Functions ..
INTRINSIC CEILING, REAL, MAX, MIN
* ..
* .. Scalars in Common ..
CHARACTER(LEN=32) SRNAMT
* ..
* .. Common blocks ..
COMMON / SRMNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEED / 1988, 1989, 1990, 1991 /
*
* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
*
TESTZEROS = .FALSE.
*
EPS = SLAMCH( 'Epsilon' )
K = MIN( M, N )
L = MAX( M, N, 1)
*
* Dynamically allocate local arrays
*
ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
$ C(M,N), CF(M,N),
$ D(N,M), DF(N,M) )
*
* Put random numbers into A and copy to AF
*
DO J = 1, N
CALL CLARNV( 2, ISEED, M, A( 1, J ) )
END DO
IF( TESTZEROS ) THEN
IF( M.GE.4 ) THEN
DO J = 1, N
CALL CLARNV( 2, ISEED, M/2, A( M/4, J ) )
END DO
END IF
END IF
CALL CLACPY( 'Full', M, N, A, M, AF, M )
*
* Number of row blocks in CLATSQR
*
NRB = MAX( 1, CEILING( REAL( M - N ) / REAL( MB1 - N ) ) )
*
ALLOCATE ( T1( NB1, N * NRB ) )
ALLOCATE ( T2( NB2, N ) )
ALLOCATE ( DIAG( N ) )
*
* Begin determine LWORK for the array WORK and allocate memory.
*
* CGEMQRT requires NB2 to be bounded by N.
*
NB2_UB = MIN( NB2, N)
*
*
CALL CGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORKQUERY, -1, INFO )
*
LWORK = INT( WORKQUERY( 1 ) )
*
* In CGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
* or M*NB2_UB if SIDE = 'R'.
*
LWORK = MAX( LWORK, NB2_UB * N, NB2_UB * M )
*
ALLOCATE ( WORK( LWORK ) )
*
* End allocate memory for WORK.
*
*
* Begin Householder reconstruction routines
*
* Factor the matrix A in the array AF.
*
SRNAMT = 'CGETSQRHRT'
CALL CGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORK, LWORK, INFO )
*
* End Householder reconstruction routines.
*
*
* Generate the m-by-m matrix Q
*
CALL CLASET( 'Full', M, M, CZERO, CONE, Q, M )
*
SRNAMT = 'CGEMQRT'
CALL CGEMQRT( 'L', 'N', M, M, K, NB2_UB, AF, M, T2, NB2, Q, M,
$ WORK, INFO )
*
* Copy R
*
CALL CLASET( 'Full', M, N, CZERO, CZERO, R, M )
*
CALL CLACPY( 'Upper', M, N, AF, M, R, M )
*
* TEST 1
* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
*
CALL CGEMM( 'C', 'N', M, N, M, -CONE, Q, M, A, M, CONE, R, M )
*
ANORM = CLANGE( '1', M, N, A, M, RWORK )
RESID = CLANGE( '1', M, N, R, M, RWORK )
IF( ANORM.GT.ZERO ) THEN
RESULT( 1 ) = RESID / ( EPS * MAX( 1, M ) * ANORM )
ELSE
RESULT( 1 ) = ZERO
END IF
*
* TEST 2
* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
*
CALL CLASET( 'Full', M, M, CZERO, CONE, R, M )
CALL CHERK( 'U', 'C', M, M, -CONE, Q, M, CONE, R, M )
RESID = CLANSY( '1', 'Upper', M, R, M, RWORK )
RESULT( 2 ) = RESID / ( EPS * MAX( 1, M ) )
*
* Generate random m-by-n matrix C
*
DO J = 1, N
CALL CLARNV( 2, ISEED, M, C( 1, J ) )
END DO
CNORM = CLANGE( '1', M, N, C, M, RWORK )
CALL CLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as Q*C = CF
*
SRNAMT = 'CGEMQRT'
CALL CGEMQRT( 'L', 'N', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 3
* Compute |CF - Q*C| / ( eps * m * |C| )
*
CALL CGEMM( 'N', 'N', M, N, M, -CONE, Q, M, C, M, CONE, CF, M )
RESID = CLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 3 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 3 ) = ZERO
END IF
*
* Copy C into CF again
*
CALL CLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as (Q**T)*C = CF
*
SRNAMT = 'CGEMQRT'
CALL CGEMQRT( 'L', 'C', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 4
* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
*
CALL CGEMM( 'C', 'N', M, N, M, -CONE, Q, M, C, M, CONE, CF, M )
RESID = CLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 4 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 4 ) = ZERO
END IF
*
* Generate random n-by-m matrix D and a copy DF
*
DO J = 1, M
CALL CLARNV( 2, ISEED, N, D( 1, J ) )
END DO
DNORM = CLANGE( '1', N, M, D, N, RWORK )
CALL CLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*Q = DF
*
SRNAMT = 'CGEMQRT'
CALL CGEMQRT( 'R', 'N', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 5
* Compute |DF - D*Q| / ( eps * m * |D| )
*
CALL CGEMM( 'N', 'N', N, M, M, -CONE, D, N, Q, M, CONE, DF, N )
RESID = CLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 5 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 5 ) = ZERO
END IF
*
* Copy D into DF again
*
CALL CLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*QT = DF
*
SRNAMT = 'CGEMQRT'
CALL CGEMQRT( 'R', 'C', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 6
* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
*
CALL CGEMM( 'N', 'C', N, M, M, -CONE, D, N, Q, M, CONE, DF, N )
RESID = CLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 6 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 6 ) = ZERO
END IF
*
* Deallocate all arrays
*
DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T1, T2, DIAG,
$ C, D, CF, DF )
*
RETURN
*
* End of CUNHR_COL02
*
END

View File

@ -24,9 +24,12 @@
*>
*> \verbatim
*>
*> DCHKORHR_COL tests DORHR_COL using DLATSQR and DGEMQRT. Therefore, DLATSQR
*> (used in DGEQR) and DGEMQRT (used in DGEMQR) have to be tested
*> before this test.
*> DCHKORHR_COL tests:
*> 1) DORGTSQR and DORHR_COL using DLATSQR, DGEMQRT,
*> 2) DORGTSQR_ROW and DORHR_COL inside DGETSQRHRT
*> (which calls DLATSQR, DORGTSQR_ROW and DORHR_COL) using DGEMQRT.
*> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
@ -97,19 +100,16 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup double_lin
*
* =====================================================================
SUBROUTINE DCHKORHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
$ NBVAL, NOUT )
SUBROUTINE DCHKORHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL,
$ NNB, NBVAL, NOUT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.7.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
* .. Scalar Arguments ..
LOGICAL TSTERR
@ -135,10 +135,11 @@
DOUBLE PRECISION RESULT( NTESTS )
* ..
* .. External Subroutines ..
EXTERNAL ALAHD, ALASUM, DERRORHR_COL, DORHR_COL01
EXTERNAL ALAHD, ALASUM, DERRORHR_COL, DORHR_COL01,
$ DORHR_COL02
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
INTRINSIC MAX, MIN
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
@ -201,8 +202,8 @@
*
* Test DORHR_COL
*
CALL DORHR_COL01( M, N, MB1, NB1, NB2,
$ RESULT )
CALL DORHR_COL01( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
@ -226,12 +227,78 @@
END DO
END DO
*
* Do for each value of M in MVAL.
*
DO I = 1, NM
M = MVAL( I )
*
* Do for each value of N in NVAL.
*
DO J = 1, NN
N = NVAL( J )
*
* Only for M >= N
*
IF ( MIN( M, N ).GT.0 .AND. M.GE.N ) THEN
*
* Do for each possible value of MB1
*
DO IMB1 = 1, NNB
MB1 = NBVAL( IMB1 )
*
* Only for MB1 > N
*
IF ( MB1.GT.N ) THEN
*
* Do for each possible value of NB1
*
DO INB1 = 1, NNB
NB1 = NBVAL( INB1 )
*
* Do for each possible value of NB2
*
DO INB2 = 1, NNB
NB2 = NBVAL( INB2 )
*
IF( NB1.GT.0 .AND. NB2.GT.0 ) THEN
*
* Test DORHR_COL
*
CALL DORHR_COL02( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
*
DO T = 1, NTESTS
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 ) M, N, MB1,
$ NB1, NB2, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + NTESTS
END IF
END DO
END DO
END IF
END DO
END IF
END DO
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 'M=', I5, ', N=', I5, ', MB1=', I5,
$ ', NB1=', I5, ', NB2=', I5,' test(', I2, ')=', G12.5 )
9999 FORMAT( 'DORGTSQR and DORHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
9998 FORMAT( 'DORGTSQR_ROW and DORHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
RETURN
*
* End of DCHKORHR_COL

View File

@ -21,8 +21,8 @@
*>
*> \verbatim
*>
*> DORHR_COL01 tests DORHR_COL using DLATSQR, DGEMQRT and DORGTSQR.
*> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part DGEMQR), DORGTSQR
*> DORHR_COL01 tests DORGTSQR and DORHR_COL using DLATSQR, DGEMQRT.
*> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
@ -62,14 +62,46 @@
*> \verbatim
*> RESULT is DOUBLE PRECISION array, dimension (6)
*> Results of each of the six tests below.
*> ( C is a M-by-N random matrix, D is a N-by-M random matrix )
*>
*> RESULT(1) = | A - Q * R | / (eps * m * |A|)
*> RESULT(2) = | I - (Q**H) * Q | / (eps * m )
*> RESULT(3) = | Q * C - Q * C | / (eps * m * |C|)
*> RESULT(4) = | (Q**H) * C - (Q**H) * C | / (eps * m * |C|)
*> RESULT(5) = | (D * Q) - D * Q | / (eps * m * |D|)
*> RESULT(6) = | D * (Q**H) - D * (Q**H) | / (eps * m * |D|)
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m orthogonal Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in ZGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using DGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using DGEMM.
*> \endverbatim
*
* Authors:
@ -80,18 +112,15 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup single_lin
*> \ingroup double_lin
*
* =====================================================================
SUBROUTINE DORHR_COL01( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.9.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2

View File

@ -0,0 +1,377 @@
*> \brief \b DORHR_COL02
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
*
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* DOUBLE PRECISION RESULT(6)
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DORHR_COL02 tests DORGTSQR_ROW and DORHR_COL inside DGETSQRHRT
*> (which calls DLATSQR, DORGTSQR_ROW and DORHR_COL) using DGEMQRT.
*> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> Number of rows in test matrix.
*> \endverbatim
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> Number of columns in test matrix.
*> \endverbatim
*> \param[in] MB1
*> \verbatim
*> MB1 is INTEGER
*> Number of row in row block in an input test matrix.
*> \endverbatim
*>
*> \param[in] NB1
*> \verbatim
*> NB1 is INTEGER
*> Number of columns in column block an input test matrix.
*> \endverbatim
*>
*> \param[in] NB2
*> \verbatim
*> NB2 is INTEGER
*> Number of columns in column block in an output test matrix.
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is DOUBLE PRECISION array, dimension (6)
*> Results of each of the six tests below.
*>
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m orthogonal Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in ZGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using DGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using DGEMM.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_lin
*
* =====================================================================
SUBROUTINE DORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
DOUBLE PRECISION RESULT(6)
*
* =====================================================================
*
* ..
* .. Local allocatable arrays
DOUBLE PRECISION, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL TESTZEROS
INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 )
DOUBLE PRECISION WORKQUERY( 1 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
EXTERNAL DLAMCH, DLANGE, DLANSY
* ..
* .. External Subroutines ..
EXTERNAL DLACPY, DLARNV, DLASET, DGETSQRHRT,
$ DSCAL, DGEMM, DGEMQRT, DSYRK
* ..
* .. Intrinsic Functions ..
INTRINSIC CEILING, DBLE, MAX, MIN
* ..
* .. Scalars in Common ..
CHARACTER(LEN=32) SRNAMT
* ..
* .. Common blocks ..
COMMON / SRMNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEED / 1988, 1989, 1990, 1991 /
*
* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
*
TESTZEROS = .FALSE.
*
EPS = DLAMCH( 'Epsilon' )
K = MIN( M, N )
L = MAX( M, N, 1)
*
* Dynamically allocate local arrays
*
ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
$ C(M,N), CF(M,N),
$ D(N,M), DF(N,M) )
*
* Put random numbers into A and copy to AF
*
DO J = 1, N
CALL DLARNV( 2, ISEED, M, A( 1, J ) )
END DO
IF( TESTZEROS ) THEN
IF( M.GE.4 ) THEN
DO J = 1, N
CALL DLARNV( 2, ISEED, M/2, A( M/4, J ) )
END DO
END IF
END IF
CALL DLACPY( 'Full', M, N, A, M, AF, M )
*
* Number of row blocks in DLATSQR
*
NRB = MAX( 1, CEILING( DBLE( M - N ) / DBLE( MB1 - N ) ) )
*
ALLOCATE ( T1( NB1, N * NRB ) )
ALLOCATE ( T2( NB2, N ) )
ALLOCATE ( DIAG( N ) )
*
* Begin determine LWORK for the array WORK and allocate memory.
*
* DGEMQRT requires NB2 to be bounded by N.
*
NB2_UB = MIN( NB2, N)
*
*
CALL DGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORKQUERY, -1, INFO )
*
LWORK = INT( WORKQUERY( 1 ) )
*
* In DGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
* or M*NB2_UB if SIDE = 'R'.
*
LWORK = MAX( LWORK, NB2_UB * N, NB2_UB * M )
*
ALLOCATE ( WORK( LWORK ) )
*
* End allocate memory for WORK.
*
*
* Begin Householder reconstruction routines
*
* Factor the matrix A in the array AF.
*
SRNAMT = 'DGETSQRHRT'
CALL DGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORK, LWORK, INFO )
*
* End Householder reconstruction routines.
*
*
* Generate the m-by-m matrix Q
*
CALL DLASET( 'Full', M, M, ZERO, ONE, Q, M )
*
SRNAMT = 'DGEMQRT'
CALL DGEMQRT( 'L', 'N', M, M, K, NB2_UB, AF, M, T2, NB2, Q, M,
$ WORK, INFO )
*
* Copy R
*
CALL DLASET( 'Full', M, N, ZERO, ZERO, R, M )
*
CALL DLACPY( 'Upper', M, N, AF, M, R, M )
*
* TEST 1
* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
*
CALL DGEMM( 'T', 'N', M, N, M, -ONE, Q, M, A, M, ONE, R, M )
*
ANORM = DLANGE( '1', M, N, A, M, RWORK )
RESID = DLANGE( '1', M, N, R, M, RWORK )
IF( ANORM.GT.ZERO ) THEN
RESULT( 1 ) = RESID / ( EPS * MAX( 1, M ) * ANORM )
ELSE
RESULT( 1 ) = ZERO
END IF
*
* TEST 2
* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
*
CALL DLASET( 'Full', M, M, ZERO, ONE, R, M )
CALL DSYRK( 'U', 'T', M, M, -ONE, Q, M, ONE, R, M )
RESID = DLANSY( '1', 'Upper', M, R, M, RWORK )
RESULT( 2 ) = RESID / ( EPS * MAX( 1, M ) )
*
* Generate random m-by-n matrix C
*
DO J = 1, N
CALL DLARNV( 2, ISEED, M, C( 1, J ) )
END DO
CNORM = DLANGE( '1', M, N, C, M, RWORK )
CALL DLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as Q*C = CF
*
SRNAMT = 'DGEMQRT'
CALL DGEMQRT( 'L', 'N', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 3
* Compute |CF - Q*C| / ( eps * m * |C| )
*
CALL DGEMM( 'N', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
RESID = DLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 3 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 3 ) = ZERO
END IF
*
* Copy C into CF again
*
CALL DLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as (Q**T)*C = CF
*
SRNAMT = 'DGEMQRT'
CALL DGEMQRT( 'L', 'T', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 4
* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
*
CALL DGEMM( 'T', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
RESID = DLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 4 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 4 ) = ZERO
END IF
*
* Generate random n-by-m matrix D and a copy DF
*
DO J = 1, M
CALL DLARNV( 2, ISEED, N, D( 1, J ) )
END DO
DNORM = DLANGE( '1', N, M, D, N, RWORK )
CALL DLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*Q = DF
*
SRNAMT = 'DGEMQRT'
CALL DGEMQRT( 'R', 'N', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 5
* Compute |DF - D*Q| / ( eps * m * |D| )
*
CALL DGEMM( 'N', 'N', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
RESID = DLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 5 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 5 ) = ZERO
END IF
*
* Copy D into DF again
*
CALL DLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*QT = DF
*
SRNAMT = 'DGEMQRT'
CALL DGEMQRT( 'R', 'T', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 6
* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
*
CALL DGEMM( 'N', 'T', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
RESID = DLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 6 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 6 ) = ZERO
END IF
*
* Deallocate all arrays
*
DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T1, T2, DIAG,
$ C, D, CF, DF )
*
RETURN
*
* End of DORHR_COL02
*
END

View File

@ -24,8 +24,11 @@
*>
*> \verbatim
*>
*> SCHKORHR_COL tests SORHR_COL using SLATSQR, SGEMQRT and SORGTSQR.
*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part SGEMQR), SORGTSQR
*> SCHKORHR_COL tests:
*> 1) SORGTSQR and SORHR_COL using SLATSQR, SGEMQRT,
*> 2) SORGTSQR_ROW and SORHR_COL inside DGETSQRHRT
*> (which calls SLATSQR, SORGTSQR_ROW and SORHR_COL) using SGEMQRT.
*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part of SGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
@ -97,19 +100,16 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup sigle_lin
*> \ingroup single_lin
*
* =====================================================================
SUBROUTINE SCHKORHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
$ NBVAL, NOUT )
SUBROUTINE SCHKORHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL,
$ NNB, NBVAL, NOUT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.9.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* June 2019
*
* .. Scalar Arguments ..
LOGICAL TSTERR
@ -135,7 +135,8 @@
REAL RESULT( NTESTS )
* ..
* .. External Subroutines ..
EXTERNAL ALAHD, ALASUM, SERRORHR_COL, SORHR_COL01
EXTERNAL ALAHD, ALASUM, SERRORHR_COL, SORHR_COL01,
$ SORHR_COL02
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
@ -201,8 +202,8 @@
*
* Test SORHR_COL
*
CALL SORHR_COL01( M, N, MB1, NB1, NB2,
$ RESULT )
CALL SORHR_COL01( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
@ -226,12 +227,78 @@
END DO
END DO
*
* Do for each value of M in MVAL.
*
DO I = 1, NM
M = MVAL( I )
*
* Do for each value of N in NVAL.
*
DO J = 1, NN
N = NVAL( J )
*
* Only for M >= N
*
IF ( MIN( M, N ).GT.0 .AND. M.GE.N ) THEN
*
* Do for each possible value of MB1
*
DO IMB1 = 1, NNB
MB1 = NBVAL( IMB1 )
*
* Only for MB1 > N
*
IF ( MB1.GT.N ) THEN
*
* Do for each possible value of NB1
*
DO INB1 = 1, NNB
NB1 = NBVAL( INB1 )
*
* Do for each possible value of NB2
*
DO INB2 = 1, NNB
NB2 = NBVAL( INB2 )
*
IF( NB1.GT.0 .AND. NB2.GT.0 ) THEN
*
* Test SORHR_COL
*
CALL SORHR_COL02( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
*
DO T = 1, NTESTS
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 ) M, N, MB1,
$ NB1, NB2, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + NTESTS
END IF
END DO
END DO
END IF
END DO
END IF
END DO
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 'M=', I5, ', N=', I5, ', MB1=', I5,
$ ', NB1=', I5, ', NB2=', I5,' test(', I2, ')=', G12.5 )
9999 FORMAT( 'SORGTSQR and SORHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
9998 FORMAT( 'SORGTSQR_ROW and SORHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
RETURN
*
* End of SCHKORHR_COL

View File

@ -8,12 +8,12 @@
* Definition:
* ===========
*
* SUBROUTINE SORHR_COL01( M, N, MB1, NB1, NB2, RESULT)
* SUBROUTINE SORHR_COL01( M, N, MB1, NB1, NB2, RESULT )
*
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* REAL RESULT(6)
* REAL RESULT(6)
*
*
*> \par Purpose:
@ -21,8 +21,8 @@
*>
*> \verbatim
*>
*> SORHR_COL01 tests SORHR_COL using SLATSQR, SGEMQRT and SORGTSQR.
*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part SGEMQR), SORGTSQR
*> SORHR_COL01 tests SORGTSQR and SORHR_COL using SLATSQR, SGEMQRT.
*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part of SGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
@ -62,14 +62,46 @@
*> \verbatim
*> RESULT is REAL array, dimension (6)
*> Results of each of the six tests below.
*> ( C is a M-by-N random matrix, D is a N-by-M random matrix )
*>
*> RESULT(1) = | A - Q * R | / (eps * m * |A|)
*> RESULT(2) = | I - (Q**H) * Q | / (eps * m )
*> RESULT(3) = | Q * C - Q * C | / (eps * m * |C|)
*> RESULT(4) = | (Q**H) * C - (Q**H) * C | / (eps * m * |C|)
*> RESULT(5) = | (D * Q) - D * Q | / (eps * m * |D|)
*> RESULT(6) = | D * (Q**H) - D * (Q**H) | / (eps * m * |D|)
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m orthogonal Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in SGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using SGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using SGEMM.
*> \endverbatim
*
* Authors:
@ -80,18 +112,15 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup single_lin
*
* =====================================================================
SUBROUTINE SORHR_COL01( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.9.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
@ -102,7 +131,7 @@
*
* ..
* .. Local allocatable arrays
REAL, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
REAL , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
*
@ -128,7 +157,7 @@
$ SORGTSQR, SSCAL, SGEMM, SGEMQRT, SSYRK
* ..
* .. Intrinsic Functions ..
INTRINSIC CEILING, MAX, MIN, REAL
INTRINSIC CEILING, REAL, MAX, MIN
* ..
* .. Scalars in Common ..
CHARACTER(LEN=32) SRNAMT
@ -230,7 +259,7 @@
*
* Compute the factor R_hr corresponding to the Householder
* reconstructed Q_hr and place it in the upper triangle of AF to
* match the Q storage format in DGEQRT. R_hr = R_tsqr * S,
* match the Q storage format in SGEQRT. R_hr = R_tsqr * S,
* this means changing the sign of I-th row of the matrix R_tsqr
* according to sign of of I-th diagonal element DIAG(I) of the
* matrix S.

View File

@ -0,0 +1,376 @@
*> \brief \b SORHR_COL02
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE SORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
*
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* REAL RESULT(6)
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SORHR_COL02 tests SORGTSQR_ROW and SORHR_COL inside SGETSQRHRT
*> (which calls SLATSQR, SORGTSQR_ROW and SORHR_COL) using SGEMQRT.
*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part of SGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> Number of rows in test matrix.
*> \endverbatim
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> Number of columns in test matrix.
*> \endverbatim
*> \param[in] MB1
*> \verbatim
*> MB1 is INTEGER
*> Number of row in row block in an input test matrix.
*> \endverbatim
*>
*> \param[in] NB1
*> \verbatim
*> NB1 is INTEGER
*> Number of columns in column block an input test matrix.
*> \endverbatim
*>
*> \param[in] NB2
*> \verbatim
*> NB2 is INTEGER
*> Number of columns in column block in an output test matrix.
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is REAL array, dimension (6)
*> Results of each of the six tests below.
*>
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m orthogonal Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in SGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using SGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using SGEMM.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup single_lin
*
* =====================================================================
SUBROUTINE SORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
REAL RESULT(6)
*
* =====================================================================
*
* ..
* .. Local allocatable arrays
REAL , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
* ..
* .. Local Scalars ..
LOGICAL TESTZEROS
INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
REAL ANORM, EPS, RESID, CNORM, DNORM
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 )
REAL WORKQUERY( 1 )
* ..
* .. External Functions ..
REAL SLAMCH, SLANGE, SLANSY
EXTERNAL SLAMCH, SLANGE, SLANSY
* ..
* .. External Subroutines ..
EXTERNAL SLACPY, SLARNV, SLASET, SGETSQRHRT,
$ SSCAL, SGEMM, SGEMQRT, SSYRK
* ..
* .. Intrinsic Functions ..
INTRINSIC CEILING, REAL, MAX, MIN
* ..
* .. Scalars in Common ..
CHARACTER(LEN=32) SRNAMT
* ..
* .. Common blocks ..
COMMON / SRMNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEED / 1988, 1989, 1990, 1991 /
*
* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
*
TESTZEROS = .FALSE.
*
EPS = SLAMCH( 'Epsilon' )
K = MIN( M, N )
L = MAX( M, N, 1)
*
* Dynamically allocate local arrays
*
ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
$ C(M,N), CF(M,N),
$ D(N,M), DF(N,M) )
*
* Put random numbers into A and copy to AF
*
DO J = 1, N
CALL SLARNV( 2, ISEED, M, A( 1, J ) )
END DO
IF( TESTZEROS ) THEN
IF( M.GE.4 ) THEN
DO J = 1, N
CALL SLARNV( 2, ISEED, M/2, A( M/4, J ) )
END DO
END IF
END IF
CALL SLACPY( 'Full', M, N, A, M, AF, M )
*
* Number of row blocks in SLATSQR
*
NRB = MAX( 1, CEILING( REAL( M - N ) / REAL( MB1 - N ) ) )
*
ALLOCATE ( T1( NB1, N * NRB ) )
ALLOCATE ( T2( NB2, N ) )
ALLOCATE ( DIAG( N ) )
*
* Begin determine LWORK for the array WORK and allocate memory.
*
* SGEMQRT requires NB2 to be bounded by N.
*
NB2_UB = MIN( NB2, N)
*
CALL SGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORKQUERY, -1, INFO )
*
LWORK = INT( WORKQUERY( 1 ) )
*
* In SGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
* or M*NB2_UB if SIDE = 'R'.
*
LWORK = MAX( LWORK, NB2_UB * N, NB2_UB * M )
*
ALLOCATE ( WORK( LWORK ) )
*
* End allocate memory for WORK.
*
*
* Begin Householder reconstruction routines
*
* Factor the matrix A in the array AF.
*
SRNAMT = 'SGETSQRHRT'
CALL SGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORK, LWORK, INFO )
*
* End Householder reconstruction routines.
*
*
* Generate the m-by-m matrix Q
*
CALL SLASET( 'Full', M, M, ZERO, ONE, Q, M )
*
SRNAMT = 'SGEMQRT'
CALL SGEMQRT( 'L', 'N', M, M, K, NB2_UB, AF, M, T2, NB2, Q, M,
$ WORK, INFO )
*
* Copy R
*
CALL SLASET( 'Full', M, N, ZERO, ZERO, R, M )
*
CALL SLACPY( 'Upper', M, N, AF, M, R, M )
*
* TEST 1
* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
*
CALL SGEMM( 'T', 'N', M, N, M, -ONE, Q, M, A, M, ONE, R, M )
*
ANORM = SLANGE( '1', M, N, A, M, RWORK )
RESID = SLANGE( '1', M, N, R, M, RWORK )
IF( ANORM.GT.ZERO ) THEN
RESULT( 1 ) = RESID / ( EPS * MAX( 1, M ) * ANORM )
ELSE
RESULT( 1 ) = ZERO
END IF
*
* TEST 2
* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
*
CALL SLASET( 'Full', M, M, ZERO, ONE, R, M )
CALL SSYRK( 'U', 'T', M, M, -ONE, Q, M, ONE, R, M )
RESID = SLANSY( '1', 'Upper', M, R, M, RWORK )
RESULT( 2 ) = RESID / ( EPS * MAX( 1, M ) )
*
* Generate random m-by-n matrix C
*
DO J = 1, N
CALL SLARNV( 2, ISEED, M, C( 1, J ) )
END DO
CNORM = SLANGE( '1', M, N, C, M, RWORK )
CALL SLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as Q*C = CF
*
SRNAMT = 'SGEMQRT'
CALL SGEMQRT( 'L', 'N', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 3
* Compute |CF - Q*C| / ( eps * m * |C| )
*
CALL SGEMM( 'N', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
RESID = SLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 3 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 3 ) = ZERO
END IF
*
* Copy C into CF again
*
CALL SLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as (Q**T)*C = CF
*
SRNAMT = 'SGEMQRT'
CALL SGEMQRT( 'L', 'T', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 4
* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
*
CALL SGEMM( 'T', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
RESID = SLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 4 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 4 ) = ZERO
END IF
*
* Generate random n-by-m matrix D and a copy DF
*
DO J = 1, M
CALL SLARNV( 2, ISEED, N, D( 1, J ) )
END DO
DNORM = SLANGE( '1', N, M, D, N, RWORK )
CALL SLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*Q = DF
*
SRNAMT = 'SGEMQRT'
CALL SGEMQRT( 'R', 'N', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 5
* Compute |DF - D*Q| / ( eps * m * |D| )
*
CALL SGEMM( 'N', 'N', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
RESID = SLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 5 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 5 ) = ZERO
END IF
*
* Copy D into DF again
*
CALL SLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*QT = DF
*
SRNAMT = 'SGEMQRT'
CALL SGEMQRT( 'R', 'T', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 6
* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
*
CALL SGEMM( 'N', 'T', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
RESID = SLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 6 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 6 ) = ZERO
END IF
*
* Deallocate all arrays
*
DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T1, T2, DIAG,
$ C, D, CF, DF )
*
RETURN
*
* End of SORHR_COL02
*
END

View File

@ -24,9 +24,12 @@
*>
*> \verbatim
*>
*> ZCHKUNHR_COL tests ZUNHR_COL using ZLATSQR and ZGEMQRT. Therefore, ZLATSQR
*> (used in ZGEQR) and ZGEMQRT (used in ZGEMQR) have to be tested
*> before this test.
*> ZCHKUNHR_COL tests:
*> 1) ZUNGTSQR and ZUNHR_COL using ZLATSQR, ZGEMQRT,
*> 2) ZUNGTSQR_ROW and ZUNHR_COL inside ZGETSQRHRT
*> (which calls ZLATSQR, ZUNGTSQR_ROW and ZUNHR_COL) using ZGEMQRT.
*> Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part of ZGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
@ -97,19 +100,16 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup complex16_lin
*
* =====================================================================
SUBROUTINE ZCHKUNHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
$ NBVAL, NOUT )
SUBROUTINE ZCHKUNHR_COL( THRESH, TSTERR, NM, MVAL, NN, NVAL,
$ NNB, NBVAL, NOUT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.7.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
* .. Scalar Arguments ..
LOGICAL TSTERR
@ -135,10 +135,11 @@
DOUBLE PRECISION RESULT( NTESTS )
* ..
* .. External Subroutines ..
EXTERNAL ALAHD, ALASUM, ZERRUNHR_COL, ZUNHR_COL01
EXTERNAL ALAHD, ALASUM, ZERRUNHR_COL, ZUNHR_COL01,
$ ZUNHR_COL02
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
INTRINSIC MAX, MIN
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
@ -201,8 +202,8 @@
*
* Test ZUNHR_COL
*
CALL ZUNHR_COL01( M, N, MB1, NB1, NB2,
$ RESULT )
CALL ZUNHR_COL01( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
@ -226,12 +227,78 @@
END DO
END DO
*
* Do for each value of M in MVAL.
*
DO I = 1, NM
M = MVAL( I )
*
* Do for each value of N in NVAL.
*
DO J = 1, NN
N = NVAL( J )
*
* Only for M >= N
*
IF ( MIN( M, N ).GT.0 .AND. M.GE.N ) THEN
*
* Do for each possible value of MB1
*
DO IMB1 = 1, NNB
MB1 = NBVAL( IMB1 )
*
* Only for MB1 > N
*
IF ( MB1.GT.N ) THEN
*
* Do for each possible value of NB1
*
DO INB1 = 1, NNB
NB1 = NBVAL( INB1 )
*
* Do for each possible value of NB2
*
DO INB2 = 1, NNB
NB2 = NBVAL( INB2 )
*
IF( NB1.GT.0 .AND. NB2.GT.0 ) THEN
*
* Test ZUNHR_COL
*
CALL ZUNHR_COL02( M, N, MB1, NB1,
$ NB2, RESULT )
*
* Print information about the tests that did
* not pass the threshold.
*
DO T = 1, NTESTS
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 ) M, N, MB1,
$ NB1, NB2, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + NTESTS
END IF
END DO
END DO
END IF
END DO
END IF
END DO
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 'M=', I5, ', N=', I5, ', MB1=', I5,
$ ', NB1=', I5, ', NB2=', I5,' test(', I2, ')=', G12.5 )
9999 FORMAT( 'ZUNGTSQR and ZUNHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
9998 FORMAT( 'ZUNGTSQR_ROW and ZUNHR_COL: M=', I5, ', N=', I5,
$ ', MB1=', I5, ', NB1=', I5, ', NB2=', I5,
$ ' test(', I2, ')=', G12.5 )
RETURN
*
* End of ZCHKUNHR_COL

View File

@ -21,8 +21,8 @@
*>
*> \verbatim
*>
*> ZUNHR_COL01 tests ZUNHR_COL using ZLATSQR, ZGEMQRT and ZUNGTSQR.
*> Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part ZGEMQR), ZUNGTSQR
*> ZUNHR_COL01 tests ZUNGTSQR and ZUNHR_COL using ZLATSQR, ZGEMQRT.
*> Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part of ZGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
@ -62,14 +62,46 @@
*> \verbatim
*> RESULT is DOUBLE PRECISION array, dimension (6)
*> Results of each of the six tests below.
*> ( C is a M-by-N random matrix, D is a N-by-M random matrix )
*>
*> RESULT(1) = | A - Q * R | / (eps * m * |A|)
*> RESULT(2) = | I - (Q**H) * Q | / (eps * m )
*> RESULT(3) = | Q * C - Q * C | / (eps * m * |C|)
*> RESULT(4) = | (Q**H) * C - (Q**H) * C | / (eps * m * |C|)
*> RESULT(5) = | (D * Q) - D * Q | / (eps * m * |D|)
*> RESULT(6) = | D * (Q**H) - D * (Q**H) | / (eps * m * |D|)
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m unitary Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in ZGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m unitary matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using ZGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using ZGEMM.
*> \endverbatim
*
* Authors:
@ -80,18 +112,15 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup complex16_lin
*
* =====================================================================
SUBROUTINE ZUNHR_COL01( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.9.0) --
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
@ -102,7 +131,7 @@
*
* ..
* .. Local allocatable arrays
COMPLEX*16, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
COMPLEX*16 , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
@ -218,7 +247,7 @@
* Copy the factor R into the array R.
*
SRNAMT = 'ZLACPY'
CALL ZLACPY( 'U', M, N, AF, M, R, M )
CALL ZLACPY( 'U', N, N, AF, M, R, M )
*
* Reconstruct the orthogonal matrix Q.
*
@ -240,7 +269,7 @@
* matrix S.
*
SRNAMT = 'ZLACPY'
CALL ZLACPY( 'U', M, N, R, M, AF, M )
CALL ZLACPY( 'U', N, N, R, M, AF, M )
*
DO I = 1, N
IF( DIAG( I ).EQ.-CONE ) THEN

View File

@ -0,0 +1,381 @@
*> \brief \b ZUNHR_COL02
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE ZUNHR_COL02( M, N, MB1, NB1, NB2, RESULT )
*
* .. Scalar Arguments ..
* INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
* DOUBLE PRECISION RESULT(6)
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZUNHR_COL02 tests ZUNGTSQR_ROW and ZUNHR_COL inside ZGETSQRHRT
*> (which calls ZLATSQR, ZUNGTSQR_ROW and ZUNHR_COL) using ZGEMQRT.
*> Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part of ZGEMQR)
*> have to be tested before this test.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> Number of rows in test matrix.
*> \endverbatim
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> Number of columns in test matrix.
*> \endverbatim
*> \param[in] MB1
*> \verbatim
*> MB1 is INTEGER
*> Number of row in row block in an input test matrix.
*> \endverbatim
*>
*> \param[in] NB1
*> \verbatim
*> NB1 is INTEGER
*> Number of columns in column block an input test matrix.
*> \endverbatim
*>
*> \param[in] NB2
*> \verbatim
*> NB2 is INTEGER
*> Number of columns in column block in an output test matrix.
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is DOUBLE PRECISION array, dimension (6)
*> Results of each of the six tests below.
*>
*> A is a m-by-n test input matrix to be factored.
*> so that A = Q_gr * ( R )
*> ( 0 ),
*>
*> Q_qr is an implicit m-by-m unitary Q matrix, the result
*> of factorization in blocked WY-representation,
*> stored in ZGEQRT output format.
*>
*> R is a n-by-n upper-triangular matrix,
*>
*> 0 is a (m-n)-by-n zero matrix,
*>
*> Q is an explicit m-by-m unitary matrix Q = Q_gr * I
*>
*> C is an m-by-n random matrix,
*>
*> D is an n-by-m random matrix.
*>
*> The six tests are:
*>
*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
*>
*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
*>
*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
*>
*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
*>
*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
*>
*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
*>
*> where:
*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
*> computed using ZGEMQRT,
*>
*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
*> computed using ZGEMM.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex16_lin
*
* =====================================================================
SUBROUTINE ZUNHR_COL02( M, N, MB1, NB1, NB2, RESULT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER M, N, MB1, NB1, NB2
* .. Return values ..
DOUBLE PRECISION RESULT(6)
*
* =====================================================================
*
* ..
* .. Local allocatable arrays
COMPLEX*16 , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
$ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
$ C(:,:), CF(:,:), D(:,:), DF(:,:)
DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
COMPLEX*16 CONE, CZERO
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
$ CZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
LOGICAL TESTZEROS
INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 )
COMPLEX*16 WORKQUERY( 1 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
EXTERNAL DLAMCH, ZLANGE, ZLANSY
* ..
* .. External Subroutines ..
EXTERNAL ZLACPY, ZLARNV, ZLASET, ZGETSQRHRT,
$ ZSCAL, ZGEMM, ZGEMQRT, ZHERK
* ..
* .. Intrinsic Functions ..
INTRINSIC CEILING, DBLE, MAX, MIN
* ..
* .. Scalars in Common ..
CHARACTER(LEN=32) SRNAMT
* ..
* .. Common blocks ..
COMMON / SRMNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEED / 1988, 1989, 1990, 1991 /
*
* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
*
TESTZEROS = .FALSE.
*
EPS = DLAMCH( 'Epsilon' )
K = MIN( M, N )
L = MAX( M, N, 1)
*
* Dynamically allocate local arrays
*
ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
$ C(M,N), CF(M,N),
$ D(N,M), DF(N,M) )
*
* Put random numbers into A and copy to AF
*
DO J = 1, N
CALL ZLARNV( 2, ISEED, M, A( 1, J ) )
END DO
IF( TESTZEROS ) THEN
IF( M.GE.4 ) THEN
DO J = 1, N
CALL ZLARNV( 2, ISEED, M/2, A( M/4, J ) )
END DO
END IF
END IF
CALL ZLACPY( 'Full', M, N, A, M, AF, M )
*
* Number of row blocks in ZLATSQR
*
NRB = MAX( 1, CEILING( DBLE( M - N ) / DBLE( MB1 - N ) ) )
*
ALLOCATE ( T1( NB1, N * NRB ) )
ALLOCATE ( T2( NB2, N ) )
ALLOCATE ( DIAG( N ) )
*
* Begin determine LWORK for the array WORK and allocate memory.
*
* ZGEMQRT requires NB2 to be bounded by N.
*
NB2_UB = MIN( NB2, N)
*
*
CALL ZGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORKQUERY, -1, INFO )
*
LWORK = INT( WORKQUERY( 1 ) )
*
* In ZGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
* or M*NB2_UB if SIDE = 'R'.
*
LWORK = MAX( LWORK, NB2_UB * N, NB2_UB * M )
*
ALLOCATE ( WORK( LWORK ) )
*
* End allocate memory for WORK.
*
*
* Begin Householder reconstruction routines
*
* Factor the matrix A in the array AF.
*
SRNAMT = 'ZGETSQRHRT'
CALL ZGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2,
$ WORK, LWORK, INFO )
*
* End Householder reconstruction routines.
*
*
* Generate the m-by-m matrix Q
*
CALL ZLASET( 'Full', M, M, CZERO, CONE, Q, M )
*
SRNAMT = 'ZGEMQRT'
CALL ZGEMQRT( 'L', 'N', M, M, K, NB2_UB, AF, M, T2, NB2, Q, M,
$ WORK, INFO )
*
* Copy R
*
CALL ZLASET( 'Full', M, N, CZERO, CZERO, R, M )
*
CALL ZLACPY( 'Upper', M, N, AF, M, R, M )
*
* TEST 1
* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
*
CALL ZGEMM( 'C', 'N', M, N, M, -CONE, Q, M, A, M, CONE, R, M )
*
ANORM = ZLANGE( '1', M, N, A, M, RWORK )
RESID = ZLANGE( '1', M, N, R, M, RWORK )
IF( ANORM.GT.ZERO ) THEN
RESULT( 1 ) = RESID / ( EPS * MAX( 1, M ) * ANORM )
ELSE
RESULT( 1 ) = ZERO
END IF
*
* TEST 2
* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
*
CALL ZLASET( 'Full', M, M, CZERO, CONE, R, M )
CALL ZHERK( 'U', 'C', M, M, -CONE, Q, M, CONE, R, M )
RESID = ZLANSY( '1', 'Upper', M, R, M, RWORK )
RESULT( 2 ) = RESID / ( EPS * MAX( 1, M ) )
*
* Generate random m-by-n matrix C
*
DO J = 1, N
CALL ZLARNV( 2, ISEED, M, C( 1, J ) )
END DO
CNORM = ZLANGE( '1', M, N, C, M, RWORK )
CALL ZLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as Q*C = CF
*
SRNAMT = 'ZGEMQRT'
CALL ZGEMQRT( 'L', 'N', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 3
* Compute |CF - Q*C| / ( eps * m * |C| )
*
CALL ZGEMM( 'N', 'N', M, N, M, -CONE, Q, M, C, M, CONE, CF, M )
RESID = ZLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 3 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 3 ) = ZERO
END IF
*
* Copy C into CF again
*
CALL ZLACPY( 'Full', M, N, C, M, CF, M )
*
* Apply Q to C as (Q**T)*C = CF
*
SRNAMT = 'ZGEMQRT'
CALL ZGEMQRT( 'L', 'C', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M,
$ WORK, INFO )
*
* TEST 4
* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
*
CALL ZGEMM( 'C', 'N', M, N, M, -CONE, Q, M, C, M, CONE, CF, M )
RESID = ZLANGE( '1', M, N, CF, M, RWORK )
IF( CNORM.GT.ZERO ) THEN
RESULT( 4 ) = RESID / ( EPS * MAX( 1, M ) * CNORM )
ELSE
RESULT( 4 ) = ZERO
END IF
*
* Generate random n-by-m matrix D and a copy DF
*
DO J = 1, M
CALL ZLARNV( 2, ISEED, N, D( 1, J ) )
END DO
DNORM = ZLANGE( '1', N, M, D, N, RWORK )
CALL ZLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*Q = DF
*
SRNAMT = 'ZGEMQRT'
CALL ZGEMQRT( 'R', 'N', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 5
* Compute |DF - D*Q| / ( eps * m * |D| )
*
CALL ZGEMM( 'N', 'N', N, M, M, -CONE, D, N, Q, M, CONE, DF, N )
RESID = ZLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 5 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 5 ) = ZERO
END IF
*
* Copy D into DF again
*
CALL ZLACPY( 'Full', N, M, D, N, DF, N )
*
* Apply Q to D as D*QT = DF
*
SRNAMT = 'ZGEMQRT'
CALL ZGEMQRT( 'R', 'C', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N,
$ WORK, INFO )
*
* TEST 6
* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
*
CALL ZGEMM( 'N', 'C', N, M, M, -CONE, D, N, Q, M, CONE, DF, N )
RESID = ZLANGE( '1', N, M, DF, N, RWORK )
IF( DNORM.GT.ZERO ) THEN
RESULT( 6 ) = RESID / ( EPS * MAX( 1, M ) * DNORM )
ELSE
RESULT( 6 ) = ZERO
END IF
*
* Deallocate all arrays
*
DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T1, T2, DIAG,
$ C, D, CF, DF )
*
RETURN
*
* End of ZUNHR_COL02
*
END