LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zbdt01()

subroutine zbdt01 ( integer m,
integer n,
integer kd,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
complex*16, dimension( ldpt, * ) pt,
integer ldpt,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
double precision resid )

ZBDT01

Purpose:
!>
!> ZBDT01 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*16 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*16 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 DOUBLE PRECISION array, dimension (min(M,N))
!>          The diagonal elements of the bidiagonal matrix B.
!>
[in]E
!>          E is DOUBLE PRECISION 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*16 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*16 array, dimension (M+N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (M)
!>
[out]RESID
!>          RESID is DOUBLE PRECISION
!>          The test ratio:
!>          norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 145 of file zbdt01.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 DOUBLE PRECISION RESID
155* ..
156* .. Array Arguments ..
157 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
158 COMPLEX*16 A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
159 $ WORK( * )
160* ..
161*
162* =====================================================================
163*
164* .. Parameters ..
165 DOUBLE PRECISION ZERO, ONE
166 parameter( zero = 0.0d+0, one = 1.0d+0 )
167* ..
168* .. Local Scalars ..
169 INTEGER I, J
170 DOUBLE PRECISION ANORM, EPS
171* ..
172* .. External Functions ..
173 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
174 EXTERNAL dlamch, dzasum, zlange
175* ..
176* .. External Subroutines ..
177 EXTERNAL zcopy, zgemv
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC dble, dcmplx, max, min
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 zcopy( 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 zgemv( 'No transpose', m, n, -dcmplx( one ), q, ldq,
209 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
210 resid = max( resid, dzasum( 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 zcopy( 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 zgemv( 'No transpose', m, m, -dcmplx( one ), q, ldq,
223 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
224 resid = max( resid, dzasum( m, work, 1 ) )
225 40 CONTINUE
226 ELSE
227*
228* B is lower bidiagonal.
229*
230 DO 60 j = 1, n
231 CALL zcopy( 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 zgemv( 'No transpose', m, m, -dcmplx( one ), q, ldq,
238 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
239 resid = max( resid, dzasum( 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 zcopy( 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 zgemv( 'No transpose', m, n, -dcmplx( one ), q, ldq,
253 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
254 resid = max( resid, dzasum( m, work, 1 ) )
255 80 CONTINUE
256 ELSE
257 DO 100 j = 1, n
258 CALL zcopy( 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 zgemv( 'No transpose', m, m, -dcmplx( one ), q, ldq,
263 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
264 resid = max( resid, dzasum( 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 = zlange( '1', m, n, a, lda, rwork )
272 eps = dlamch( '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 ) / ( dble( n )*eps )
280 ELSE
281 IF( anorm.LT.one ) THEN
282 resid = ( min( resid, dble( n )*anorm ) / anorm ) /
283 $ ( dble( n )*eps )
284 ELSE
285 resid = min( resid / anorm, dble( n ) ) /
286 $ ( dble( n )*eps )
287 END IF
288 END IF
289 END IF
290*
291 RETURN
292*
293* End of ZBDT01
294*
dzasum
double precision function dzasum(n, zx, incx)
DZASUM
Definition dzasum.f:72
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:113
Here is the call graph for this function:
Here is the caller graph for this function: