LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dbdt04.f
Go to the documentation of this file.
1*> \brief \b DBDT04
2* =========== DOCUMENTATION ===========
3*
4* Online html documentation available at
5* http://www.netlib.org/lapack/explore-html/
6*
7* Definition:
8* ===========
9*
10* SUBROUTINE DBDT04( UPLO, N, D, E, S, NS, U, LDU, VT, LDVT,
11* WORK, RESID )
12*
13* .. Scalar Arguments ..
14* CHARACTER UPLO
15* INTEGER LDU, LDVT, N, NS
16* DOUBLE PRECISION RESID
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ),
20* $ VT( LDVT, * ), WORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DBDT04 reconstructs a bidiagonal matrix B from its (partial) SVD:
30*> S = U' * B * V
31*> where U and V are orthogonal matrices and S is diagonal.
32*>
33*> The test ratio to test the singular value decomposition is
34*> RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
35*> where VT = V' and EPS is the machine precision.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] UPLO
42*> \verbatim
43*> UPLO is CHARACTER*1
44*> Specifies whether the matrix B is upper or lower bidiagonal.
45*> = 'U': Upper bidiagonal
46*> = 'L': Lower bidiagonal
47*> \endverbatim
48*>
49*> \param[in] N
50*> \verbatim
51*> N is INTEGER
52*> The order of the matrix B.
53*> \endverbatim
54*>
55*> \param[in] D
56*> \verbatim
57*> D is DOUBLE PRECISION array, dimension (N)
58*> The n diagonal elements of the bidiagonal matrix B.
59*> \endverbatim
60*>
61*> \param[in] E
62*> \verbatim
63*> E is DOUBLE PRECISION array, dimension (N-1)
64*> The (n-1) superdiagonal elements of the bidiagonal matrix B
65*> if UPLO = 'U', or the (n-1) subdiagonal elements of B if
66*> UPLO = 'L'.
67*> \endverbatim
68*>
69*> \param[in] S
70*> \verbatim
71*> S is DOUBLE PRECISION array, dimension (NS)
72*> The singular values from the (partial) SVD of B, sorted in
73*> decreasing order.
74*> \endverbatim
75*>
76*> \param[in] NS
77*> \verbatim
78*> NS is INTEGER
79*> The number of singular values/vectors from the (partial)
80*> SVD of B.
81*> \endverbatim
82*>
83*> \param[in] U
84*> \verbatim
85*> U is DOUBLE PRECISION array, dimension (LDU,NS)
86*> The n by ns orthogonal matrix U in S = U' * B * V.
87*> \endverbatim
88*>
89*> \param[in] LDU
90*> \verbatim
91*> LDU is INTEGER
92*> The leading dimension of the array U. LDU >= max(1,N)
93*> \endverbatim
94*>
95*> \param[in] VT
96*> \verbatim
97*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
98*> The n by ns orthogonal matrix V in S = U' * B * V.
99*> \endverbatim
100*>
101*> \param[in] LDVT
102*> \verbatim
103*> LDVT is INTEGER
104*> The leading dimension of the array VT.
105*> \endverbatim
106*>
107*> \param[out] WORK
108*> \verbatim
109*> WORK is DOUBLE PRECISION array, dimension (2*N)
110*> \endverbatim
111*>
112*> \param[out] RESID
113*> \verbatim
114*> RESID is DOUBLE PRECISION
115*> The test ratio: norm(S - U' * B * V) / ( n * norm(B) * EPS )
116*> \endverbatim
117*
118* Authors:
119* ========
120*
121*> \author Univ. of Tennessee
122*> \author Univ. of California Berkeley
123*> \author Univ. of Colorado Denver
124*> \author NAG Ltd.
125*
126*> \ingroup double_eig
127*
128* =====================================================================
129 SUBROUTINE dbdt04( UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK,
130 $ RESID )
131*
132* -- LAPACK test routine --
133* -- LAPACK is a software package provided by Univ. of Tennessee, --
134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*
136* .. Scalar Arguments ..
137 CHARACTER UPLO
138 INTEGER LDU, LDVT, N, NS
139 DOUBLE PRECISION RESID
140* ..
141* .. Array Arguments ..
142 DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ),
143 $ vt( ldvt, * ), work( * )
144* ..
145*
146* ======================================================================
147*
148* .. Parameters ..
149 DOUBLE PRECISION ZERO, ONE
150 parameter( zero = 0.0d+0, one = 1.0d+0 )
151* ..
152* .. Local Scalars ..
153 INTEGER I, J, K
154 DOUBLE PRECISION BNORM, EPS
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 INTEGER IDAMAX
159 DOUBLE PRECISION DASUM, DLAMCH
160 EXTERNAL lsame, idamax, dasum, dlamch
161* ..
162* .. External Subroutines ..
163 EXTERNAL dgemm
164* ..
165* .. Intrinsic Functions ..
166 INTRINSIC abs, dble, max, min
167* ..
168* .. Executable Statements ..
169*
170* Quick return if possible.
171*
172 resid = zero
173 IF( n.LE.0 .OR. ns.LE.0 )
174 $ RETURN
175*
176 eps = dlamch( 'Precision' )
177*
178* Compute S - U' * B * V.
179*
180 bnorm = zero
181*
182 IF( lsame( uplo, 'U' ) ) THEN
183*
184* B is upper bidiagonal.
185*
186 k = 0
187 DO 20 i = 1, ns
188 DO 10 j = 1, n-1
189 k = k + 1
190 work( k ) = d( j )*vt( i, j ) + e( j )*vt( i, j+1 )
191 10 CONTINUE
192 k = k + 1
193 work( k ) = d( n )*vt( i, n )
194 20 CONTINUE
195 bnorm = abs( d( 1 ) )
196 DO 30 i = 2, n
197 bnorm = max( bnorm, abs( d( i ) )+abs( e( i-1 ) ) )
198 30 CONTINUE
199 ELSE
200*
201* B is lower bidiagonal.
202*
203 k = 0
204 DO 50 i = 1, ns
205 k = k + 1
206 work( k ) = d( 1 )*vt( i, 1 )
207 DO 40 j = 1, n-1
208 k = k + 1
209 work( k ) = e( j )*vt( i, j ) + d( j+1 )*vt( i, j+1 )
210 40 CONTINUE
211 50 CONTINUE
212 bnorm = abs( d( n ) )
213 DO 60 i = 1, n-1
214 bnorm = max( bnorm, abs( d( i ) )+abs( e( i ) ) )
215 60 CONTINUE
216 END IF
217*
218 CALL dgemm( 'T', 'N', ns, ns, n, -one, u, ldu, work( 1 ),
219 $ n, zero, work( 1+n*ns ), ns )
220*
221* norm(S - U' * B * V)
222*
223 k = n*ns
224 DO 70 i = 1, ns
225 work( k+i ) = work( k+i ) + s( i )
226 resid = max( resid, dasum( ns, work( k+1 ), 1 ) )
227 k = k + ns
228 70 CONTINUE
229*
230 IF( bnorm.LE.zero ) THEN
231 IF( resid.NE.zero )
232 $ resid = one / eps
233 ELSE
234 IF( bnorm.GE.resid ) THEN
235 resid = ( resid / bnorm ) / ( dble( n )*eps )
236 ELSE
237 IF( bnorm.LT.one ) THEN
238 resid = ( min( resid, dble( n )*bnorm ) / bnorm ) /
239 $ ( dble( n )*eps )
240 ELSE
241 resid = min( resid / bnorm, dble( n ) ) /
242 $ ( dble( n )*eps )
243 END IF
244 END IF
245 END IF
246*
247 RETURN
248*
249* End of DBDT04
250*
251 END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dbdt04(UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT04
Definition: dbdt04.f:131