zgbbrd.f

Go to the documentation of this file.

*> \brief \b ZGBBRD
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download ZGBBRD + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbbrd.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbbrd.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbbrd.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
*                          LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          VECT
*       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
*       COMPLEX*16         AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
*      $                   Q( LDQ, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZGBBRD reduces a complex general m-by-n band matrix A to real upper
*> bidiagonal form B by a unitary transformation: Q**H * A * P = B.
*>
*> The routine computes B, and optionally forms Q or P**H, or computes
*> Q**H*C for a given matrix C.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] VECT
*> \verbatim
*>          VECT is CHARACTER*1
*>          Specifies whether or not the matrices Q and P**H are to be
*>          formed.
*>          = 'N': do not form Q or P**H;
*>          = 'Q': form Q only;
*>          = 'P': form P**H only;
*>          = 'B': form both.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] NCC
*> \verbatim
*>          NCC is INTEGER
*>          The number of columns of the matrix C.  NCC >= 0.
*> \endverbatim
*>
*> \param[in] KL
*> \verbatim
*>          KL is INTEGER
*>          The number of subdiagonals of the matrix A. KL >= 0.
*> \endverbatim
*>
*> \param[in] KU
*> \verbatim
*>          KU is INTEGER
*>          The number of superdiagonals of the matrix A. KU >= 0.
*> \endverbatim
*>
*> \param[in,out] AB
*> \verbatim
*>          AB is COMPLEX*16 array, dimension (LDAB,N)
*>          On entry, the m-by-n band matrix A, stored in rows 1 to
*>          KL+KU+1. The j-th column of A is stored in the j-th column of
*>          the array AB as follows:
*>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
*>          On exit, A is overwritten by values generated during the
*>          reduction.
*> \endverbatim
*>
*> \param[in] LDAB
*> \verbatim
*>          LDAB is INTEGER
*>          The leading dimension of the array A. LDAB >= KL+KU+1.
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (min(M,N))
*>          The diagonal elements of the bidiagonal matrix B.
*> \endverbatim
*>
*> \param[out] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
*>          The superdiagonal elements of the bidiagonal matrix B.
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is COMPLEX*16 array, dimension (LDQ,M)
*>          If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
*>          If VECT = 'N' or 'P', the array Q is not referenced.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
*> \endverbatim
*>
*> \param[out] PT
*> \verbatim
*>          PT is COMPLEX*16 array, dimension (LDPT,N)
*>          If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
*>          If VECT = 'N' or 'Q', the array PT is not referenced.
*> \endverbatim
*>
*> \param[in] LDPT
*> \verbatim
*>          LDPT is INTEGER
*>          The leading dimension of the array PT.
*>          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is COMPLEX*16 array, dimension (LDC,NCC)
*>          On entry, an m-by-ncc matrix C.
*>          On exit, C is overwritten by Q**H*C.
*>          C is not referenced if NCC = 0.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C.
*>          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension (max(M,N))
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION 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.
*
*> \ingroup gbbrd
*
*  =====================================================================
      SUBROUTINE zgbbrd( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
     $                   LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          VECT
      INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
      COMPLEX*16         AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
     $                   q( ldq, * ), work( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d+0 )
      COMPLEX*16         CZERO, CONE
      parameter( czero = ( 0.0d+0, 0.0d+0 ),
     $                   cone = ( 1.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            WANTB, WANTC, WANTPT, WANTQ
      INTEGER            I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
     $                   kun, l, minmn, ml, ml0, mu, mu0, nr, nrt
      DOUBLE PRECISION   ABST, RC
      COMPLEX*16         RA, RB, RS, T
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zlargv, zlartg, zlartv, zlaset,
     $                   zrot,
     $                   zscal
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dconjg, max, min
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      wantb = lsame( vect, 'B' )
      wantq = lsame( vect, 'Q' ) .OR. wantb
      wantpt = lsame( vect, 'P' ) .OR. wantb
      wantc = ncc.GT.0
      klu1 = kl + ku + 1
      info = 0
      IF( .NOT.wantq .AND.
     $    .NOT.wantpt .AND.
     $    .NOT.lsame( vect, 'N' ) )
     $     THEN
         info = -1
      ELSE IF( m.LT.0 ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( ncc.LT.0 ) THEN
         info = -4
      ELSE IF( kl.LT.0 ) THEN
         info = -5
      ELSE IF( ku.LT.0 ) THEN
         info = -6
      ELSE IF( ldab.LT.klu1 ) THEN
         info = -8
      ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
         info = -12
      ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
         info = -14
      ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
         info = -16
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZGBBRD', -info )
         RETURN
      END IF
*
*     Initialize Q and P**H to the unit matrix, if needed
*
      IF( wantq )
     $   CALL zlaset( 'Full', m, m, czero, cone, q, ldq )
      IF( wantpt )
     $   CALL zlaset( 'Full', n, n, czero, cone, pt, ldpt )
*
*     Quick return if possible.
*
      IF( m.EQ.0 .OR. n.EQ.0 )
     $   RETURN
*
      minmn = min( m, n )
*
      IF( kl+ku.GT.1 ) THEN
*
*        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
*        first to lower bidiagonal form and then transform to upper
*        bidiagonal
*
         IF( ku.GT.0 ) THEN
            ml0 = 1
            mu0 = 2
         ELSE
            ml0 = 2
            mu0 = 1
         END IF
*
*        Wherever possible, plane rotations are generated and applied in
*        vector operations of length NR over the index set J1:J2:KLU1.
*
*        The complex sines of the plane rotations are stored in WORK,
*        and the real cosines in RWORK.
*
         klm = min( m-1, kl )
         kun = min( n-1, ku )
         kb = klm + kun
         kb1 = kb + 1
         inca = kb1*ldab
         nr = 0
         j1 = klm + 2
         j2 = 1 - kun
*
         DO 90 i = 1, minmn
*
*           Reduce i-th column and i-th row of matrix to bidiagonal form
*
            ml = klm + 1
            mu = kun + 1
            DO 80 kk = 1, kb
               j1 = j1 + kb
               j2 = j2 + kb
*
*              generate plane rotations to annihilate nonzero elements
*              which have been created below the band
*
               IF( nr.GT.0 )
     $            CALL zlargv( nr, ab( klu1, j1-klm-1 ), inca,
     $                         work( j1 ), kb1, rwork( j1 ), kb1 )
*
*              apply plane rotations from the left
*
               DO 10 l = 1, kb
                  IF( j2-klm+l-1.GT.n ) THEN
                     nrt = nr - 1
                  ELSE
                     nrt = nr
                  END IF
                  IF( nrt.GT.0 )
     $               CALL zlartv( nrt, ab( klu1-l, j1-klm+l-1 ),
     $                            inca,
     $                            ab( klu1-l+1, j1-klm+l-1 ), inca,
     $                            rwork( j1 ), work( j1 ), kb1 )
   10          CONTINUE
*
               IF( ml.GT.ml0 ) THEN
                  IF( ml.LE.m-i+1 ) THEN
*
*                    generate plane rotation to annihilate a(i+ml-1,i)
*                    within the band, and apply rotation from the left
*
                     CALL zlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
     $                            rwork( i+ml-1 ), work( i+ml-1 ), ra )
                     ab( ku+ml-1, i ) = ra
                     IF( i.LT.n )
     $                  CALL zrot( min( ku+ml-2, n-i ),
     $                             ab( ku+ml-2, i+1 ), ldab-1,
     $                             ab( ku+ml-1, i+1 ), ldab-1,
     $                             rwork( i+ml-1 ), work( i+ml-1 ) )
                  END IF
                  nr = nr + 1
                  j1 = j1 - kb1
               END IF
*
               IF( wantq ) THEN
*
*                 accumulate product of plane rotations in Q
*
                  DO 20 j = j1, j2, kb1
                     CALL zrot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
     $                          rwork( j ), dconjg( work( j ) ) )
   20             CONTINUE
               END IF
*
               IF( wantc ) THEN
*
*                 apply plane rotations to C
*
                  DO 30 j = j1, j2, kb1
                     CALL zrot( ncc, c( j-1, 1 ), ldc, c( j, 1 ),
     $                          ldc,
     $                          rwork( j ), work( j ) )
   30             CONTINUE
               END IF
*
               IF( j2+kun.GT.n ) THEN
*
*                 adjust J2 to keep within the bounds of the matrix
*
                  nr = nr - 1
                  j2 = j2 - kb1
               END IF
*
               DO 40 j = j1, j2, kb1
*
*                 create nonzero element a(j-1,j+ku) above the band
*                 and store it in WORK(n+1:2*n)
*
                  work( j+kun ) = work( j )*ab( 1, j+kun )
                  ab( 1, j+kun ) = rwork( j )*ab( 1, j+kun )
   40          CONTINUE
*
*              generate plane rotations to annihilate nonzero elements
*              which have been generated above the band
*
               IF( nr.GT.0 )
     $            CALL zlargv( nr, ab( 1, j1+kun-1 ), inca,
     $                         work( j1+kun ), kb1, rwork( j1+kun ),
     $                         kb1 )
*
*              apply plane rotations from the right
*
               DO 50 l = 1, kb
                  IF( j2+l-1.GT.m ) THEN
                     nrt = nr - 1
                  ELSE
                     nrt = nr
                  END IF
                  IF( nrt.GT.0 )
     $               CALL zlartv( nrt, ab( l+1, j1+kun-1 ), inca,
     $                            ab( l, j1+kun ), inca,
     $                            rwork( j1+kun ), work( j1+kun ), kb1 )
   50          CONTINUE
*
               IF( ml.EQ.ml0 .AND. mu.GT.mu0 ) THEN
                  IF( mu.LE.n-i+1 ) THEN
*
*                    generate plane rotation to annihilate a(i,i+mu-1)
*                    within the band, and apply rotation from the right
*
                     CALL zlartg( ab( ku-mu+3, i+mu-2 ),
     $                            ab( ku-mu+2, i+mu-1 ),
     $                            rwork( i+mu-1 ), work( i+mu-1 ), ra )
                     ab( ku-mu+3, i+mu-2 ) = ra
                     CALL zrot( min( kl+mu-2, m-i ),
     $                          ab( ku-mu+4, i+mu-2 ), 1,
     $                          ab( ku-mu+3, i+mu-1 ), 1,
     $                          rwork( i+mu-1 ), work( i+mu-1 ) )
                  END IF
                  nr = nr + 1
                  j1 = j1 - kb1
               END IF
*
               IF( wantpt ) THEN
*
*                 accumulate product of plane rotations in P**H
*
                  DO 60 j = j1, j2, kb1
                     CALL zrot( n, pt( j+kun-1, 1 ), ldpt,
     $                          pt( j+kun, 1 ), ldpt, rwork( j+kun ),
     $                          dconjg( work( j+kun ) ) )
   60             CONTINUE
               END IF
*
               IF( j2+kb.GT.m ) THEN
*
*                 adjust J2 to keep within the bounds of the matrix
*
                  nr = nr - 1
                  j2 = j2 - kb1
               END IF
*
               DO 70 j = j1, j2, kb1
*
*                 create nonzero element a(j+kl+ku,j+ku-1) below the
*                 band and store it in WORK(1:n)
*
                  work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                  ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun )
   70          CONTINUE
*
               IF( ml.GT.ml0 ) THEN
                  ml = ml - 1
               ELSE
                  mu = mu - 1
               END IF
   80       CONTINUE
   90    CONTINUE
      END IF
*
      IF( ku.EQ.0 .AND. kl.GT.0 ) THEN
*
*        A has been reduced to complex lower bidiagonal form
*
*        Transform lower bidiagonal form to upper bidiagonal by applying
*        plane rotations from the left, overwriting superdiagonal
*        elements on subdiagonal elements
*
         DO 100 i = 1, min( m-1, n )
            CALL zlartg( ab( 1, i ), ab( 2, i ), rc, rs, ra )
            ab( 1, i ) = ra
            IF( i.LT.n ) THEN
               ab( 2, i ) = rs*ab( 1, i+1 )
               ab( 1, i+1 ) = rc*ab( 1, i+1 )
            END IF
            IF( wantq )
     $         CALL zrot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc,
     $                    dconjg( rs ) )
            IF( wantc )
     $         CALL zrot( ncc, c( i, 1 ), ldc, c( i+1, 1 ), ldc, rc,
     $                    rs )
  100    CONTINUE
      ELSE
