LAPACK 3.3.0

chst01.f

Go to the documentation of this file.
00001       SUBROUTINE CHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
00002      $                   LWORK, RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            IHI, ILO, LDA, LDH, LDQ, LWORK, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       REAL               RESULT( 2 ), RWORK( * )
00013       COMPLEX            A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
00014      $                   WORK( LWORK )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CHST01 tests the reduction of a general matrix A to upper Hessenberg
00021 *  form:  A = Q*H*Q'.  Two test ratios are computed;
00022 *
00023 *  RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00024 *  RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
00025 *
00026 *  The matrix Q is assumed to be given explicitly as it would be
00027 *  following CGEHRD + CUNGHR.
00028 *
00029 *  In this version, ILO and IHI are not used, but they could be used
00030 *  to save some work if this is desired.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  N       (input) INTEGER
00036 *          The order of the matrix A.  N >= 0.
00037 *
00038 *  ILO     (input) INTEGER
00039 *  IHI     (input) INTEGER
00040 *          A is assumed to be upper triangular in rows and columns
00041 *          1:ILO-1 and IHI+1:N, so Q differs from the identity only in
00042 *          rows and columns ILO+1:IHI.
00043 *
00044 *  A       (input) COMPLEX array, dimension (LDA,N)
00045 *          The original n by n matrix A.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A.  LDA >= max(1,N).
00049 *
00050 *  H       (input) COMPLEX array, dimension (LDH,N)
00051 *          The upper Hessenberg matrix H from the reduction A = Q*H*Q'
00052 *          as computed by CGEHRD.  H is assumed to be zero below the
00053 *          first subdiagonal.
00054 *
00055 *  LDH     (input) INTEGER
00056 *          The leading dimension of the array H.  LDH >= max(1,N).
00057 *
00058 *  Q       (input) COMPLEX array, dimension (LDQ,N)
00059 *          The orthogonal matrix Q from the reduction A = Q*H*Q' as
00060 *          computed by CGEHRD + CUNGHR.
00061 *
00062 *  LDQ     (input) INTEGER
00063 *          The leading dimension of the array Q.  LDQ >= max(1,N).
00064 *
00065 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
00066 *
00067 *  LWORK   (input) INTEGER
00068 *          The length of the array WORK.  LWORK >= 2*N*N.
00069 *
00070 *  RWORK   (workspace) REAL array, dimension (N)
00071 *
00072 *  RESULT  (output) REAL array, dimension (2)
00073 *          RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00074 *          RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
00075 *
00076 *  =====================================================================
00077 *
00078 *     .. Parameters ..
00079       REAL               ONE, ZERO
00080       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       INTEGER            LDWORK
00084       REAL               ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
00085 *     ..
00086 *     .. External Functions ..
00087       REAL               CLANGE, SLAMCH
00088       EXTERNAL           CLANGE, SLAMCH
00089 *     ..
00090 *     .. External Subroutines ..
00091       EXTERNAL           CGEMM, CLACPY, CUNT01, SLABAD
00092 *     ..
00093 *     .. Intrinsic Functions ..
00094       INTRINSIC          CMPLX, MAX, MIN
00095 *     ..
00096 *     .. Executable Statements ..
00097 *
00098 *     Quick return if possible
00099 *
00100       IF( N.LE.0 ) THEN
00101          RESULT( 1 ) = ZERO
00102          RESULT( 2 ) = ZERO
00103          RETURN
00104       END IF
00105 *
00106       UNFL = SLAMCH( 'Safe minimum' )
00107       EPS = SLAMCH( 'Precision' )
00108       OVFL = ONE / UNFL
00109       CALL SLABAD( UNFL, OVFL )
00110       SMLNUM = UNFL*N / EPS
00111 *
00112 *     Test 1:  Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00113 *
00114 *     Copy A to WORK
00115 *
00116       LDWORK = MAX( 1, N )
00117       CALL CLACPY( ' ', N, N, A, LDA, WORK, LDWORK )
00118 *
00119 *     Compute Q*H
00120 *
00121       CALL CGEMM( 'No transpose', 'No transpose', N, N, N, CMPLX( ONE ),
00122      $            Q, LDQ, H, LDH, CMPLX( ZERO ), WORK( LDWORK*N+1 ),
00123      $            LDWORK )
00124 *
00125 *     Compute A - Q*H*Q'
00126 *
00127       CALL CGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
00128      $            CMPLX( -ONE ), WORK( LDWORK*N+1 ), LDWORK, Q, LDQ,
00129      $            CMPLX( ONE ), WORK, LDWORK )
00130 *
00131       ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
00132       WNORM = CLANGE( '1', N, N, WORK, LDWORK, RWORK )
00133 *
00134 *     Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
00135 *
00136       RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N
00137 *
00138 *     Test 2:  Compute norm( I - Q'*Q ) / ( N * EPS )
00139 *
00140       CALL CUNT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RWORK,
00141      $             RESULT( 2 ) )
00142 *
00143       RETURN
00144 *
00145 *     End of CHST01
00146 *
00147       END
 All Files Functions