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
      SUBROUTINE DBDT03UPLONKDDEULDUSVTLDVTWORK,
     $                   RESID )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            KDLDULDVTN
      DOUBLE PRECISION   RESID
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D* ), E* ), S* ), ULDU* ),
     $                   VTLDVT* ), WORK* )
*     ..
*
*  Purpose
*  =======
*
*  DBDT03 reconstructs a bidiagonal matrix B from its SVD:
*     S = U' * B * V
*  where U and V are orthogonal matrices and S is diagonal.
*
*  The test ratio to test the singular value decomposition is
*     RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
*  where VT = V' and EPS is the machine precision.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          Specifies whether the matrix B is upper or lower bidiagonal.
*          = 'U':  Upper bidiagonal
*          = 'L':  Lower bidiagonal
*
*  N       (input) INTEGER
*          The order of the matrix B.
*
*  KD      (input) INTEGER
*          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
*          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
*          not referenced.  If KD is greater than 1, it is assumed to be
*          1, and if KD is less than 0, it is assumed to be 0.
*
*  D       (input) DOUBLE PRECISION array, dimension (N)
*          The n diagonal elements of the bidiagonal matrix B.
*
*  E       (input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) superdiagonal elements of the bidiagonal matrix B
*          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
*          UPLO = 'L'.
*
*  U       (input) DOUBLE PRECISION array, dimension (LDU,N)
*          The n by n orthogonal matrix U in the reduction B = U'*A*P.
*
*  LDU     (input) INTEGER
*          The leading dimension of the array U.  LDU >= max(1,N)
*
*  S       (input) DOUBLE PRECISION array, dimension (N)
*          The singular values from the SVD of B, sorted in decreasing
*          order.
*
*  VT      (input) DOUBLE PRECISION array, dimension (LDVT,N)
*          The n by n orthogonal matrix V' in the reduction
*          B = U * S * V'.
*
*  LDVT    (input) INTEGER
*          The leading dimension of the array VT.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
*
*  RESID   (output) DOUBLE PRECISION
*          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
*
* ======================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZEROONE
      PARAMETER          ( ZERO = 0.0D+0ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            IJ
      DOUBLE PRECISION   BNORMEPS
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IDAMAX
      DOUBLE PRECISION   DASUMDLAMCH
      EXTERNAL           LSAMEIDAMAXDASUMDLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEMV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABSDBLEMAXMIN
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      RESID = ZERO
      IFN.LE.0 )
     $   RETURN
*
*     Compute B - U * S * V' one column at a time.
*
      BNORM = ZERO
      IFKD.GE.1 ) THEN
*
*        B is bidiagonal.
*
         IFLSAMEUPLO'U' ) ) THEN
*
*           B is upper bidiagonal.
*
            DO 20 J = 1N
               DO 10 I = 1N
                  WORKN+I ) = SI )*VTIJ )
   10          CONTINUE
               CALL DGEMV'No transpose'NN-ONEULDU,
     $                     WORKN+1 ), 1ZEROWORK1 )
               WORKJ ) = WORKJ ) + DJ )
               IFJ.GT.1 ) THEN
                  WORKJ-1 ) = WORKJ-1 ) + EJ-1 )
                  BNORM = MAXBNORMABSDJ ) )+ABSEJ-1 ) ) )
               ELSE
                  BNORM = MAXBNORMABSDJ ) ) )
               END IF
               RESID = MAXRESIDDASUMNWORK1 ) )
   20       CONTINUE
         ELSE
*
*           B is lower bidiagonal.
*
            DO 40 J = 1N
               DO 30 I = 1N
                  WORKN+I ) = SI )*VTIJ )
   30          CONTINUE
               CALL DGEMV'No transpose'NN-ONEULDU,
     $                     WORKN+1 ), 1ZEROWORK1 )
               WORKJ ) = WORKJ ) + DJ )
               IFJ.LT.N ) THEN
                  WORKJ+1 ) = WORKJ+1 ) + EJ )
                  BNORM = MAXBNORMABSDJ ) )+ABSEJ ) ) )
               ELSE
                  BNORM = MAXBNORMABSDJ ) ) )
               END IF
               RESID = MAXRESIDDASUMNWORK1 ) )
   40       CONTINUE
         END IF
      ELSE
*
*        B is diagonal.
*
         DO 60 J = 1N
            DO 50 I = 1N
               WORKN+I ) = SI )*VTIJ )
   50       CONTINUE
            CALL DGEMV'No transpose'NN-ONEULDUWORKN+1 ),
     $                  1ZEROWORK1 )
            WORKJ ) = WORKJ ) + DJ )
            RESID = MAXRESIDDASUMNWORK1 ) )
   60    CONTINUE
         J = IDAMAXND1 )
         BNORM = ABSDJ ) )
      END IF
*
*     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
*
      EPS = DLAMCH'Precision' )
*
      IFBNORM.LE.ZERO ) THEN
         IFRESID.NE.ZERO )
     $      RESID = ONE / EPS
      ELSE
         IFBNORM.GE.RESID ) THEN
            RESID = ( RESID / BNORM ) / ( DBLEN )*EPS )
         ELSE
            IFBNORM.LT.ONE ) THEN
               RESID = ( MINRESIDDBLEN )*BNORM ) / BNORM ) /
     $                 ( DBLEN )*EPS )
            ELSE
               RESID = MINRESID / BNORMDBLEN ) ) /
     $                 ( DBLEN )*EPS )
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of DBDT03
*
      END