Files
lapack/TESTING/LIN/cqrt01.f
Julie e18d437924 Making LAPACK's code eternal... no more version and date in source files.
GitHub is now enabling us to track accurately version and date.
No need for this anymore.
2021-03-25 10:16:58 -07:00

229 lines
6.0 KiB
FortranFixed

*> \brief \b CQRT01
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE CQRT01( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK,
* RWORK, RESULT )
*
* .. Scalar Arguments ..
* INTEGER LDA, LWORK, M, N
* ..
* .. Array Arguments ..
* REAL RESULT( * ), RWORK( * )
* COMPLEX A( LDA, * ), AF( LDA, * ), Q( LDA, * ),
* $ R( LDA, * ), TAU( * ), WORK( LWORK )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CQRT01 tests CGEQRF, which computes the QR factorization of an m-by-n
*> matrix A, and partially tests CUNGQR which forms the m-by-m
*> orthogonal matrix Q.
*>
*> CQRT01 compares R with Q'*A, and checks that Q is orthogonal.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of columns of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is COMPLEX array, dimension (LDA,N)
*> The m-by-n matrix A.
*> \endverbatim
*>
*> \param[out] AF
*> \verbatim
*> AF is COMPLEX array, dimension (LDA,N)
*> Details of the QR factorization of A, as returned by CGEQRF.
*> See CGEQRF for further details.
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*> Q is COMPLEX array, dimension (LDA,M)
*> The m-by-m orthogonal matrix Q.
*> \endverbatim
*>
*> \param[out] R
*> \verbatim
*> R is COMPLEX array, dimension (LDA,max(M,N))
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the arrays A, AF, Q and R.
*> LDA >= max(M,N).
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is COMPLEX array, dimension (min(M,N))
*> The scalar factors of the elementary reflectors, as returned
*> by CGEQRF.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*> RWORK is REAL array, dimension (M)
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is REAL array, dimension (2)
*> The test ratios:
*> RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS )
*> RESULT(2) = norm( I - Q'*Q ) / ( M * EPS )
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex_lin
*
* =====================================================================
SUBROUTINE CQRT01( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK,
$ RWORK, RESULT )
*
* -- 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 LDA, LWORK, M, N
* ..
* .. Array Arguments ..
REAL RESULT( * ), RWORK( * )
COMPLEX A( LDA, * ), AF( LDA, * ), Q( LDA, * ),
$ R( LDA, * ), TAU( * ), WORK( LWORK )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
COMPLEX ROGUE
PARAMETER ( ROGUE = ( -1.0E+10, -1.0E+10 ) )
* ..
* .. Local Scalars ..
INTEGER INFO, MINMN
REAL ANORM, EPS, RESID
* ..
* .. External Functions ..
REAL CLANGE, CLANSY, SLAMCH
EXTERNAL CLANGE, CLANSY, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CGEMM, CGEQRF, CHERK, CLACPY, CLASET, CUNGQR
* ..
* .. Intrinsic Functions ..
INTRINSIC CMPLX, MAX, MIN, REAL
* ..
* .. Scalars in Common ..
CHARACTER*32 SRNAMT
* ..
* .. Common blocks ..
COMMON / SRNAMC / SRNAMT
* ..
* .. Executable Statements ..
*
MINMN = MIN( M, N )
EPS = SLAMCH( 'Epsilon' )
*
* Copy the matrix A to the array AF.
*
CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
*
* Factorize the matrix A in the array AF.
*
SRNAMT = 'CGEQRF'
CALL CGEQRF( M, N, AF, LDA, TAU, WORK, LWORK, INFO )
*
* Copy details of Q
*
CALL CLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA )
CALL CLACPY( 'Lower', M-1, N, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA )
*
* Generate the m-by-m matrix Q
*
SRNAMT = 'CUNGQR'
CALL CUNGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO )
*
* Copy R
*
CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), R, LDA )
CALL CLACPY( 'Upper', M, N, AF, LDA, R, LDA )
*
* Compute R - Q'*A
*
CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M,
$ CMPLX( -ONE ), Q, LDA, A, LDA, CMPLX( ONE ), R, LDA )
*
* Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) .
*
ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
RESID = CLANGE( '1', M, N, R, LDA, RWORK )
IF( ANORM.GT.ZERO ) THEN
RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M ) ) ) / ANORM ) / EPS
ELSE
RESULT( 1 ) = ZERO
END IF
*
* Compute I - Q'*Q
*
CALL CLASET( 'Full', M, M, CMPLX( ZERO ), CMPLX( ONE ), R, LDA )
CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, Q, LDA,
$ ONE, R, LDA )
*
* Compute norm( I - Q'*Q ) / ( M * EPS ) .
*
RESID = CLANSY( '1', 'Upper', M, R, LDA, RWORK )
*
RESULT( 2 ) = ( RESID / REAL( MAX( 1, M ) ) ) / EPS
*
RETURN
*
* End of CQRT01
*
END