TESTING/EIG/dbdt02.f

   1 *> \brief \b DBDT02
   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 DBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RESID )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            LDB, LDC, LDU, M, N
  15 *       DOUBLE PRECISION   RESID
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       DOUBLE PRECISION   B( LDB, * ), C( LDC, * ), U( LDU, * ),
  19 *      $                   WORK( * )
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> DBDT02 tests the change of basis C = U' * B by computing the residual
  29 *>
  30 *>    RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
  31 *>
  32 *> where B and C are M by N matrices, U is an M by M orthogonal matrix,
  33 *> 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 B and C and the order of
  43 *>          the matrix Q.
  44 *> \endverbatim
  45 *>
  46 *> \param[in] N
  47 *> \verbatim
  48 *>          N is INTEGER
  49 *>          The number of columns of the matrices B and C.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] B
  53 *> \verbatim
  54 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
  55 *>          The m by n matrix B.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] LDB
  59 *> \verbatim
  60 *>          LDB is INTEGER
  61 *>          The leading dimension of the array B.  LDB >= max(1,M).
  62 *> \endverbatim
  63 *>
  64 *> \param[in] C
  65 *> \verbatim
  66 *>          C is DOUBLE PRECISION array, dimension (LDC,N)
  67 *>          The m by n matrix C, assumed to contain U' * B.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] LDC
  71 *> \verbatim
  72 *>          LDC is INTEGER
  73 *>          The leading dimension of the array C.  LDC >= max(1,M).
  74 *> \endverbatim
  75 *>
  76 *> \param[in] U
  77 *> \verbatim
  78 *>          U is DOUBLE PRECISION array, dimension (LDU,M)
  79 *>          The m by m orthogonal matrix U.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LDU
  83 *> \verbatim
  84 *>          LDU is INTEGER
  85 *>          The leading dimension of the array U.  LDU >= max(1,M).
  86 *> \endverbatim
  87 *>
  88 *> \param[out] WORK
  89 *> \verbatim
  90 *>          WORK is DOUBLE PRECISION array, dimension (M)
  91 *> \endverbatim
  92 *>
  93 *> \param[out] RESID
  94 *> \verbatim
  95 *>          RESID is DOUBLE PRECISION
  96 *>          RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
  97 *> \endverbatim
  98 *
  99 *  Authors:
 100 *  ========
 101 *
 102 *> \author Univ. of Tennessee
 103 *> \author Univ. of California Berkeley
 104 *> \author Univ. of Colorado Denver
 105 *> \author NAG Ltd.
 106 *
 107 *> \date November 2011
 108 *
 109 *> \ingroup double_eig
 110 *
 111 *  =====================================================================
 112       SUBROUTINE DBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RESID )
 113 *
 114 *  -- LAPACK test routine (version 3.4.0) --
 115 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 116 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 117 *     November 2011
 118 *
 119 *     .. Scalar Arguments ..
 120       INTEGER            LDB, LDC, LDU, M, N
 121       DOUBLE PRECISION   RESID
 122 *     ..
 123 *     .. Array Arguments ..
 124       DOUBLE PRECISION   B( LDB, * ), C( LDC, * ), U( LDU, * ),
 125      $                   WORK( * )
 126 *     ..
 127 *
 128 * ======================================================================
 129 *
 130 *     .. Parameters ..
 131       DOUBLE PRECISION   ZERO, ONE
 132       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 133 *     ..
 134 *     .. Local Scalars ..
 135       INTEGER            J
 136       DOUBLE PRECISION   BNORM, EPS, REALMN
 137 *     ..
 138 *     .. External Functions ..
 139       DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
 140       EXTERNAL           DASUM, DLAMCH, DLANGE
 141 *     ..
 142 *     .. External Subroutines ..
 143       EXTERNAL           DCOPY, DGEMV
 144 *     ..
 145 *     .. Intrinsic Functions ..
 146       INTRINSIC          DBLE, MAX, MIN
 147 *     ..
 148 *     .. Executable Statements ..
 149 *
 150 *     Quick return if possible
 151 *
 152       RESID = ZERO
 153       IF( M.LE.0 .OR. N.LE.0 )
 154      $   RETURN
 155       REALMN = DBLE( MAX( M, N ) )
 156       EPS = DLAMCH( 'Precision' )
 157 *
 158 *     Compute norm( B - U * C )
 159 *
 160       DO 10 J = 1, N
 161          CALL DCOPY( M, B( 1, J ), 1, WORK, 1 )
 162          CALL DGEMV( 'No transpose', M, M, -ONE, U, LDU, C( 1, J ), 1,
 163      $               ONE, WORK, 1 )
 164          RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
 165    10 CONTINUE
 166 *
 167 *     Compute norm of B.
 168 *
 169       BNORM = DLANGE( '1', M, N, B, LDB, WORK )
 170 *
 171       IF( BNORM.LE.ZERO ) THEN
 172          IF( RESID.NE.ZERO )
 173      $      RESID = ONE / EPS
 174       ELSE
 175          IF( BNORM.GE.RESID ) THEN
 176             RESID = ( RESID / BNORM ) / ( REALMN*EPS )
 177          ELSE
 178             IF( BNORM.LT.ONE ) THEN
 179                RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) /
 180      $                 ( REALMN*EPS )
 181             ELSE
 182                RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS )
 183             END IF
 184          END IF
 185       END IF
 186       RETURN
 187 *
 188 *     End of DBDT02
 189 *
 190       END