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