LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dhst01.f
Go to the documentation of this file.
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
 All Files Functions