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
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
     157
     158
     159
     160
     161
     162
     163
     164
     165
     166
     167
     168
     169
     170
     171
     172
     173
     174
     175
     176
     177
     178
     179
     180
     181
     182
     183
     184
     185
     186
     187
     188
     189
     190
     191
     192
     193
     194
     195
     196
     197
     198
     199
     200
     201
     202
     203
     204
     205
     206
     207
     208
     209
      SUBROUTINE DSBT21UPLONKAKSALDADEULDUWORK,
     $                   RESULT )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            KAKSLDALDUN
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   ALDA* ), D* ), E* ), RESULT2 ),
     $                   ULDU* ), WORK* )
*     ..
*
*  Purpose
*  =======
*
*  DSBT21  generally checks a decomposition of the form
*
*          A = U S U'
*
*  where ' means transpose, A is symmetric banded, U is
*  orthogonal, and S is diagonal (if KS=0) or symmetric
*  tridiagonal (if KS=1).
*
*  Specifically:
*
*          RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
*          RESULT(2) = | I - UU' | / ( n ulp )
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER
*          If UPLO='U', the upper triangle of A and V will be used and
*          the (strictly) lower triangle will not be referenced.
*          If UPLO='L', the lower triangle of A and V will be used and
*          the (strictly) upper triangle will not be referenced.
*
*  N       (input) INTEGER
*          The size of the matrix.  If it is zero, DSBT21 does nothing.
*          It must be at least zero.
*
*  KA      (input) INTEGER
*          The bandwidth of the matrix A.  It must be at least zero.  If
*          it is larger than N-1, then max( 0, N-1 ) will be used.
*
*  KS      (input) INTEGER
*          The bandwidth of the matrix S.  It may only be zero or one.
*          If zero, then S is diagonal, and E is not referenced.  If
*          one, then S is symmetric tri-diagonal.
*
*  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
*          The original (unfactored) matrix.  It is assumed to be
*          symmetric, and only the upper (UPLO='U') or only the lower
*          (UPLO='L') will be referenced.
*
*  LDA     (input) INTEGER
*          The leading dimension of A.  It must be at least 1
*          and at least min( KA, N-1 ).
*
*  D       (input) DOUBLE PRECISION array, dimension (N)
*          The diagonal of the (symmetric tri-) diagonal matrix S.
*
*  E       (input) DOUBLE PRECISION array, dimension (N-1)
*          The off-diagonal of the (symmetric tri-) diagonal matrix S.
*          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
*          (3,2) element, etc.
*          Not referenced if KS=0.
*
*  U       (input) DOUBLE PRECISION array, dimension (LDU, N)
*          The orthogonal matrix in the decomposition, expressed as a
*          dense matrix (i.e., not as a product of Householder
*          transformations, Givens transformations, etc.)
*
*  LDU     (input) INTEGER
*          The leading dimension of U.  LDU must be at least N and
*          at least 1.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (N**2+N)
*
*  RESULT  (output) DOUBLE PRECISION array, dimension (2)
*          The values computed by the two tests described above.  The
*          values are currently limited to 1/ulp, to avoid overflow.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZEROONE
      PARAMETER          ( ZERO = 0.0D0ONE = 1.0D0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LOWER
      CHARACTER          CUPLO
      INTEGER            IKAJJCJRLW
      DOUBLE PRECISION   ANORMULPUNFLWNORM
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCHDLANGEDLANSBDLANSP
      EXTERNAL           LSAMEDLAMCHDLANGEDLANSBDLANSP
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEMMDSPRDSPR2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DBLEMAXMIN
*     ..
*     .. Executable Statements ..
*
*     Constants
*
      RESULT1 ) = ZERO
      RESULT2 ) = ZERO
      IFN.LE.0 )
     $   RETURN
*
      IKA = MAX0MINN-1KA ) )
      LW = ( N*N+1 ) ) / 2
*
      IFLSAMEUPLO'U' ) ) THEN
         LOWER = .FALSE.
         CUPLO = 'U'
      ELSE
         LOWER = .TRUE.
         CUPLO = 'L'
      END IF
*
      UNFL = DLAMCH'Safe minimum' )
      ULP = DLAMCH'Epsilon' )*DLAMCH'Base' )
*
*     Some Error Checks
*
*     Do Test 1
*
*     Norm of A:
*
      ANORM = MAXDLANSB'1'CUPLONIKAALDAWORK ), UNFL )
*
*     Compute error matrix:    Error = A - U S U'
*
*     Copy A from SB to SP storage format.
*
      J = 0
      DO 50 JC = 1N
         IFLOWER ) THEN
            DO 10 JR = 1MINIKA+1N+1-JC )
               J = J + 1
               WORKJ ) = AJRJC )
   10       CONTINUE
            DO 20 JR = IKA + 2N + 1 - JC
               J = J + 1
               WORKJ ) = ZERO
   20       CONTINUE
         ELSE
            DO 30 JR = IKA + 2JC
               J = J + 1
               WORKJ ) = ZERO
   30       CONTINUE
            DO 40 JR = MINIKAJC-1 ), 0-1
               J = J + 1
               WORKJ ) = AIKA+1-JRJC )
   40       CONTINUE
         END IF
   50 CONTINUE
*
      DO 60 J = 1N
         CALL DSPRCUPLON-DJ ), U1J ), 1WORK )
   60 CONTINUE
*
      IFN.GT.1 .AND. KS.EQ.1 ) THEN
         DO 70 J = 1N - 1
            CALL DSPR2CUPLON-EJ ), U1J ), 1U1J+1 ), 1,
     $                  WORK )
   70    CONTINUE
      END IF
      WNORM = DLANSP'1'CUPLONWORKWORKLW+1 ) )
*
      IFANORM.GT.WNORM ) THEN
         RESULT1 ) = ( WNORM / ANORM ) / ( N*ULP )
      ELSE
         IFANORM.LT.ONE ) THEN
            RESULT1 ) = ( MINWNORMN*ANORM ) / ANORM ) / ( N*ULP )
         ELSE
            RESULT1 ) = MINWNORM / ANORMDBLEN ) ) / ( N*ULP )
         END IF
      END IF
*
*     Do Test 2
*
*     Compute  UU' - I
*
      CALL DGEMM'N''C'NNNONEULDUULDUZEROWORK,
     $            N )
*
      DO 80 J = 1N
         WORK( ( N+1 )*J-1 )+1 ) = WORK( ( N+1 )*J-1 )+1 ) - ONE
   80 CONTINUE
*
      RESULT2 ) = MINDLANGE'1'NNWORKNWORKN**2+1 ) ),
     $              DBLEN ) ) / ( N*ULP )
*
      RETURN
*
*     End of DSBT21
*
      END