Move LAPACK trunk into position.

1 files changed, 177 insertions, 0 deletions

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