Browse Files
       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
      SUBROUTINE SHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
     $                   LWORK, RESULT )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            IHI, ILO, LDA, LDH, LDQ, LWORK, N
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
     $                   RESULT( 2 ), WORK( LWORK )
*     ..
*
*  Purpose
*  =======
*
*  SHST01 tests the reduction of a general matrix A to upper Hessenberg
*  form:  A = Q*H*Q'.  Two test ratios are computed;
*
*  RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*  RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
*
*  The matrix Q is assumed to be given explicitly as it would be
*  following SGEHRD + SORGHR.
*
*  In this version, ILO and IHI are not used and are assumed to be 1 and
*  N, respectively.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  ILO     (input) INTEGER
*  IHI     (input) INTEGER
*          A is assumed to be upper triangular in rows and columns
*          1:ILO-1 and IHI+1:N, so Q differs from the identity only in
*          rows and columns ILO+1:IHI.
*
*  A       (input) REAL array, dimension (LDA,N)
*          The original n by n matrix A.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  H       (input) REAL array, dimension (LDH,N)
*          The upper Hessenberg matrix H from the reduction A = Q*H*Q'
*          as computed by SGEHRD.  H is assumed to be zero below the
*          first subdiagonal.
*
*  LDH     (input) INTEGER
*          The leading dimension of the array H.  LDH >= max(1,N).
*
*  Q       (input) REAL array, dimension (LDQ,N)
*          The orthogonal matrix Q from the reduction A = Q*H*Q' as
*          computed by SGEHRD + SORGHR.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q.  LDQ >= max(1,N).
*
*  WORK    (workspace) REAL array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The length of the array WORK.  LWORK >= 2*N*N.
*
*  RESULT  (output) REAL array, dimension (2)
*          RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*          RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            LDWORK
      REAL               ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLANGE
      EXTERNAL           SLAMCH, SLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMM, SLABAD, SLACPY, SORT01
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.LE.0 ) THEN
         RESULT( 1 ) = ZERO
         RESULT( 2 ) = ZERO
         RETURN
      END IF
*
      UNFL = SLAMCH( 'Safe minimum' )
      EPS = SLAMCH( 'Precision' )
      OVFL = ONE / UNFL
      CALL SLABAD( UNFL, OVFL )
      SMLNUM = UNFL*N / EPS
*
*     Test 1:  Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*
*     Copy A to WORK
*
      LDWORK = MAX( 1, N )
      CALL SLACPY( ' ', N, N, A, LDA, WORK, LDWORK )
*
*     Compute Q*H
*
      CALL SGEMM( 'No transpose', 'No transpose', N, N, N, ONE, Q, LDQ,
     $            H, LDH, ZERO, WORK( LDWORK*N+1 ), LDWORK )
*
*     Compute A - Q*H*Q'
*
      CALL SGEMM( 'No transpose', 'Transpose', N, N, N, -ONE,
     $            WORK( LDWORK*N+1 ), LDWORK, Q, LDQ, ONE, WORK,
     $            LDWORK )
*
      ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK( LDWORK*N+1 ) ),
     $        UNFL )
      WNORM = SLANGE( '1', N, N, WORK, LDWORK, WORK( LDWORK*N+1 ) )
*
*     Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
*
      RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N
*
*     Test 2:  Compute norm( I - Q'*Q ) / ( N * EPS )
*
      CALL SORT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RESULT( 2 ) )
*
      RETURN
*
*     End of SHST01
*
      END