![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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