LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cbdt01.f
Go to the documentation of this file.
1*> \brief \b CBDT01
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 CBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
12* RWORK, RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER KD, LDA, LDPT, LDQ, M, N
16* REAL RESID
17* ..
18* .. Array Arguments ..
19* REAL D( * ), E( * ), RWORK( * )
20* COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
21* $ WORK( * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> CBDT01 reconstructs a general matrix A from its bidiagonal form
31*> A = Q * B * P**H
32*> where Q (m by min(m,n)) and P**H (min(m,n) by n) are unitary
33*> matrices and B is bidiagonal.
34*>
35*> The test ratio to test the reduction is
36*> RESID = norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
37*> where EPS is the machine precision.
38*> \endverbatim
39*
40* Arguments:
41* ==========
42*
43*> \param[in] M
44*> \verbatim
45*> M is INTEGER
46*> The number of rows of the matrices A and Q.
47*> \endverbatim
48*>
49*> \param[in] N
50*> \verbatim
51*> N is INTEGER
52*> The number of columns of the matrices A and P**H.
53*> \endverbatim
54*>
55*> \param[in] KD
56*> \verbatim
57*> KD is INTEGER
58*> If KD = 0, B is diagonal and the array E is not referenced.
59*> If KD = 1, the reduction was performed by xGEBRD; B is upper
60*> bidiagonal if M >= N, and lower bidiagonal if M < N.
61*> If KD = -1, the reduction was performed by xGBBRD; B is
62*> always upper bidiagonal.
63*> \endverbatim
64*>
65*> \param[in] A
66*> \verbatim
67*> A is COMPLEX array, dimension (LDA,N)
68*> The m by n matrix A.
69*> \endverbatim
70*>
71*> \param[in] LDA
72*> \verbatim
73*> LDA is INTEGER
74*> The leading dimension of the array A. LDA >= max(1,M).
75*> \endverbatim
76*>
77*> \param[in] Q
78*> \verbatim
79*> Q is COMPLEX array, dimension (LDQ,N)
80*> The m by min(m,n) unitary matrix Q in the reduction
81*> A = Q * B * P**H.
82*> \endverbatim
83*>
84*> \param[in] LDQ
85*> \verbatim
86*> LDQ is INTEGER
87*> The leading dimension of the array Q. LDQ >= max(1,M).
88*> \endverbatim
89*>
90*> \param[in] D
91*> \verbatim
92*> D is REAL array, dimension (min(M,N))
93*> The diagonal elements of the bidiagonal matrix B.
94*> \endverbatim
95*>
96*> \param[in] E
97*> \verbatim
98*> E is REAL array, dimension (min(M,N)-1)
99*> The superdiagonal elements of the bidiagonal matrix B if
100*> m >= n, or the subdiagonal elements of B if m < n.
101*> \endverbatim
102*>
103*> \param[in] PT
104*> \verbatim
105*> PT is COMPLEX array, dimension (LDPT,N)
106*> The min(m,n) by n unitary matrix P**H in the reduction
107*> A = Q * B * P**H.
108*> \endverbatim
109*>
110*> \param[in] LDPT
111*> \verbatim
112*> LDPT is INTEGER
113*> The leading dimension of the array PT.
114*> LDPT >= max(1,min(M,N)).
115*> \endverbatim
116*>
117*> \param[out] WORK
118*> \verbatim
119*> WORK is COMPLEX array, dimension (M+N)
120*> \endverbatim
121*>
122*> \param[out] RWORK
123*> \verbatim
124*> RWORK is REAL array, dimension (M)
125*> \endverbatim
126*>
127*> \param[out] RESID
128*> \verbatim
129*> RESID is REAL
130*> The test ratio:
131*> norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
132*> \endverbatim
133*
134* Authors:
135* ========
136*
137*> \author Univ. of Tennessee
138*> \author Univ. of California Berkeley
139*> \author Univ. of Colorado Denver
140*> \author NAG Ltd.
141*
142*> \ingroup complex_eig
143*
144* =====================================================================
145 SUBROUTINE cbdt01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
146 $ RWORK, RESID )
147*
148* -- LAPACK test routine --
149* -- LAPACK is a software package provided by Univ. of Tennessee, --
150* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151*
152* .. Scalar Arguments ..
153 INTEGER KD, LDA, LDPT, LDQ, M, N
154 REAL RESID
155* ..
156* .. Array Arguments ..
157 REAL D( * ), E( * ), RWORK( * )
158 COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
159 $ work( * )
160* ..
161*
162* =====================================================================
163*
164* .. Parameters ..
165 REAL ZERO, ONE
166 parameter( zero = 0.0e+0, one = 1.0e+0 )
167* ..
168* .. Local Scalars ..
169 INTEGER I, J
170 REAL ANORM, EPS
171* ..
172* .. External Functions ..
173 REAL CLANGE, SCASUM, SLAMCH
174 EXTERNAL clange, scasum, slamch
175* ..
176* .. External Subroutines ..
177 EXTERNAL ccopy, cgemv
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC cmplx, max, min, real
181* ..
182* .. Executable Statements ..
183*
184* Quick return if possible
185*
186 IF( m.LE.0 .OR. n.LE.0 ) THEN
187 resid = zero
188 RETURN
189 END IF
190*
191* Compute A - Q * B * P**H one column at a time.
192*
193 resid = zero
194 IF( kd.NE.0 ) THEN
195*
196* B is bidiagonal.
197*
198 IF( kd.NE.0 .AND. m.GE.n ) THEN
199*
200* B is upper bidiagonal and M >= N.
201*
202 DO 20 j = 1, n
203 CALL ccopy( m, a( 1, j ), 1, work, 1 )
204 DO 10 i = 1, n - 1
205 work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
206 10 CONTINUE
207 work( m+n ) = d( n )*pt( n, j )
208 CALL cgemv( 'No transpose', m, n, -cmplx( one ), q, ldq,
209 $ work( m+1 ), 1, cmplx( one ), work, 1 )
210 resid = max( resid, scasum( m, work, 1 ) )
211 20 CONTINUE
212 ELSE IF( kd.LT.0 ) THEN
213*
214* B is upper bidiagonal and M < N.
215*
216 DO 40 j = 1, n
217 CALL ccopy( m, a( 1, j ), 1, work, 1 )
218 DO 30 i = 1, m - 1
219 work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
220 30 CONTINUE
221 work( m+m ) = d( m )*pt( m, j )
222 CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
223 $ work( m+1 ), 1, cmplx( one ), work, 1 )
224 resid = max( resid, scasum( m, work, 1 ) )
225 40 CONTINUE
226 ELSE
227*
228* B is lower bidiagonal.
229*
230 DO 60 j = 1, n
231 CALL ccopy( m, a( 1, j ), 1, work, 1 )
232 work( m+1 ) = d( 1 )*pt( 1, j )
233 DO 50 i = 2, m
234 work( m+i ) = e( i-1 )*pt( i-1, j ) +
235 $ d( i )*pt( i, j )
236 50 CONTINUE
237 CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
238 $ work( m+1 ), 1, cmplx( one ), work, 1 )
239 resid = max( resid, scasum( m, work, 1 ) )
240 60 CONTINUE
241 END IF
242 ELSE
243*
244* B is diagonal.
245*
246 IF( m.GE.n ) THEN
247 DO 80 j = 1, n
248 CALL ccopy( m, a( 1, j ), 1, work, 1 )
249 DO 70 i = 1, n
250 work( m+i ) = d( i )*pt( i, j )
251 70 CONTINUE
252 CALL cgemv( 'No transpose', m, n, -cmplx( one ), q, ldq,
253 $ work( m+1 ), 1, cmplx( one ), work, 1 )
254 resid = max( resid, scasum( m, work, 1 ) )
255 80 CONTINUE
256 ELSE
257 DO 100 j = 1, n
258 CALL ccopy( m, a( 1, j ), 1, work, 1 )
259 DO 90 i = 1, m
260 work( m+i ) = d( i )*pt( i, j )
261 90 CONTINUE
262 CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
263 $ work( m+1 ), 1, cmplx( one ), work, 1 )
264 resid = max( resid, scasum( m, work, 1 ) )
265 100 CONTINUE
266 END IF
267 END IF
268*
269* Compute norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
270*
271 anorm = clange( '1', m, n, a, lda, rwork )
272 eps = slamch( 'Precision' )
273*
274 IF( anorm.LE.zero ) THEN
275 IF( resid.NE.zero )
276 $ resid = one / eps
277 ELSE
278 IF( anorm.GE.resid ) THEN
279 resid = ( resid / anorm ) / ( real( n )*eps )
280 ELSE
281 IF( anorm.LT.one ) THEN
282 resid = ( min( resid, real( n )*anorm ) / anorm ) /
283 $ ( real( n )*eps )
284 ELSE
285 resid = min( resid / anorm, real( n ) ) /
286 $ ( real( n )*eps )
287 END IF
288 END IF
289 END IF
290*
291 RETURN
292*
293* End of CBDT01
294*
295 END
subroutine cbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
CBDT01
Definition cbdt01.f:147
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160