sbdt01()

subroutine sbdt01	(	integer	m,
		integer	n,
		integer	kd,
		real, dimension( lda, * )	a,
		integer	lda,
		real, dimension( ldq, * )	q,
		integer	ldq,
		real, dimension( * )	d,
		real, dimension( * )	e,
		real, dimension( ldpt, * )	pt,
		integer	ldpt,
		real, dimension( * )	work,
		real	resid )

SBDT01

Purpose:

!>
!> SBDT01 reconstructs a general matrix A from its bidiagonal form
!>    A = Q * B * P**T
!> where Q (m by min(m,n)) and P**T (min(m,n) by n) are orthogonal
!> matrices and B is bidiagonal.
!>
!> The test ratio to test the reduction is
!>    RESID = norm(A - Q * B * P**T) / ( 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**T. !>
[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 REAL 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 REAL array, dimension (LDQ,N) !> The m by min(m,n) orthogonal matrix Q in the reduction !> A = Q * B * P**T. !>
[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 REAL array, dimension (LDPT,N) !> The min(m,n) by n orthogonal matrix P**T in the reduction !> A = Q * B * P**T. !>
[in]	LDPT	!> LDPT is INTEGER !> The leading dimension of the array PT. !> LDPT >= max(1,min(M,N)). !>
[out]	WORK	!> WORK is REAL array, dimension (M+N) !>
[out]	RESID	!> RESID is REAL !> The test ratio: !> norm(A - Q * B * P**T) / ( n * norm(A) * EPS ) !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 139 of file sbdt01.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               A( LDA, * ), D( * ), E( * ), 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               SASUM, SLAMCH, SLANGE
      EXTERNAL           sasum, slamch, slange
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgemv
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          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**T 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 scopy( 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 sgemv( 'No transpose', m, n, -one, q, ldq,
     $                     work( m+1 ), 1, one, work, 1 )
               resid = max( resid, sasum( m, work, 1 ) )
   20       CONTINUE
         ELSE IF( kd.LT.0 ) THEN
*
*           B is upper bidiagonal and M < N.
*
            DO 40 j = 1, n
               CALL scopy( 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 sgemv( 'No transpose', m, m, -one, q, ldq,
     $                     work( m+1 ), 1, one, work, 1 )
               resid = max( resid, sasum( m, work, 1 ) )
   40       CONTINUE
         ELSE
*
*           B is lower bidiagonal.
*
            DO 60 j = 1, n
               CALL scopy( 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 sgemv( 'No transpose', m, m, -one, q, ldq,
     $                     work( m+1 ), 1, one, work, 1 )
               resid = max( resid, sasum( m, work, 1 ) )
   60       CONTINUE
         END IF
      ELSE
*
*        B is diagonal.
*
         IF( m.GE.n ) THEN
            DO 80 j = 1, n
               CALL scopy( m, a( 1, j ), 1, work, 1 )
               DO 70 i = 1, n
                  work( m+i ) = d( i )*pt( i, j )
   70          CONTINUE
               CALL sgemv( 'No transpose', m, n, -one, q, ldq,
     $                     work( m+1 ), 1, one, work, 1 )
               resid = max( resid, sasum( m, work, 1 ) )
   80       CONTINUE
         ELSE
            DO 100 j = 1, n
               CALL scopy( m, a( 1, j ), 1, work, 1 )
               DO 90 i = 1, m
                  work( m+i ) = d( i )*pt( i, j )
   90          CONTINUE
               CALL sgemv( 'No transpose', m, m, -one, q, ldq,
     $                     work( m+1 ), 1, one, work, 1 )
               resid = max( resid, sasum( m, work, 1 ) )
  100       CONTINUE
         END IF
      END IF
*
*     Compute norm(A - Q * B * P**T) / ( n * norm(A) * EPS )
*
      anorm = slange( '1', m, n, a, lda, work )
      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 SBDT01
*

Here is the call graph for this function:

Here is the caller graph for this function: