LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DBDT03( 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 DOUBLE PRECISION RESID 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ), 00015 $ VT( LDVT, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DBDT03 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) DOUBLE PRECISION array, dimension (N) 00047 * The n diagonal elements of the bidiagonal matrix B. 00048 * 00049 * E (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N) 00061 * The singular values from the SVD of B, sorted in decreasing 00062 * order. 00063 * 00064 * VT (input) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (2*N) 00072 * 00073 * RESID (output) DOUBLE PRECISION 00074 * The test ratio: norm(B - U * S * V') / ( n * norm(A) * EPS ) 00075 * 00076 * ====================================================================== 00077 * 00078 * .. Parameters .. 00079 DOUBLE PRECISION ZERO, ONE 00080 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00081 * .. 00082 * .. Local Scalars .. 00083 INTEGER I, J 00084 DOUBLE PRECISION BNORM, EPS 00085 * .. 00086 * .. External Functions .. 00087 LOGICAL LSAME 00088 INTEGER IDAMAX 00089 DOUBLE PRECISION DASUM, DLAMCH 00090 EXTERNAL LSAME, IDAMAX, DASUM, DLAMCH 00091 * .. 00092 * .. External Subroutines .. 00093 EXTERNAL DGEMV 00094 * .. 00095 * .. Intrinsic Functions .. 00096 INTRINSIC ABS, DBLE, MAX, MIN 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, 00122 $ WORK( N+1 ), 1, 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, DASUM( 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, 00141 $ WORK( N+1 ), 1, 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, DASUM( 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 DGEMV( 'No transpose', N, N, -ONE, U, LDU, WORK( N+1 ), 00161 $ 1, ZERO, WORK, 1 ) 00162 WORK( J ) = WORK( J ) + D( J ) 00163 RESID = MAX( RESID, DASUM( N, WORK, 1 ) ) 00164 60 CONTINUE 00165 J = IDAMAX( 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 = DLAMCH( '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 ) / ( DBLE( N )*EPS ) 00179 ELSE 00180 IF( BNORM.LT.ONE ) THEN 00181 RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) / 00182 $ ( DBLE( N )*EPS ) 00183 ELSE 00184 RESID = MIN( RESID / BNORM, DBLE( N ) ) / 00185 $ ( DBLE( N )*EPS ) 00186 END IF 00187 END IF 00188 END IF 00189 * 00190 RETURN 00191 * 00192 * End of DBDT03 00193 * 00194 END