*
*        A has been reduced to complex upper bidiagonal form or is
*        diagonal
*
         IF( ku.GT.0 .AND. m.LT.n ) THEN
*
*           Annihilate a(m,m+1) by applying plane rotations from the
*           right
*
            rb = ab( ku, m+1 )
            DO 110 i = m, 1, -1
               CALL zlartg( ab( ku+1, i ), rb, rc, rs, ra )
               ab( ku+1, i ) = ra
               IF( i.GT.1 ) THEN
                  rb = -dconjg( rs )*ab( ku, i )
                  ab( ku, i ) = rc*ab( ku, i )
               END IF
               IF( wantpt )
     $            CALL zrot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
     $                       rc, dconjg( rs ) )
  110       CONTINUE
         END IF
      END IF
*
*     Make diagonal and superdiagonal elements real, storing them in D
*     and E
*
      t = ab( ku+1, 1 )
      DO 120 i = 1, minmn
         abst = abs( t )
         d( i ) = abst
         IF( abst.NE.zero ) THEN
            t = t / abst
         ELSE
            t = cone
         END IF
         IF( wantq )
     $      CALL zscal( m, t, q( 1, i ), 1 )
         IF( wantc )
     $      CALL zscal( ncc, dconjg( t ), c( i, 1 ), ldc )
         IF( i.LT.minmn ) THEN
            IF( ku.EQ.0 .AND. kl.EQ.0 ) THEN
               e( i ) = zero
               t = ab( 1, i+1 )
            ELSE
               IF( ku.EQ.0 ) THEN
                  t = ab( 2, i )*dconjg( t )
               ELSE
                  t = ab( ku, i+1 )*dconjg( t )
               END IF
               abst = abs( t )
               e( i ) = abst
               IF( abst.NE.zero ) THEN
                  t = t / abst
               ELSE
                  t = cone
               END IF
               IF( wantpt )
     $            CALL zscal( n, t, pt( i+1, 1 ), ldpt )
               t = ab( ku+1, i+1 )*dconjg( t )
            END IF
         END IF
  120 CONTINUE
      RETURN
*
*     End of ZGBBRD
*
      END