1       SUBROUTINE SBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RESID )
  2 *
  3 *  -- LAPACK test routine (version 3.1) --
  4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  5 *     November 2006
  6 *
  7 *     .. Scalar Arguments ..
  8       INTEGER            LDB, LDC, LDU, M, N
  9       REAL               RESID
 10 *     ..
 11 *     .. Array Arguments ..
 12       REAL               B( LDB, * ), C( LDC, * ), U( LDU, * ),
 13      $                   WORK( * )
 14 *     ..
 15 *
 16 *  Purpose
 17 *  =======
 18 *
 19 *  SBDT02 tests the change of basis C = U' * B by computing the residual
 20 *
 21 *     RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
 22 *
 23 *  where B and C are M by N matrices, U is an M by M orthogonal matrix,
 24 *  and EPS is the machine precision.
 25 *
 26 *  Arguments
 27 *  =========
 28 *
 29 *  M       (input) INTEGER
 30 *          The number of rows of the matrices B and C and the order of
 31 *          the matrix Q.
 32 *
 33 *  N       (input) INTEGER
 34 *          The number of columns of the matrices B and C.
 35 *
 36 *  B       (input) REAL array, dimension (LDB,N)
 37 *          The m by n matrix B.
 38 *
 39 *  LDB     (input) INTEGER
 40 *          The leading dimension of the array B.  LDB >= max(1,M).
 41 *
 42 *  C       (input) REAL array, dimension (LDC,N)
 43 *          The m by n matrix C, assumed to contain U' * B.
 44 *
 45 *  LDC     (input) INTEGER
 46 *          The leading dimension of the array C.  LDC >= max(1,M).
 47 *
 48 *  U       (input) REAL array, dimension (LDU,M)
 49 *          The m by m orthogonal matrix U.
 50 *
 51 *  LDU     (input) INTEGER
 52 *          The leading dimension of the array U.  LDU >= max(1,M).
 53 *
 54 *  WORK    (workspace) REAL array, dimension (M)
 55 *
 56 *  RESID   (output) REAL
 57 *          RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
 58 *
 59 * ======================================================================
 60 *
 61 *     .. Parameters ..
 62       REAL               ZERO, ONE
 63       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 64 *     ..
 65 *     .. Local Scalars ..
 66       INTEGER            J
 67       REAL               BNORM, EPS, REALMN
 68 *     ..
 69 *     .. External Functions ..
 70       REAL               SASUM, SLAMCH, SLANGE
 71       EXTERNAL           SASUM, SLAMCH, SLANGE
 72 *     ..
 73 *     .. External Subroutines ..
 74       EXTERNAL           SCOPY, SGEMV
 75 *     ..
 76 *     .. Intrinsic Functions ..
 77       INTRINSIC          MAXMIN, REAL
 78 *     ..
 79 *     .. Executable Statements ..
 80 *
 81 *     Quick return if possible
 82 *
 83       RESID = ZERO
 84       IF( M.LE.0 .OR. N.LE.0 )
 85      $   RETURN
 86       REALMN = REALMAX( M, N ) )
 87       EPS = SLAMCH( 'Precision' )
 88 *
 89 *     Compute norm( B - U * C )
 90 *
 91       DO 10 J = 1, N
 92          CALL SCOPY( M, B( 1, J ), 1, WORK, 1 )
 93          CALL SGEMV( 'No transpose', M, M, -ONE, U, LDU, C( 1, J ), 1,
 94      $               ONE, WORK, 1 )
 95          RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
 96    10 CONTINUE
 97 *
 98 *     Compute norm of B.
 99 *
100       BNORM = SLANGE( '1', M, N, B, LDB, WORK )
101 *
102       IF( BNORM.LE.ZERO ) THEN
103          IF( RESID.NE.ZERO )
104      $      RESID = ONE / EPS
105       ELSE
106          IF( BNORM.GE.RESID ) THEN
107             RESID = ( RESID / BNORM ) / ( REALMN*EPS )
108          ELSE
109             IF( BNORM.LT.ONE ) THEN
110                RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) /
111      $                 ( REALMN*EPS )
112             ELSE
113                RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS )
114             END IF
115          END IF
116       END IF
117       RETURN
118 *
119 *     End of SBDT02
120 *
121       END