LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cbdt03.f
Go to the documentation of this file.
00001 *> \brief \b CBDT03
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            KD, LDU, LDVT, N
00017 *       REAL               RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               D( * ), E( * ), S( * )
00021 *       COMPLEX            U( LDU, * ), VT( LDVT, * ), WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CBDT03 reconstructs a bidiagonal matrix B from its SVD:
00031 *>    S = U' * B * V
00032 *> where U and V are orthogonal matrices and S is diagonal.
00033 *>
00034 *> The test ratio to test the singular value decomposition is
00035 *>    RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
00036 *> where VT = V' and EPS is the machine precision.
00037 *> \endverbatim
00038 *
00039 *  Arguments:
00040 *  ==========
00041 *
00042 *> \param[in] UPLO
00043 *> \verbatim
00044 *>          UPLO is CHARACTER*1
00045 *>          Specifies whether the matrix B is upper or lower bidiagonal.
00046 *>          = 'U':  Upper bidiagonal
00047 *>          = 'L':  Lower bidiagonal
00048 *> \endverbatim
00049 *>
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The order of the matrix B.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] KD
00057 *> \verbatim
00058 *>          KD is INTEGER
00059 *>          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
00060 *>          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
00061 *>          not referenced.  If KD is greater than 1, it is assumed to be
00062 *>          1, and if KD is less than 0, it is assumed to be 0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] D
00066 *> \verbatim
00067 *>          D is REAL array, dimension (N)
00068 *>          The n diagonal elements of the bidiagonal matrix B.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] E
00072 *> \verbatim
00073 *>          E is REAL array, dimension (N-1)
00074 *>          The (n-1) superdiagonal elements of the bidiagonal matrix B
00075 *>          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
00076 *>          UPLO = 'L'.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] U
00080 *> \verbatim
00081 *>          U is COMPLEX array, dimension (LDU,N)
00082 *>          The n by n orthogonal matrix U in the reduction B = U'*A*P.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] LDU
00086 *> \verbatim
00087 *>          LDU is INTEGER
00088 *>          The leading dimension of the array U.  LDU >= max(1,N)
00089 *> \endverbatim
00090 *>
00091 *> \param[in] S
00092 *> \verbatim
00093 *>          S is REAL array, dimension (N)
00094 *>          The singular values from the SVD of B, sorted in decreasing
00095 *>          order.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] VT
00099 *> \verbatim
00100 *>          VT is COMPLEX array, dimension (LDVT,N)
00101 *>          The n by n orthogonal matrix V' in the reduction
00102 *>          B = U * S * V'.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] LDVT
00106 *> \verbatim
00107 *>          LDVT is INTEGER
00108 *>          The leading dimension of the array VT.
00109 *> \endverbatim
00110 *>
00111 *> \param[out] WORK
00112 *> \verbatim
00113 *>          WORK is COMPLEX array, dimension (2*N)
00114 *> \endverbatim
00115 *>
00116 *> \param[out] RESID
00117 *> \verbatim
00118 *>          RESID is REAL
00119 *>          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
00120 *> \endverbatim
00121 *
00122 *  Authors:
00123 *  ========
00124 *
00125 *> \author Univ. of Tennessee 
00126 *> \author Univ. of California Berkeley 
00127 *> \author Univ. of Colorado Denver 
00128 *> \author NAG Ltd. 
00129 *
00130 *> \date November 2011
00131 *
00132 *> \ingroup complex_eig
00133 *
00134 *  =====================================================================
00135       SUBROUTINE CBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
00136      $                   RESID )
00137 *
00138 *  -- LAPACK test routine (version 3.4.0) --
00139 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00140 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00141 *     November 2011
00142 *
00143 *     .. Scalar Arguments ..
00144       CHARACTER          UPLO
00145       INTEGER            KD, LDU, LDVT, N
00146       REAL               RESID
00147 *     ..
00148 *     .. Array Arguments ..
00149       REAL               D( * ), E( * ), S( * )
00150       COMPLEX            U( LDU, * ), VT( LDVT, * ), WORK( * )
00151 *     ..
00152 *
00153 * ======================================================================
00154 *
00155 *     .. Parameters ..
00156       REAL               ZERO, ONE
00157       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00158 *     ..
00159 *     .. Local Scalars ..
00160       INTEGER            I, J
00161       REAL               BNORM, EPS
00162 *     ..
00163 *     .. External Functions ..
00164       LOGICAL            LSAME
00165       INTEGER            ISAMAX
00166       REAL               SCASUM, SLAMCH
00167       EXTERNAL           LSAME, ISAMAX, SCASUM, SLAMCH
00168 *     ..
00169 *     .. External Subroutines ..
00170       EXTERNAL           CGEMV
00171 *     ..
00172 *     .. Intrinsic Functions ..
00173       INTRINSIC          ABS, CMPLX, MAX, MIN, REAL
00174 *     ..
00175 *     .. Executable Statements ..
00176 *
00177 *     Quick return if possible
00178 *
00179       RESID = ZERO
00180       IF( N.LE.0 )
00181      $   RETURN
00182 *
00183 *     Compute B - U * S * V' one column at a time.
00184 *
00185       BNORM = ZERO
00186       IF( KD.GE.1 ) THEN
00187 *
00188 *        B is bidiagonal.
00189 *
00190          IF( LSAME( UPLO, 'U' ) ) THEN
00191 *
00192 *           B is upper bidiagonal.
00193 *
00194             DO 20 J = 1, N
00195                DO 10 I = 1, N
00196                   WORK( N+I ) = S( I )*VT( I, J )
00197    10          CONTINUE
00198                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00199      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00200                WORK( J ) = WORK( J ) + D( J )
00201                IF( J.GT.1 ) THEN
00202                   WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
00203                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
00204                ELSE
00205                   BNORM = MAX( BNORM, ABS( D( J ) ) )
00206                END IF
00207                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00208    20       CONTINUE
00209          ELSE
00210 *
00211 *           B is lower bidiagonal.
00212 *
00213             DO 40 J = 1, N
00214                DO 30 I = 1, N
00215                   WORK( N+I ) = S( I )*VT( I, J )
00216    30          CONTINUE
00217                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00218      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00219                WORK( J ) = WORK( J ) + D( J )
00220                IF( J.LT.N ) THEN
00221                   WORK( J+1 ) = WORK( J+1 ) + E( J )
00222                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
00223                ELSE
00224                   BNORM = MAX( BNORM, ABS( D( J ) ) )
00225                END IF
00226                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00227    40       CONTINUE
00228          END IF
00229       ELSE
00230 *
00231 *        B is diagonal.
00232 *
00233          DO 60 J = 1, N
00234             DO 50 I = 1, N
00235                WORK( N+I ) = S( I )*VT( I, J )
00236    50       CONTINUE
00237             CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
00238      $                  WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
00239             WORK( J ) = WORK( J ) + D( J )
00240             RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
00241    60    CONTINUE
00242          J = ISAMAX( N, D, 1 )
00243          BNORM = ABS( D( J ) )
00244       END IF
00245 *
00246 *     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
00247 *
00248       EPS = SLAMCH( 'Precision' )
00249 *
00250       IF( BNORM.LE.ZERO ) THEN
00251          IF( RESID.NE.ZERO )
00252      $      RESID = ONE / EPS
00253       ELSE
00254          IF( BNORM.GE.RESID ) THEN
00255             RESID = ( RESID / BNORM ) / ( REAL( N )*EPS )
00256          ELSE
00257             IF( BNORM.LT.ONE ) THEN
00258                RESID = ( MIN( RESID, REAL( N )*BNORM ) / BNORM ) /
00259      $                 ( REAL( N )*EPS )
00260             ELSE
00261                RESID = MIN( RESID / BNORM, REAL( N ) ) /
00262      $                 ( REAL( N )*EPS )
00263             END IF
00264          END IF
00265       END IF
00266 *
00267       RETURN
00268 *
00269 *     End of CBDT03
00270 *
00271       END
 All Files Functions