LAPACK 3.12.0
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:115
Here is the call graph for this function:
Here is the caller graph for this function: