![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DBDT03 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 DBDT03( 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 * DOUBLE PRECISION RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ), 00021 * $ VT( LDVT, * ), WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> DBDT03 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (2*N) 00114 *> \endverbatim 00115 *> 00116 *> \param[out] RESID 00117 *> \verbatim 00118 *> RESID is DOUBLE PRECISION 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 double_eig 00133 * 00134 * ===================================================================== 00135 SUBROUTINE DBDT03( 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 DOUBLE PRECISION RESID 00147 * .. 00148 * .. Array Arguments .. 00149 DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ), 00150 $ VT( LDVT, * ), WORK( * ) 00151 * .. 00152 * 00153 * ====================================================================== 00154 * 00155 * .. Parameters .. 00156 DOUBLE PRECISION ZERO, ONE 00157 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 INTEGER I, J 00161 DOUBLE PRECISION BNORM, EPS 00162 * .. 00163 * .. External Functions .. 00164 LOGICAL LSAME 00165 INTEGER IDAMAX 00166 DOUBLE PRECISION DASUM, DLAMCH 00167 EXTERNAL LSAME, IDAMAX, DASUM, DLAMCH 00168 * .. 00169 * .. External Subroutines .. 00170 EXTERNAL DGEMV 00171 * .. 00172 * .. Intrinsic Functions .. 00173 INTRINSIC ABS, DBLE, MAX, MIN 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, 00199 $ WORK( N+1 ), 1, 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, DASUM( 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, 00218 $ WORK( N+1 ), 1, 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, DASUM( 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, WORK( N+1 ), 00238 $ 1, ZERO, WORK, 1 ) 00239 WORK( J ) = WORK( J ) + D( J ) 00240 RESID = MAX( RESID, DASUM( N, WORK, 1 ) ) 00241 60 CONTINUE 00242 J = IDAMAX( 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 = DLAMCH( '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 ) / ( DBLE( N )*EPS ) 00256 ELSE 00257 IF( BNORM.LT.ONE ) THEN 00258 RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) / 00259 $ ( DBLE( N )*EPS ) 00260 ELSE 00261 RESID = MIN( RESID / BNORM, DBLE( N ) ) / 00262 $ ( DBLE( N )*EPS ) 00263 END IF 00264 END IF 00265 END IF 00266 * 00267 RETURN 00268 * 00269 * End of DBDT03 00270 * 00271 END