From baba851215b44ac3b60b9248eb02bcce7eb76247 Mon Sep 17 00:00:00 2001 From: jason Date: Tue, 28 Oct 2008 01:38:50 +0000 Subject: Move LAPACK trunk into position. --- TESTING/LIN/zrqt01.f | 177 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 TESTING/LIN/zrqt01.f (limited to 'TESTING/LIN/zrqt01.f') diff --git a/TESTING/LIN/zrqt01.f b/TESTING/LIN/zrqt01.f new file mode 100644 index 00000000..2a3a588e --- /dev/null +++ b/TESTING/LIN/zrqt01.f @@ -0,0 +1,177 @@ + SUBROUTINE ZRQT01( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK, + $ RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, LWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION RESULT( * ), RWORK( * ) + COMPLEX*16 A( LDA, * ), AF( LDA, * ), Q( LDA, * ), + $ R( LDA, * ), TAU( * ), WORK( LWORK ) +* .. +* +* Purpose +* ======= +* +* ZRQT01 tests ZGERQF, which computes the RQ factorization of an m-by-n +* matrix A, and partially tests ZUNGRQ which forms the n-by-n +* orthogonal matrix Q. +* +* ZRQT01 compares R with A*Q', and checks that Q is orthogonal. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input) COMPLEX*16 array, dimension (LDA,N) +* The m-by-n matrix A. +* +* AF (output) COMPLEX*16 array, dimension (LDA,N) +* Details of the RQ factorization of A, as returned by ZGERQF. +* See ZGERQF for further details. +* +* Q (output) COMPLEX*16 array, dimension (LDA,N) +* The n-by-n orthogonal matrix Q. +* +* R (workspace) COMPLEX*16 array, dimension (LDA,max(M,N)) +* +* LDA (input) INTEGER +* The leading dimension of the arrays A, AF, Q and L. +* LDA >= max(M,N). +* +* TAU (output) COMPLEX*16 array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors, as returned +* by ZGERQF. +* +* WORK (workspace) COMPLEX*16 array, dimension (LWORK) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (max(M,N)) +* +* RESULT (output) DOUBLE PRECISION array, dimension (2) +* The test ratios: +* RESULT(1) = norm( R - A*Q' ) / ( N * norm(A) * EPS ) +* RESULT(2) = norm( I - Q*Q' ) / ( N * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX*16 ROGUE + PARAMETER ( ROGUE = ( -1.0D+10, -1.0D+10 ) ) +* .. +* .. Local Scalars .. + INTEGER INFO, MINMN + DOUBLE PRECISION ANORM, EPS, RESID +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY + EXTERNAL DLAMCH, ZLANGE, ZLANSY +* .. +* .. External Subroutines .. + EXTERNAL ZGEMM, ZGERQF, ZHERK, ZLACPY, ZLASET, ZUNGRQ +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, DCMPLX, MAX, MIN +* .. +* .. Scalars in Common .. + CHARACTER(32) SRNAMT +* .. +* .. Common blocks .. + COMMON / SRNAMC / SRNAMT +* .. +* .. Executable Statements .. +* + MINMN = MIN( M, N ) + EPS = DLAMCH( 'Epsilon' ) +* +* Copy the matrix A to the array AF. +* + CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) +* +* Factorize the matrix A in the array AF. +* + SRNAMT = 'ZGERQF' + CALL ZGERQF( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy details of Q +* + CALL ZLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA ) + IF( M.LE.N ) THEN + IF( M.GT.0 .AND. M.LT.N ) + $ CALL ZLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA ) + IF( M.GT.1 ) + $ CALL ZLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA, + $ Q( N-M+2, N-M+1 ), LDA ) + ELSE + IF( N.GT.1 ) + $ CALL ZLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA, + $ Q( 2, 1 ), LDA ) + END IF +* +* Generate the n-by-n matrix Q +* + SRNAMT = 'ZUNGRQ' + CALL ZUNGRQ( N, N, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy R +* + CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), R, + $ LDA ) + IF( M.LE.N ) THEN + IF( M.GT.0 ) + $ CALL ZLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, + $ R( 1, N-M+1 ), LDA ) + ELSE + IF( M.GT.N .AND. N.GT.0 ) + $ CALL ZLACPY( 'Full', M-N, N, AF, LDA, R, LDA ) + IF( N.GT.0 ) + $ CALL ZLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, + $ R( M-N+1, 1 ), LDA ) + END IF +* +* Compute R - A*Q' +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', M, N, N, + $ DCMPLX( -ONE ), A, LDA, Q, LDA, DCMPLX( ONE ), R, + $ LDA ) +* +* Compute norm( R - Q'*A ) / ( N * norm(A) * EPS ) . +* + ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + RESID = ZLANGE( '1', M, N, R, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, N ) ) ) / ANORM ) / EPS + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute I - Q*Q' +* + CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ), DCMPLX( ONE ), R, LDA ) + CALL ZHERK( 'Upper', 'No transpose', N, N, -ONE, Q, LDA, ONE, R, + $ LDA ) +* +* Compute norm( I - Q*Q' ) / ( N * EPS ) . +* + RESID = ZLANSY( '1', 'Upper', N, R, LDA, RWORK ) +* + RESULT( 2 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / EPS +* + RETURN +* +* End of ZRQT01 +* + END -- cgit v1.2.3