TESTING/EIG/sbdt01.f

   1 *> \brief \b SBDT01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE SBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
  12 *                          RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            KD, LDA, LDPT, LDQ, M, N
  16 *       REAL               RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
  20 *      $                   Q( LDQ, * ), WORK( * )
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *> SBDT01 reconstructs a general matrix A from its bidiagonal form
  30 *>    A = Q * B * P'
  31 *> where Q (m by min(m,n)) and P' (min(m,n) by n) are orthogonal
  32 *> matrices and B is bidiagonal.
  33 *>
  34 *> The test ratio to test the reduction is
  35 *>    RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
  36 *> where PT = P' and EPS is the machine precision.
  37 *> \endverbatim
  38 *
  39 *  Arguments:
  40 *  ==========
  41 *
  42 *> \param[in] M
  43 *> \verbatim
  44 *>          M is INTEGER
  45 *>          The number of rows of the matrices A and Q.
  46 *> \endverbatim
  47 *>
  48 *> \param[in] N
  49 *> \verbatim
  50 *>          N is INTEGER
  51 *>          The number of columns of the matrices A and P'.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] KD
  55 *> \verbatim
  56 *>          KD is INTEGER
  57 *>          If KD = 0, B is diagonal and the array E is not referenced.
  58 *>          If KD = 1, the reduction was performed by xGEBRD; B is upper
  59 *>          bidiagonal if M >= N, and lower bidiagonal if M < N.
  60 *>          If KD = -1, the reduction was performed by xGBBRD; B is
  61 *>          always upper bidiagonal.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] A
  65 *> \verbatim
  66 *>          A is REAL array, dimension (LDA,N)
  67 *>          The m by n matrix A.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] LDA
  71 *> \verbatim
  72 *>          LDA is INTEGER
  73 *>          The leading dimension of the array A.  LDA >= max(1,M).
  74 *> \endverbatim
  75 *>
  76 *> \param[in] Q
  77 *> \verbatim
  78 *>          Q is REAL array, dimension (LDQ,N)
  79 *>          The m by min(m,n) orthogonal matrix Q in the reduction
  80 *>          A = Q * B * P'.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] LDQ
  84 *> \verbatim
  85 *>          LDQ is INTEGER
  86 *>          The leading dimension of the array Q.  LDQ >= max(1,M).
  87 *> \endverbatim
  88 *>
  89 *> \param[in] D
  90 *> \verbatim
  91 *>          D is REAL array, dimension (min(M,N))
  92 *>          The diagonal elements of the bidiagonal matrix B.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] E
  96 *> \verbatim
  97 *>          E is REAL array, dimension (min(M,N)-1)
  98 *>          The superdiagonal elements of the bidiagonal matrix B if
  99 *>          m >= n, or the subdiagonal elements of B if m < n.
 100 *> \endverbatim
 101 *>
 102 *> \param[in] PT
 103 *> \verbatim
 104 *>          PT is REAL array, dimension (LDPT,N)
 105 *>          The min(m,n) by n orthogonal matrix P' in the reduction
 106 *>          A = Q * B * P'.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LDPT
 110 *> \verbatim
 111 *>          LDPT is INTEGER
 112 *>          The leading dimension of the array PT.
 113 *>          LDPT >= max(1,min(M,N)).
 114 *> \endverbatim
 115 *>
 116 *> \param[out] WORK
 117 *> \verbatim
 118 *>          WORK is REAL array, dimension (M+N)
 119 *> \endverbatim
 120 *>
 121 *> \param[out] RESID
 122 *> \verbatim
 123 *>          RESID is REAL
 124 *>          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS )
 125 *> \endverbatim
 126 *
 127 *  Authors:
 128 *  ========
 129 *
 130 *> \author Univ. of Tennessee
 131 *> \author Univ. of California Berkeley
 132 *> \author Univ. of Colorado Denver
 133 *> \author NAG Ltd.
 134 *
 135 *> \date November 2011
 136 *
 137 *> \ingroup single_eig
 138 *
 139 *  =====================================================================
 140       SUBROUTINE SBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
 141      $                   RESID )
 142 *
 143 *  -- LAPACK test routine (version 3.4.0) --
 144 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 145 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 146 *     November 2011
 147 *
 148 *     .. Scalar Arguments ..
 149       INTEGER            KD, LDA, LDPT, LDQ, M, N
 150       REAL               RESID
 151 *     ..
 152 *     .. Array Arguments ..
 153       REAL               A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
 154      $                   Q( LDQ, * ), WORK( * )
 155 *     ..
 156 *
 157 *  =====================================================================
 158 *
 159 *     .. Parameters ..
 160       REAL               ZERO, ONE
 161       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 162 *     ..
 163 *     .. Local Scalars ..
 164       INTEGER            I, J
 165       REAL               ANORM, EPS
 166 *     ..
 167 *     .. External Functions ..
 168       REAL               SASUM, SLAMCH, SLANGE
 169       EXTERNAL           SASUM, SLAMCH, SLANGE
 170 *     ..
 171 *     .. External Subroutines ..
 172       EXTERNAL           SCOPY, SGEMV
 173 *     ..
 174 *     .. Intrinsic Functions ..
 175       INTRINSIC          MAX, MIN, REAL
 176 *     ..
 177 *     .. Executable Statements ..
 178 *
 179 *     Quick return if possible
 180 *
 181       IF( M.LE.0 .OR. N.LE.0 ) THEN
 182          RESID = ZERO
 183          RETURN
 184       END IF
 185 *
 186 *     Compute A - Q * B * P' one column at a time.
 187 *
 188       RESID = ZERO
 189       IF( KD.NE.0 ) THEN
 190 *
 191 *        B is bidiagonal.
 192 *
 193          IF( KD.NE.0 .AND. M.GE.N ) THEN
 194 *
 195 *           B is upper bidiagonal and M >= N.
 196 *
 197             DO 20 J = 1, N
 198                CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
 199                DO 10 I = 1, N - 1
 200                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
 201    10          CONTINUE
 202                WORK( M+N ) = D( N )*PT( N, J )
 203                CALL SGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
 204      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
 205                RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 206    20       CONTINUE
 207          ELSE IF( KD.LT.0 ) THEN
 208 *
 209 *           B is upper bidiagonal and M < N.
 210 *
 211             DO 40 J = 1, N
 212                CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
 213                DO 30 I = 1, M - 1
 214                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
 215    30          CONTINUE
 216                WORK( M+M ) = D( M )*PT( M, J )
 217                CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
 218      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
 219                RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 220    40       CONTINUE
 221          ELSE
 222 *
 223 *           B is lower bidiagonal.
 224 *
 225             DO 60 J = 1, N
 226                CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
 227                WORK( M+1 ) = D( 1 )*PT( 1, J )
 228                DO 50 I = 2, M
 229                   WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
 230      $                          D( I )*PT( I, J )
 231    50          CONTINUE
 232                CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
 233      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
 234                RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 235    60       CONTINUE
 236          END IF
 237       ELSE
 238 *
 239 *        B is diagonal.
 240 *
 241          IF( M.GE.N ) THEN
 242             DO 80 J = 1, N
 243                CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
 244                DO 70 I = 1, N
 245                   WORK( M+I ) = D( I )*PT( I, J )
 246    70          CONTINUE
 247                CALL SGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
 248      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
 249                RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 250    80       CONTINUE
 251          ELSE
 252             DO 100 J = 1, N
 253                CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
 254                DO 90 I = 1, M
 255                   WORK( M+I ) = D( I )*PT( I, J )
 256    90          CONTINUE
 257                CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
 258      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
 259                RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 260   100       CONTINUE
 261          END IF
 262       END IF
 263 *
 264 *     Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
 265 *
 266       ANORM = SLANGE( '1', M, N, A, LDA, WORK )
 267       EPS = SLAMCH( 'Precision' )
 268 *
 269       IF( ANORM.LE.ZERO ) THEN
 270          IF( RESID.NE.ZERO )
 271      $      RESID = ONE / EPS
 272       ELSE
 273          IF( ANORM.GE.RESID ) THEN
 274             RESID = ( RESID / ANORM ) / ( REAL( N )*EPS )
 275          ELSE
 276             IF( ANORM.LT.ONE ) THEN
 277                RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) /
 278      $                 ( REAL( N )*EPS )
 279             ELSE
 280                RESID = MIN( RESID / ANORM, REAL( N ) ) /
 281      $                 ( REAL( N )*EPS )
 282             END IF
 283          END IF
 284       END IF
 285 *
 286       RETURN
 287 *
 288 *     End of SBDT01
 289 *
 290       END