LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chst01.f
Go to the documentation of this file.
00001 *> \brief \b CHST01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
00012 *                          LWORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            IHI, ILO, LDA, LDH, LDQ, LWORK, N
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       REAL               RESULT( 2 ), RWORK( * )
00019 *       COMPLEX            A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
00020 *      $                   WORK( LWORK )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> CHST01 tests the reduction of a general matrix A to upper Hessenberg
00030 *> form:  A = Q*H*Q'.  Two test ratios are computed;
00031 *>
00032 *> RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00033 *> RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
00034 *>
00035 *> The matrix Q is assumed to be given explicitly as it would be
00036 *> following CGEHRD + CUNGHR.
00037 *>
00038 *> In this version, ILO and IHI are not used, but they could be used
00039 *> to save some work if this is desired.
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] N
00046 *> \verbatim
00047 *>          N is INTEGER
00048 *>          The order of the matrix A.  N >= 0.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] ILO
00052 *> \verbatim
00053 *>          ILO is INTEGER
00054 *> \endverbatim
00055 *>
00056 *> \param[in] IHI
00057 *> \verbatim
00058 *>          IHI is INTEGER
00059 *>
00060 *>          A is assumed to be upper triangular in rows and columns
00061 *>          1:ILO-1 and IHI+1:N, so Q differs from the identity only in
00062 *>          rows and columns ILO+1:IHI.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] A
00066 *> \verbatim
00067 *>          A is COMPLEX array, dimension (LDA,N)
00068 *>          The original n by n matrix A.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] LDA
00072 *> \verbatim
00073 *>          LDA is INTEGER
00074 *>          The leading dimension of the array A.  LDA >= max(1,N).
00075 *> \endverbatim
00076 *>
00077 *> \param[in] H
00078 *> \verbatim
00079 *>          H is COMPLEX array, dimension (LDH,N)
00080 *>          The upper Hessenberg matrix H from the reduction A = Q*H*Q'
00081 *>          as computed by CGEHRD.  H is assumed to be zero below the
00082 *>          first subdiagonal.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] LDH
00086 *> \verbatim
00087 *>          LDH is INTEGER
00088 *>          The leading dimension of the array H.  LDH >= max(1,N).
00089 *> \endverbatim
00090 *>
00091 *> \param[in] Q
00092 *> \verbatim
00093 *>          Q is COMPLEX array, dimension (LDQ,N)
00094 *>          The orthogonal matrix Q from the reduction A = Q*H*Q' as
00095 *>          computed by CGEHRD + CUNGHR.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] LDQ
00099 *> \verbatim
00100 *>          LDQ is INTEGER
00101 *>          The leading dimension of the array Q.  LDQ >= max(1,N).
00102 *> \endverbatim
00103 *>
00104 *> \param[out] WORK
00105 *> \verbatim
00106 *>          WORK is COMPLEX array, dimension (LWORK)
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LWORK
00110 *> \verbatim
00111 *>          LWORK is INTEGER
00112 *>          The length of the array WORK.  LWORK >= 2*N*N.
00113 *> \endverbatim
00114 *>
00115 *> \param[out] RWORK
00116 *> \verbatim
00117 *>          RWORK is REAL array, dimension (N)
00118 *> \endverbatim
00119 *>
00120 *> \param[out] RESULT
00121 *> \verbatim
00122 *>          RESULT is REAL array, dimension (2)
00123 *>          RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00124 *>          RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
00125 *> \endverbatim
00126 *
00127 *  Authors:
00128 *  ========
00129 *
00130 *> \author Univ. of Tennessee 
00131 *> \author Univ. of California Berkeley 
00132 *> \author Univ. of Colorado Denver 
00133 *> \author NAG Ltd. 
00134 *
00135 *> \date November 2011
00136 *
00137 *> \ingroup complex_eig
00138 *
00139 *  =====================================================================
00140       SUBROUTINE CHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
00141      $                   LWORK, RWORK, RESULT )
00142 *
00143 *  -- LAPACK test routine (version 3.4.0) --
00144 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00145 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00146 *     November 2011
00147 *
00148 *     .. Scalar Arguments ..
00149       INTEGER            IHI, ILO, LDA, LDH, LDQ, LWORK, N
00150 *     ..
00151 *     .. Array Arguments ..
00152       REAL               RESULT( 2 ), RWORK( * )
00153       COMPLEX            A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
00154      $                   WORK( LWORK )
00155 *     ..
00156 *
00157 *  =====================================================================
00158 *
00159 *     .. Parameters ..
00160       REAL               ONE, ZERO
00161       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00162 *     ..
00163 *     .. Local Scalars ..
00164       INTEGER            LDWORK
00165       REAL               ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
00166 *     ..
00167 *     .. External Functions ..
00168       REAL               CLANGE, SLAMCH
00169       EXTERNAL           CLANGE, SLAMCH
00170 *     ..
00171 *     .. External Subroutines ..
00172       EXTERNAL           CGEMM, CLACPY, CUNT01, SLABAD
00173 *     ..
00174 *     .. Intrinsic Functions ..
00175       INTRINSIC          CMPLX, MAX, MIN
00176 *     ..
00177 *     .. Executable Statements ..
00178 *
00179 *     Quick return if possible
00180 *
00181       IF( N.LE.0 ) THEN
00182          RESULT( 1 ) = ZERO
00183          RESULT( 2 ) = ZERO
00184          RETURN
00185       END IF
00186 *
00187       UNFL = SLAMCH( 'Safe minimum' )
00188       EPS = SLAMCH( 'Precision' )
00189       OVFL = ONE / UNFL
00190       CALL SLABAD( UNFL, OVFL )
00191       SMLNUM = UNFL*N / EPS
00192 *
00193 *     Test 1:  Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
00194 *
00195 *     Copy A to WORK
00196 *
00197       LDWORK = MAX( 1, N )
00198       CALL CLACPY( ' ', N, N, A, LDA, WORK, LDWORK )
00199 *
00200 *     Compute Q*H
00201 *
00202       CALL CGEMM( 'No transpose', 'No transpose', N, N, N, CMPLX( ONE ),
00203      $            Q, LDQ, H, LDH, CMPLX( ZERO ), WORK( LDWORK*N+1 ),
00204      $            LDWORK )
00205 *
00206 *     Compute A - Q*H*Q'
00207 *
00208       CALL CGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
00209      $            CMPLX( -ONE ), WORK( LDWORK*N+1 ), LDWORK, Q, LDQ,
00210      $            CMPLX( ONE ), WORK, LDWORK )
00211 *
00212       ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
00213       WNORM = CLANGE( '1', N, N, WORK, LDWORK, RWORK )
00214 *
00215 *     Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
00216 *
00217       RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N
00218 *
00219 *     Test 2:  Compute norm( I - Q'*Q ) / ( N * EPS )
00220 *
00221       CALL CUNT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RWORK,
00222      $             RESULT( 2 ) )
00223 *
00224       RETURN
00225 *
00226 *     End of CHST01
00227 *
00228       END
 All Files Functions