diff --git a/lapack-netlib/TESTING/LIN/CMakeLists.txt b/lapack-netlib/TESTING/LIN/CMakeLists.txt index 0d0bb5418..309ed7e77 100644 --- a/lapack-netlib/TESTING/LIN/CMakeLists.txt +++ b/lapack-netlib/TESTING/LIN/CMakeLists.txt @@ -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 diff --git a/lapack-netlib/TESTING/LIN/Makefile b/lapack-netlib/TESTING/LIN/Makefile index 6e790aa93..674265816 100644 --- a/lapack-netlib/TESTING/LIN/Makefile +++ b/lapack-netlib/TESTING/LIN/Makefile @@ -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 \ diff --git a/lapack-netlib/TESTING/LIN/cchkunhr_col.f b/lapack-netlib/TESTING/LIN/cchkunhr_col.f index 00077ddd9..0d6a9063d 100644 --- a/lapack-netlib/TESTING/LIN/cchkunhr_col.f +++ b/lapack-netlib/TESTING/LIN/cchkunhr_col.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/cunhr_col01.f b/lapack-netlib/TESTING/LIN/cunhr_col01.f index d760caba5..d77d60b1a 100644 --- a/lapack-netlib/TESTING/LIN/cunhr_col01.f +++ b/lapack-netlib/TESTING/LIN/cunhr_col01.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/cunhr_col02.f b/lapack-netlib/TESTING/LIN/cunhr_col02.f new file mode 100644 index 000000000..001f291da --- /dev/null +++ b/lapack-netlib/TESTING/LIN/cunhr_col02.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/dchkorhr_col.f b/lapack-netlib/TESTING/LIN/dchkorhr_col.f index 3b3e421eb..0e2d44d8d 100644 --- a/lapack-netlib/TESTING/LIN/dchkorhr_col.f +++ b/lapack-netlib/TESTING/LIN/dchkorhr_col.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/dorhr_col01.f b/lapack-netlib/TESTING/LIN/dorhr_col01.f index 3e48de37f..979255ca9 100644 --- a/lapack-netlib/TESTING/LIN/dorhr_col01.f +++ b/lapack-netlib/TESTING/LIN/dorhr_col01.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/dorhr_col02.f b/lapack-netlib/TESTING/LIN/dorhr_col02.f new file mode 100644 index 000000000..d4c438edb --- /dev/null +++ b/lapack-netlib/TESTING/LIN/dorhr_col02.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/schkorhr_col.f b/lapack-netlib/TESTING/LIN/schkorhr_col.f index cf6d2d323..f61b74902 100644 --- a/lapack-netlib/TESTING/LIN/schkorhr_col.f +++ b/lapack-netlib/TESTING/LIN/schkorhr_col.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/sorhr_col01.f b/lapack-netlib/TESTING/LIN/sorhr_col01.f index 02429041b..dcc2c1cae 100644 --- a/lapack-netlib/TESTING/LIN/sorhr_col01.f +++ b/lapack-netlib/TESTING/LIN/sorhr_col01.f @@ -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. diff --git a/lapack-netlib/TESTING/LIN/sorhr_col02.f b/lapack-netlib/TESTING/LIN/sorhr_col02.f new file mode 100644 index 000000000..1cbe40577 --- /dev/null +++ b/lapack-netlib/TESTING/LIN/sorhr_col02.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/zchkunhr_col.f b/lapack-netlib/TESTING/LIN/zchkunhr_col.f index ef8f8bcc4..395ea178a 100644 --- a/lapack-netlib/TESTING/LIN/zchkunhr_col.f +++ b/lapack-netlib/TESTING/LIN/zchkunhr_col.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/zunhr_col01.f b/lapack-netlib/TESTING/LIN/zunhr_col01.f index 9fb3bf352..b7590a8ea 100644 --- a/lapack-netlib/TESTING/LIN/zunhr_col01.f +++ b/lapack-netlib/TESTING/LIN/zunhr_col01.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/zunhr_col02.f b/lapack-netlib/TESTING/LIN/zunhr_col02.f new file mode 100644 index 000000000..c6e7f80cd --- /dev/null +++ b/lapack-netlib/TESTING/LIN/zunhr_col02.f @@ -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