LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ cbdt01()

 subroutine cbdt01 ( integer M, integer N, integer KD, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldq, * ) Q, integer LDQ, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( ldpt, * ) PT, integer LDPT, complex, dimension( * ) WORK, real, dimension( * ) RWORK, real RESID )

CBDT01

Purpose:
``` CBDT01 reconstructs a general matrix A from its bidiagonal form
A = Q * B * P**H
where Q (m by min(m,n)) and P**H (min(m,n) by n) are unitary
matrices and B is bidiagonal.

The test ratio to test the reduction is
RESID = norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
where EPS is the machine precision.```
Parameters
 [in] M ``` M is INTEGER The number of rows of the matrices A and Q.``` [in] N ``` N is INTEGER The number of columns of the matrices A and P**H.``` [in] KD ``` KD is INTEGER If KD = 0, B is diagonal and the array E is not referenced. If KD = 1, the reduction was performed by xGEBRD; B is upper bidiagonal if M >= N, and lower bidiagonal if M < N. If KD = -1, the reduction was performed by xGBBRD; B is always upper bidiagonal.``` [in] A ``` A is COMPLEX array, dimension (LDA,N) The m by n matrix A.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).``` [in] Q ``` Q is COMPLEX array, dimension (LDQ,N) The m by min(m,n) unitary matrix Q in the reduction A = Q * B * P**H.``` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,M).``` [in] D ``` D is REAL array, dimension (min(M,N)) The diagonal elements of the bidiagonal matrix B.``` [in] E ``` E is REAL array, dimension (min(M,N)-1) The superdiagonal elements of the bidiagonal matrix B if m >= n, or the subdiagonal elements of B if m < n.``` [in] PT ``` PT is COMPLEX array, dimension (LDPT,N) The min(m,n) by n unitary matrix P**H in the reduction A = Q * B * P**H.``` [in] LDPT ``` LDPT is INTEGER The leading dimension of the array PT. LDPT >= max(1,min(M,N)).``` [out] WORK ` WORK is COMPLEX array, dimension (M+N)` [out] RWORK ` RWORK is REAL array, dimension (M)` [out] RESID ``` RESID is REAL The test ratio: norm(A - Q * B * P**H) / ( n * norm(A) * EPS )```

Definition at line 145 of file cbdt01.f.

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*
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:158
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:115
real function scasum(N, CX, INCX)
SCASUM
Definition: scasum.f:72
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: