LAPACK 3.3.0

dhst01.f

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