LAPACK 3.3.1 Linear Algebra PACKage

zbdt02.f

Go to the documentation of this file.
```00001       SUBROUTINE ZBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RWORK,
00002      \$                   RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            LDB, LDC, LDU, M, N
00010       DOUBLE PRECISION   RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   RWORK( * )
00014       COMPLEX*16         B( LDB, * ), C( LDC, * ), U( LDU, * ),
00015      \$                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZBDT02 tests the change of basis C = U' * B by computing the residual
00022 *
00023 *     RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
00024 *
00025 *  where B and C are M by N matrices, U is an M by M orthogonal matrix,
00026 *  and EPS is the machine precision.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  M       (input) INTEGER
00032 *          The number of rows of the matrices B and C and the order of
00033 *          the matrix Q.
00034 *
00035 *  N       (input) INTEGER
00036 *          The number of columns of the matrices B and C.
00037 *
00038 *  B       (input) COMPLEX*16 array, dimension (LDB,N)
00039 *          The m by n matrix B.
00040 *
00041 *  LDB     (input) INTEGER
00042 *          The leading dimension of the array B.  LDB >= max(1,M).
00043 *
00044 *  C       (input) COMPLEX*16 array, dimension (LDC,N)
00045 *          The m by n matrix C, assumed to contain U' * B.
00046 *
00047 *  LDC     (input) INTEGER
00048 *          The leading dimension of the array C.  LDC >= max(1,M).
00049 *
00050 *  U       (input) COMPLEX*16 array, dimension (LDU,M)
00051 *          The m by m orthogonal matrix U.
00052 *
00053 *  LDU     (input) INTEGER
00054 *          The leading dimension of the array U.  LDU >= max(1,M).
00055 *
00056 *  WORK    (workspace) COMPLEX*16 array, dimension (M)
00057 *
00058 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (M)
00059 *
00060 *  RESID   (output) DOUBLE PRECISION
00061 *          RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
00062 *
00063 * ======================================================================
00064 *
00065 *     .. Parameters ..
00066       DOUBLE PRECISION   ZERO, ONE
00067       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00068 *     ..
00069 *     .. Local Scalars ..
00070       INTEGER            J
00071       DOUBLE PRECISION   BNORM, EPS, REALMN
00072 *     ..
00073 *     .. External Functions ..
00074       DOUBLE PRECISION   DLAMCH, DZASUM, ZLANGE
00075       EXTERNAL           DLAMCH, DZASUM, ZLANGE
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           ZCOPY, ZGEMV
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          DBLE, DCMPLX, MAX, MIN
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Quick return if possible
00086 *
00087       RESID = ZERO
00088       IF( M.LE.0 .OR. N.LE.0 )
00089      \$   RETURN
00090       REALMN = DBLE( MAX( M, N ) )
00091       EPS = DLAMCH( 'Precision' )
00092 *
00093 *     Compute norm( B - U * C )
00094 *
00095       DO 10 J = 1, N
00096          CALL ZCOPY( M, B( 1, J ), 1, WORK, 1 )
00097          CALL ZGEMV( 'No transpose', M, M, -DCMPLX( ONE ), U, LDU,
00098      \$               C( 1, J ), 1, DCMPLX( ONE ), WORK, 1 )
00099          RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00100    10 CONTINUE
00101 *
00102 *     Compute norm of B.
00103 *
00104       BNORM = ZLANGE( '1', M, N, B, LDB, RWORK )
00105 *
00106       IF( BNORM.LE.ZERO ) THEN
00107          IF( RESID.NE.ZERO )
00108      \$      RESID = ONE / EPS
00109       ELSE
00110          IF( BNORM.GE.RESID ) THEN
00111             RESID = ( RESID / BNORM ) / ( REALMN*EPS )
00112          ELSE
00113             IF( BNORM.LT.ONE ) THEN
00114                RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) /
00115      \$                 ( REALMN*EPS )
00116             ELSE
00117                RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS )
00118             END IF
00119          END IF
00120       END IF
00121       RETURN
00122 *
00123 *     End of ZBDT02
00124 *
00125       END
```