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