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.

*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            KD, LDA, LDPT, LDQ, M, N
      REAL               RESID
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * ), RWORK( * )
      COMPLEX            A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
     $                   WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               ANORM, EPS
*     ..
*     .. External Functions ..
      REAL               CLANGE, SCASUM, SLAMCH
      EXTERNAL           clange, scasum, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, cgemv
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cmplx, max, min, real
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( m.LE.0 .OR. n.LE.0 ) THEN
         resid = zero
         RETURN
      END IF
*
*     Compute A - Q * B * P**H one column at a time.
*
      resid = zero
      IF( kd.NE.0 ) THEN
*
*        B is bidiagonal.
*
         IF( kd.NE.0 .AND. m.GE.n ) THEN
*
*           B is upper bidiagonal and M >= N.
*
            DO 20 j = 1, n
               CALL ccopy( m, a( 1, j ), 1, work, 1 )
               DO 10 i = 1, n - 1
                  work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
   10          CONTINUE
               work( m+n ) = d( n )*pt( n, j )
               CALL cgemv( 'No transpose', m, n, -cmplx( one ), q, ldq,
     $                     work( m+1 ), 1, cmplx( one ), work, 1 )
               resid = max( resid, scasum( m, work, 1 ) )
   20       CONTINUE
         ELSE IF( kd.LT.0 ) THEN
*
*           B is upper bidiagonal and M < N.
*
            DO 40 j = 1, n
               CALL ccopy( m, a( 1, j ), 1, work, 1 )
               DO 30 i = 1, m - 1
                  work( m+i ) = d( i )*pt( i, j ) + e( i )*pt( i+1, j )
   30          CONTINUE
               work( m+m ) = d( m )*pt( m, j )
               CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
     $                     work( m+1 ), 1, cmplx( one ), work, 1 )
               resid = max( resid, scasum( m, work, 1 ) )
   40       CONTINUE
         ELSE
*
*           B is lower bidiagonal.
*
            DO 60 j = 1, n
               CALL ccopy( m, a( 1, j ), 1, work, 1 )
               work( m+1 ) = d( 1 )*pt( 1, j )
               DO 50 i = 2, m
                  work( m+i ) = e( i-1 )*pt( i-1, j ) +
     $                          d( i )*pt( i, j )
   50          CONTINUE
               CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
     $                     work( m+1 ), 1, cmplx( one ), work, 1 )
               resid = max( resid, scasum( m, work, 1 ) )
   60       CONTINUE
         END IF
      ELSE
*
*        B is diagonal.
*
         IF( m.GE.n ) THEN
            DO 80 j = 1, n
               CALL ccopy( m, a( 1, j ), 1, work, 1 )
               DO 70 i = 1, n
                  work( m+i ) = d( i )*pt( i, j )
   70          CONTINUE
               CALL cgemv( 'No transpose', m, n, -cmplx( one ), q, ldq,
     $                     work( m+1 ), 1, cmplx( one ), work, 1 )
               resid = max( resid, scasum( m, work, 1 ) )
   80       CONTINUE
         ELSE
            DO 100 j = 1, n
               CALL ccopy( m, a( 1, j ), 1, work, 1 )
               DO 90 i = 1, m
                  work( m+i ) = d( i )*pt( i, j )
   90          CONTINUE
               CALL cgemv( 'No transpose', m, m, -cmplx( one ), q, ldq,
     $                     work( m+1 ), 1, cmplx( one ), work, 1 )
               resid = max( resid, scasum( m, work, 1 ) )
  100       CONTINUE
         END IF
      END IF
*
*     Compute norm(A - Q * B * P**H) / ( n * norm(A) * EPS )
*
      anorm = clange( '1', m, n, a, lda, rwork )
      eps = slamch( 'Precision' )
*
      IF( anorm.LE.zero ) THEN
         IF( resid.NE.zero )
     $      resid = one / eps
      ELSE
         IF( anorm.GE.resid ) THEN
            resid = ( resid / anorm ) / ( real( n )*eps )
         ELSE
            IF( anorm.LT.one ) THEN
               resid = ( min( resid, real( n )*anorm ) / anorm ) /
     $                 ( real( n )*eps )
            ELSE
               resid = min( resid / anorm, real( n ) ) /
     $                 ( real( n )*eps )
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of CBDT01
*

Here is the call graph for this function:

Here is the caller graph for this function: