LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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 )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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*
scasum
real function scasum(n, cx, incx)
SCASUM
Definition scasum.f:72
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
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
Here is the call graph for this function:
Here is the caller graph for this function: