TESTING/EIG/sbdt05.f

   1 *  =========== DOCUMENTATION ===========
   2 *
   3 * Online html documentation available at
   4 *            http://www.netlib.org/lapack/explore-html/
   5 *
   6 *  Definition:
   7 *  ===========
   8 *
   9 *       SUBROUTINE SBDT05( M, N, A, LDA, S, NS, U, LDU,
  10 *                          VT, LDVT, WORK, RESID )
  11 *
  12 *       .. Scalar Arguments ..
  13 *       INTEGER            LDA, LDU, LDVT, N, NS
  14 *       REAL               RESID
  15 *       ..
  16 *       .. Array Arguments ..
  17 *       REAL               D( * ), E( * ), S( * ), U( LDU, * ),
  18 *      $                   VT( LDVT, * ), WORK( * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *> SBDT05 reconstructs a bidiagonal matrix B from its (partial) SVD:
  28 *>    S = U' * B * V
  29 *> where U and V are orthogonal matrices and S is diagonal.
  30 *>
  31 *> The test ratio to test the singular value decomposition is
  32 *>    RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
  33 *> where VT = V' and EPS is the machine precision.
  34 *> \endverbatim
  35 *
  36 *  Arguments:
  37 *  ==========
  38 *
  39 *> \param[in] M
  40 *> \verbatim
  41 *>          M is INTEGER
  42 *>          The number of rows of the matrices A and U.
  43 *> \endverbatim
  44 *>
  45 *> \param[in] N
  46 *> \verbatim
  47 *>          N is INTEGER
  48 *>          The number of columns of the matrices A and VT.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] A
  52 *> \verbatim
  53 *>          A is REAL array, dimension (LDA,N)
  54 *>          The m by n matrix A.
  55 *>
  56 *> \param[in] LDA
  57 *> \verbatim
  58 *>          LDA is INTEGER
  59 *>          The leading dimension of the array A.  LDA >= max(1,M).
  60 *> \endverbatim
  61 *>
  62 *> \param[in] S
  63 *> \verbatim
  64 *>          S is REAL array, dimension (NS)
  65 *>          The singular values from the (partial) SVD of B, sorted in
  66 *>          decreasing order.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] NS
  70 *> \verbatim
  71 *>          NS is INTEGER
  72 *>          The number of singular values/vectors from the (partial)
  73 *>          SVD of B.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] U
  77 *> \verbatim
  78 *>          U is REAL array, dimension (LDU,NS)
  79 *>          The n by ns orthogonal matrix U in S = U' * B * V.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LDU
  83 *> \verbatim
  84 *>          LDU is INTEGER
  85 *>          The leading dimension of the array U.  LDU >= max(1,N)
  86 *> \endverbatim
  87 *>
  88 *> \param[in] VT
  89 *> \verbatim
  90 *>          VT is REAL array, dimension (LDVT,N)
  91 *>          The n by ns orthogonal matrix V in S = U' * B * V.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] LDVT
  95 *> \verbatim
  96 *>          LDVT is INTEGER
  97 *>          The leading dimension of the array VT.
  98 *> \endverbatim
  99 *>
 100 *> \param[out] WORK
 101 *> \verbatim
 102 *>          WORK is REAL array, dimension (M,N)
 103 *> \endverbatim
 104 *>
 105 *> \param[out] RESID
 106 *> \verbatim
 107 *>          RESID is REAL
 108 *>          The test ratio:  norm(S - U' * A * V) / ( n * norm(A) * EPS )
 109 *> \endverbatim
 110 *
 111 *  Authors:
 112 *  ========
 113 *
 114 *> \author Univ. of Tennessee
 115 *> \author Univ. of California Berkeley
 116 *> \author Univ. of Colorado Denver
 117 *> \author NAG Ltd.
 118 *
 119 *> \date November 2011
 120 *
 121 *> \ingroup double_eig
 122 *
 123 *  =====================================================================
 124       SUBROUTINE SBDT05( M, N, A, LDA, S, NS, U, LDU,
 125      $                    VT, LDVT, WORK, RESID )
 126 *
 127 *  -- LAPACK test routine (version 3.4.0) --
 128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 130 *     November 2011
 131 *
 132 *     .. Scalar Arguments ..
 133       CHARACTER          UPLO
 134       INTEGER            LDA, LDU, LDVT, M, N, NS
 135       REAL               RESID
 136 *     ..
 137 *     .. Array Arguments ..
 138       REAL               A( LDA, * ), S( * ), U( LDU, * ),
 139      $                   VT( LDVT, * ), WORK( * )
 140 *     ..
 141 *
 142 * ======================================================================
 143 *
 144 *     .. Parameters ..
 145       REAL               ZERO, ONE
 146       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 147 *     ..
 148 *     .. Local Scalars ..
 149       INTEGER            I, J
 150       REAL               ANORM, EPS
 151 *     ..
 152 *     .. External Functions ..
 153       LOGICAL            LSAME
 154       INTEGER            ISAMAX
 155       REAL               SASUM, SLAMCH, SLANGE
 156       EXTERNAL           LSAME, ISAMAX, SASUM, SLAMCH, SLANGE
 157 *     ..
 158 *     .. External Subroutines ..
 159       EXTERNAL           SGEMM
 160 *     ..
 161 *     .. Intrinsic Functions ..
 162       INTRINSIC          ABS, REAL, MAX, MIN
 163 *     ..
 164 *     .. Executable Statements ..
 165 *
 166 *     Quick return if possible.
 167 *
 168       RESID = ZERO
 169       IF( MIN( M, N ).LE.0 .OR. NS.LE.0 )
 170      $   RETURN
 171 *
 172       EPS = SLAMCH( 'Precision' )
 173       ANORM = SLANGE( 'M', M, N, A, LDA, WORK )
 174 *
 175 *     Compute U' * A * V.
 176 *
 177       CALL SGEMM( 'N', 'T', M, NS, N, ONE, A, LDA, VT,
 178      $            LDVT, ZERO, WORK( 1+NS*NS ), M )
 179       CALL SGEMM( 'T', 'N', NS, NS, M, -ONE, U, LDU, WORK( 1+NS*NS ),
 180      $            M, ZERO, WORK, NS )
 181 *
 182 *     norm(S - U' * B * V)
 183 *
 184       J = 0
 185       DO 10 I = 1, NS
 186          WORK( J+I ) =  WORK( J+I ) + S( I )
 187          RESID = MAX( RESID, SASUM( NS, WORK( J+1 ), 1 ) )
 188          J = J + NS
 189    10 CONTINUE
 190 *
 191       IF( ANORM.LE.ZERO ) THEN
 192          IF( RESID.NE.ZERO )
 193      $      RESID = ONE / EPS
 194       ELSE
 195          IF( ANORM.GE.RESID ) THEN
 196             RESID = ( RESID / ANORM ) / ( REAL( N )*EPS )
 197          ELSE
 198             IF( ANORM.LT.ONE ) THEN
 199                RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) /
 200      $                 ( REAL( N )*EPS )
 201             ELSE
 202                RESID = MIN( RESID / ANORM, REAL( N ) ) /
 203      $                 ( REAL( N )*EPS )
 204             END IF
 205          END IF
 206       END IF
 207 *
 208       RETURN
 209 *
 210 *     End of SBDT05
 211 *
 212       END