LAPACK: TESTING/EIG/cbdt03.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
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       CHARACTER          UPLO
00010       INTEGER            KD, LDU, LDVT, N
00011       REAL               RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               D( * ), E( * ), S( * )
00015       COMPLEX            U( LDU, * ), VT( LDVT, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CBDT03 reconstructs a bidiagonal matrix B from its SVD:
00022 *     S = U' * B * V
00023 *  where U and V are orthogonal matrices and S is diagonal.
00024 *
00025 *  The test ratio to test the singular value decomposition is
00026 *     RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
00027 *  where VT = V' and EPS is the machine precision.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  UPLO    (input) CHARACTER*1
00033 *          Specifies whether the matrix B is upper or lower bidiagonal.
00034 *          = 'U':  Upper bidiagonal
00035 *          = 'L':  Lower bidiagonal
00036 *
00037 *  N       (input) INTEGER
00038 *          The order of the matrix B.
00039 *
00040 *  KD      (input) INTEGER
00041 *          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
00042 *          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
00043 *          not referenced.  If KD is greater than 1, it is assumed to be
00044 *          1, and if KD is less than 0, it is assumed to be 0.
00045 *
00046 *  D       (input) REAL array, dimension (N)
00047 *          The n diagonal elements of the bidiagonal matrix B.
00048 *
00049 *  E       (input) REAL array, dimension (N-1)
00050 *          The (n-1) superdiagonal elements of the bidiagonal matrix B
00051 *          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
00052 *          UPLO = 'L'.
00053 *
00054 *  U       (input) COMPLEX array, dimension (LDU,N)
00055 *          The n by n orthogonal matrix U in the reduction B = U'*A*P.
00056 *
00057 *  LDU     (input) INTEGER
00058 *          The leading dimension of the array U.  LDU >= max(1,N)
00059 *
00060 *  S       (input) REAL array, dimension (N)
00061 *          The singular values from the SVD of B, sorted in decreasing
00062 *          order.
00063 *
00064 *  VT      (input) COMPLEX array, dimension (LDVT,N)
00065 *          The n by n orthogonal matrix V' in the reduction
00066 *          B = U * S * V'.
00067 *
00068 *  LDVT    (input) INTEGER
00069 *          The leading dimension of the array VT.
00070 *
00071 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00072 *
00073 *  RESID   (output) REAL
00074 *          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
00075 *
00076 * ======================================================================
00077 *
00078 *     .. Parameters ..
00079       REAL               ZERO, ONE
00080       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       INTEGER            I, J
00084       REAL               BNORM, EPS
00085 *     ..
00086 *     .. External Functions ..
00087       LOGICAL            LSAME
00088       INTEGER            ISAMAX
00089       REAL               SCASUM, SLAMCH
00090       EXTERNAL           LSAME, ISAMAX, SCASUM, SLAMCH
00091 *     ..
00092 *     .. External Subroutines ..
00093       EXTERNAL           CGEMV
00094 *     ..
00095 *     .. Intrinsic Functions ..
00096       INTRINSIC          ABS, CMPLX, MAX, MIN, REAL
00097 *     ..
00098 *     .. Executable Statements ..
00099 *
00100 *     Quick return if possible
00101 *
00102       RESID = ZERO
00103       IF( N.LE.0 )
00104      $   RETURN
00105 *
00106 *     Compute B - U * S * V' one column at a time.
00107 *
00108       BNORM = ZERO
00109       IF( KD.GE.1 ) THEN
00110 *
00111 *        B is bidiagonal.
00112 *
00113          IF( LSAME( UPLO, 'U' ) ) THEN
00114 *
00115 *           B is upper bidiagonal.
00116 *
00117             DO 20 J = 1, N
00118                DO 10 I = 1, N
00119                   WORK( N+I ) = S( I )*VT( I, J )
00120    10          CONTINUE
00121                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00122      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00123                WORK( J ) = WORK( J ) + D( J )
00124                IF( J.GT.1 ) THEN
00125                   WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
00126                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
00127                ELSE
00128                   BNORM = MAX( BNORM, ABS( D( J ) ) )
00129                END IF
00130                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00131    20       CONTINUE
00132          ELSE
00133 *
00134 *           B is lower bidiagonal.
00135 *
00136             DO 40 J = 1, N
00137                DO 30 I = 1, N
00138                   WORK( N+I ) = S( I )*VT( I, J )
00139    30          CONTINUE
00140                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00141      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00142                WORK( J ) = WORK( J ) + D( J )
00143                IF( J.LT.N ) THEN
00144                   WORK( J+1 ) = WORK( J+1 ) + E( J )
00145                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
00146                ELSE
00147                   BNORM = MAX( BNORM, ABS( D( J ) ) )
00148                END IF
00149                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00150    40       CONTINUE
00151          END IF
00152       ELSE
00153 *
00154 *        B is diagonal.
00155 *
00156          DO 60 J = 1, N
00157             DO 50 I = 1, N
00158                WORK( N+I ) = S( I )*VT( I, J )
00159    50       CONTINUE
00160             CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00161      $                  WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00162             WORK( J ) = WORK( J ) + D( J )
00163             RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00164    60    CONTINUE
00165          J = ISAMAX( N, D, 1 )
00166          BNORM = ABS( D( J ) )
00167       END IF
00168 *
00169 *     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
00170 *
00171       EPS = SLAMCH( 'Precision' )
00172 *
00173       IF( BNORM.LE.ZERO ) THEN
00174          IF( RESID.NE.ZERO )
00175      $      RESID = ONE / EPS
00176       ELSE
00177          IF( BNORM.GE.RESID ) THEN
00178             RESID = ( RESID / BNORM ) / ( REAL( N )*EPS )
00179          ELSE
00180             IF( BNORM.LT.ONE ) THEN
00181                RESID = ( MIN( RESID, REAL( N )*BNORM ) / BNORM ) /
00182      $                 ( REAL( N )*EPS )
00183             ELSE
00184                RESID = MIN( RESID / BNORM, REAL( N ) ) /
00185      $                 ( REAL( N )*EPS )
00186             END IF
00187          END IF
00188       END IF
00189 *
00190       RETURN
00191 *
00192 *     End of CBDT03
00193 *
00194       END