LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dbdt05.f
Go to the documentation of this file.
1*> \brief \b DBDT05
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DBDT05( M, N, A, LDA, S, NS, U, LDU,
12* VT, LDVT, WORK, RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER LDA, 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*> DBDT05 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] M
42*> \verbatim
43*> M is INTEGER
44*> The number of rows of the matrices A and U.
45*> \endverbatim
46*>
47*> \param[in] N
48*> \verbatim
49*> N is INTEGER
50*> The number of columns of the matrices A and VT.
51*> \endverbatim
52*>
53*> \param[in] A
54*> \verbatim
55*> A is DOUBLE PRECISION array, dimension (LDA,N)
56*> The m by n matrix A.
57*> \endverbatim
58*>
59*> \param[in] LDA
60*> \verbatim
61*> LDA is INTEGER
62*> The leading dimension of the array A. LDA >= max(1,M).
63*> \endverbatim
64*>
65*> \param[in] S
66*> \verbatim
67*> S is DOUBLE PRECISION array, dimension (NS)
68*> The singular values from the (partial) SVD of B, sorted in
69*> decreasing order.
70*> \endverbatim
71*>
72*> \param[in] NS
73*> \verbatim
74*> NS is INTEGER
75*> The number of singular values/vectors from the (partial)
76*> SVD of B.
77*> \endverbatim
78*>
79*> \param[in] U
80*> \verbatim
81*> U is DOUBLE PRECISION array, dimension (LDU,NS)
82*> The n by ns orthogonal matrix U in S = U' * B * V.
83*> \endverbatim
84*>
85*> \param[in] LDU
86*> \verbatim
87*> LDU is INTEGER
88*> The leading dimension of the array U. LDU >= max(1,N)
89*> \endverbatim
90*>
91*> \param[in] VT
92*> \verbatim
93*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
94*> The n by ns orthogonal matrix V in S = U' * B * V.
95*> \endverbatim
96*>
97*> \param[in] LDVT
98*> \verbatim
99*> LDVT is INTEGER
100*> The leading dimension of the array VT.
101*> \endverbatim
102*>
103*> \param[out] WORK
104*> \verbatim
105*> WORK is DOUBLE PRECISION array, dimension (M,N)
106*> \endverbatim
107*>
108*> \param[out] RESID
109*> \verbatim
110*> RESID is DOUBLE PRECISION
111*> The test ratio: norm(S - U' * A * V) / ( n * norm(A) * EPS )
112*> \endverbatim
113*
114* Authors:
115* ========
116*
117*> \author Univ. of Tennessee
118*> \author Univ. of California Berkeley
119*> \author Univ. of Colorado Denver
120*> \author NAG Ltd.
121*
122*> \ingroup double_eig
123*
124* =====================================================================
125 SUBROUTINE dbdt05( M, N, A, LDA, S, NS, U, LDU,
126 $ VT, LDVT, WORK, RESID )
127*
128* -- LAPACK test routine --
129* -- LAPACK is a software package provided by Univ. of Tennessee, --
130* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131*
132* .. Scalar Arguments ..
133 INTEGER LDA, LDU, LDVT, M, N, NS
134 DOUBLE PRECISION RESID
135* ..
136* .. Array Arguments ..
137 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
138 $ vt( ldvt, * ), work( * )
139* ..
140*
141* ======================================================================
142*
143* .. Parameters ..
144 DOUBLE PRECISION ZERO, ONE
145 parameter( zero = 0.0d+0, one = 1.0d+0 )
146* ..
147* .. Local Scalars ..
148 INTEGER I, J
149 DOUBLE PRECISION ANORM, EPS
150* ..
151* .. External Functions ..
152 LOGICAL LSAME
153 INTEGER IDAMAX
154 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
155 EXTERNAL lsame, idamax, dasum, dlamch, dlange
156* ..
157* .. External Subroutines ..
158 EXTERNAL dgemm
159* ..
160* .. Intrinsic Functions ..
161 INTRINSIC abs, dble, max, min
162* ..
163* .. Executable Statements ..
164*
165* Quick return if possible.
166*
167 resid = zero
168 IF( min( m, n ).LE.0 .OR. ns.LE.0 )
169 $ RETURN
170*
171 eps = dlamch( 'Precision' )
172 anorm = dlange( 'M', m, n, a, lda, work )
173*
174* Compute U' * A * V.
175*
176 CALL dgemm( 'N', 'T', m, ns, n, one, a, lda, vt,
177 $ ldvt, zero, work( 1+ns*ns ), m )
178 CALL dgemm( 'T', 'N', ns, ns, m, -one, u, ldu, work( 1+ns*ns ),
179 $ m, zero, work, ns )
180*
181* norm(S - U' * B * V)
182*
183 j = 0
184 DO 10 i = 1, ns
185 work( j+i ) = work( j+i ) + s( i )
186 resid = max( resid, dasum( ns, work( j+1 ), 1 ) )
187 j = j + ns
188 10 CONTINUE
189*
190 IF( anorm.LE.zero ) THEN
191 IF( resid.NE.zero )
192 $ resid = one / eps
193 ELSE
194 IF( anorm.GE.resid ) THEN
195 resid = ( resid / anorm ) / ( dble( n )*eps )
196 ELSE
197 IF( anorm.LT.one ) THEN
198 resid = ( min( resid, dble( n )*anorm ) / anorm ) /
199 $ ( dble( n )*eps )
200 ELSE
201 resid = min( resid / anorm, dble( n ) ) /
202 $ ( dble( n )*eps )
203 END IF
204 END IF
205 END IF
206*
207 RETURN
208*
209* End of DBDT05
210*
211 END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT05
Definition: dbdt05.f:127