146 SUBROUTINE zbdt01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
155 INTEGER kd, lda, ldpt, ldq, m, n
156 DOUBLE PRECISION resid
159 DOUBLE PRECISION d( * ), e( * ), rwork( * )
160 COMPLEX*16 a( lda, * ), pt( ldpt, * ), q( ldq, * ),
167 DOUBLE PRECISION zero, one
168 parameter( zero = 0.0d+0, one = 1.0d+0 )
172 DOUBLE PRECISION anorm, eps
182 INTRINSIC dble, dcmplx, max, min
188 IF( m.LE.0 .OR. n.LE.0 )
THEN
200 IF( kd.NE.0 .AND. m.GE.n )
THEN
205 CALL
zcopy( m, a( 1, j ), 1, work, 1 )
207 work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
209 work( m+n ) = d( n )*pt( n, j )
210 CALL
zgemv(
'No transpose', m, n, -dcmplx( one ), q, ldq,
211 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
212 resid = max( resid,
dzasum( m, work, 1 ) )
214 ELSE IF( kd.LT.0 )
THEN
219 CALL
zcopy( m, a( 1, j ), 1, work, 1 )
221 work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
223 work( m+m ) = d( m )*pt( m, j )
224 CALL
zgemv(
'No transpose', m, m, -dcmplx( one ), q, ldq,
225 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
226 resid = max( resid,
dzasum( m, work, 1 ) )
233 CALL
zcopy( m, a( 1, j ), 1, work, 1 )
234 work( m+1 ) = d( 1 )*pt( 1, j )
236 work( m+i ) = e( i-1 )*pt( i-1, j ) +
239 CALL
zgemv(
'No transpose', m, m, -dcmplx( one ), q, ldq,
240 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
241 resid = max( resid,
dzasum( m, work, 1 ) )
250 CALL
zcopy( m, a( 1, j ), 1, work, 1 )
252 work( m+i ) = d( i )*pt( i, j )
254 CALL
zgemv(
'No transpose', m, n, -dcmplx( one ), q, ldq,
255 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
256 resid = max( resid,
dzasum( m, work, 1 ) )
260 CALL
zcopy( m, a( 1, j ), 1, work, 1 )
262 work( m+i ) = d( i )*pt( i, j )
264 CALL
zgemv(
'No transpose', m, m, -dcmplx( one ), q, ldq,
265 $ work( m+1 ), 1, dcmplx( one ), work, 1 )
266 resid = max( resid,
dzasum( m, work, 1 ) )
273 anorm =
zlange(
'1', m, n, a, lda, rwork )
274 eps =
dlamch(
'Precision' )
276 IF( anorm.LE.zero )
THEN
280 IF( anorm.GE.resid )
THEN
281 resid = ( resid / anorm ) / ( dble( n )*eps )
283 IF( anorm.LT.one )
THEN
284 resid = ( min( resid, dble( n )*anorm ) / anorm ) /
287 resid = min( resid / anorm, dble( n ) ) /