zgebd2.f

Go to the documentation of this file.

*> \brief \b ZGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download ZGEBD2 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebd2.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebd2.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebd2.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO )
* 
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, M, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   D( * ), E( * )
*       COMPLEX*16         A( LDA, * ), TAUP( * ), TAUQ( * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZGEBD2 reduces a complex general m by n matrix A to upper or lower
*> real bidiagonal form B by a unitary transformation: Q**H * A * P = B.
*>
*> If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows in the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns in the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA,N)
*>          On entry, the m by n general matrix to be reduced.
*>          On exit,
*>          if m >= n, the diagonal and the first superdiagonal are
*>            overwritten with the upper bidiagonal matrix B; the
*>            elements below the diagonal, with the array TAUQ, represent
*>            the unitary matrix Q as a product of elementary
*>            reflectors, and the elements above the first superdiagonal,
*>            with the array TAUP, represent the unitary matrix P as
*>            a product of elementary reflectors;
*>          if m < n, the diagonal and the first subdiagonal are
*>            overwritten with the lower bidiagonal matrix B; the
*>            elements below the first subdiagonal, with the array TAUQ,
*>            represent the unitary matrix Q as a product of
*>            elementary reflectors, and the elements above the diagonal,
*>            with the array TAUP, represent the unitary matrix P as
*>            a product of elementary reflectors.
*>          See Further Details.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (min(M,N))
*>          The diagonal elements of the bidiagonal matrix B:
*>          D(i) = A(i,i).
*> \endverbatim
*>
*> \param[out] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
*>          The off-diagonal elements of the bidiagonal matrix B:
*>          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
*>          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
*> \endverbatim
*>
*> \param[out] TAUQ
*> \verbatim
*>          TAUQ is COMPLEX*16 array dimension (min(M,N))
*>          The scalar factors of the elementary reflectors which
*>          represent the unitary matrix Q. See Further Details.
*> \endverbatim
*>
*> \param[out] TAUP
*> \verbatim
*>          TAUP is COMPLEX*16 array, dimension (min(M,N))
*>          The scalar factors of the elementary reflectors which
*>          represent the unitary matrix P. See Further Details.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension (max(M,N))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup complex16GEcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  The matrices Q and P are represented as products of elementary
*>  reflectors:
*>
*>  If m >= n,
*>
*>     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
*>
*>  Each H(i) and G(i) has the form:
*>
*>     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H
*>
*>  where tauq and taup are complex scalars, and v and u are complex
*>  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
*>  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
*>  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
*>
*>  If m < n,
*>
*>     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
*>
*>  Each H(i) and G(i) has the form:
*>
*>     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H
*>
*>  where tauq and taup are complex scalars, v and u are complex vectors;
*>  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
*>  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
*>  tauq is stored in TAUQ(i) and taup in TAUP(i).
*>
*>  The contents of A on exit are illustrated by the following examples:
*>
*>  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
*>
*>    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
*>    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
*>    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
*>    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
*>    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
*>    (  v1  v2  v3  v4  v5 )
*>
*>  where d and e denote diagonal and off-diagonal elements of B, vi
*>  denotes an element of the vector defining H(i), and ui an element of
*>  the vector defining G(i).
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE zgebd2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO )
*
*  -- LAPACK computational routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            info, lda, m, n
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   d( * ), e( * )
      COMPLEX*16         a( lda, * ), taup( * ), tauq( * ), work( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         zero, one
      parameter( zero = ( 0.0d+0, 0.0d+0 ),
     $                   one = ( 1.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            i
      COMPLEX*16         alpha
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zlacgv, zlarf, zlarfg
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dconjg, max, min
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      info = 0
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -4
      END IF
      IF( info.LT.0 ) THEN
         CALL xerbla( 'ZGEBD2', -info )
         return
      END IF
*
      IF( m.GE.n ) THEN
*
*        Reduce to upper bidiagonal form
*
         DO 10 i = 1, n
*
*           Generate elementary reflector H(i) to annihilate A(i+1:m,i)
*
            alpha = a( i, i )
            CALL zlarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1,
     $                   tauq( i ) )
            d( i ) = alpha
            a( i, i ) = one
*
*           Apply H(i)**H to A(i:m,i+1:n) from the left
*
            IF( i.LT.n )
     $         CALL zlarf( 'Left', m-i+1, n-i, a( i, i ), 1,
     $                     dconjg( tauq( i ) ), a( i, i+1 ), lda, work )
            a( i, i ) = d( i )
*
            IF( i.LT.n ) THEN
*
*              Generate elementary reflector G(i) to annihilate
*              A(i,i+2:n)
*
               CALL zlacgv( n-i, a( i, i+1 ), lda )
               alpha = a( i, i+1 )
               CALL zlarfg( n-i, alpha, a( i, min( i+2, n ) ), lda,
     $                      taup( i ) )
               e( i ) = alpha
               a( i, i+1 ) = one
*
*              Apply G(i) to A(i+1:m,i+1:n) from the right
*
               CALL zlarf( 'Right', m-i, n-i, a( i, i+1 ), lda,
     $                     taup( i ), a( i+1, i+1 ), lda, work )
               CALL zlacgv( n-i, a( i, i+1 ), lda )
               a( i, i+1 ) = e( i )
            ELSE
               taup( i ) = zero
            END IF
   10    continue
      ELSE
*
*        Reduce to lower bidiagonal form
*
         DO 20 i = 1, m
*
*           Generate elementary reflector G(i) to annihilate A(i,i+1:n)
*
            CALL zlacgv( n-i+1, a( i, i ), lda )
            alpha = a( i, i )
            CALL zlarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,
     $                   taup( i ) )
            d( i ) = alpha
            a( i, i ) = one
*
*           Apply G(i) to A(i+1:m,i:n) from the right
*
            IF( i.LT.m )
     $         CALL zlarf( 'Right', m-i, n-i+1, a( i, i ), lda,
     $                     taup( i ), a( i+1, i ), lda, work )
            CALL zlacgv( n-i+1, a( i, i ), lda )
            a( i, i ) = d( i )
*
            IF( i.LT.m ) THEN
*
*              Generate elementary reflector H(i) to annihilate
*              A(i+2:m,i)
*
               alpha = a( i+1, i )
               CALL zlarfg( m-i, alpha, a( min( i+2, m ), i ), 1,
     $                      tauq( i ) )
               e( i ) = alpha
               a( i+1, i ) = one
*
*              Apply H(i)**H to A(i+1:m,i+1:n) from the left
*
               CALL zlarf( 'Left', m-i, n-i, a( i+1, i ), 1,
     $                     dconjg( tauq( i ) ), a( i+1, i+1 ), lda,
     $                     work )
               a( i+1, i ) = e( i )
            ELSE
               tauq( i ) = zero
            END IF
   20    continue
      END IF
      return
*
*     End of ZGEBD2
*
      END