zhetrf_aa_2stage.f

Go to the documentation of this file.

*> \brief \b ZHETRF_AA_2STAGE
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download ZHETRF_AA_2STAGE + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*      SUBROUTINE ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
*                                   IPIV2, WORK, LWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            N, LDA, LTB, LWORK, INFO
*       ..
*       .. Array Arguments ..
*       INTEGER            IPIV( * ), IPIV2( * )
*       COMPLEX*16         A( LDA, * ), TB( * ), WORK( * )
*       ..
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
*> using the Aasen's algorithm.  The form of the factorization is
*>
*>    A = U**H*T*U  or  A = L*T*L**H
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a hermitian band matrix with the
*> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is 
*> LU factorized with partial pivoting).
*>
*> This is the blocked version of the algorithm, calling Level 3 BLAS.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  Upper triangle of A is stored;
*>          = 'L':  Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA,N)
*>          On entry, the hermitian matrix A.  If UPLO = 'U', the leading
*>          N-by-N upper triangular part of A contains the upper
*>          triangular part of the matrix A, and the strictly lower
*>          triangular part of A is not referenced.  If UPLO = 'L', the
*>          leading N-by-N lower triangular part of A contains the lower
*>          triangular part of the matrix A, and the strictly upper
*>          triangular part of A is not referenced.
*>
*>          On exit, L is stored below (or above) the subdiagonal blocks,
*>          when UPLO  is 'L' (or 'U').
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] TB
*> \verbatim
*>          TB is COMPLEX*16 array, dimension (MAX(1,LTB))
*>          On exit, details of the LU factorization of the band matrix.
*> \endverbatim
*>
*> \param[in] LTB
*> \verbatim
*>          LTB is INTEGER
*>          The size of the array TB. LTB >= MAX(1,4*N), internally
*>          used to select NB such that LTB >= (3*NB+1)*N.
*>
*>          If LTB = -1, then a workspace query is assumed; the
*>          routine only calculates the optimal size of LTB, 
*>          returns this value as the first entry of TB, and
*>          no error message related to LTB is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension (N)
*>          On exit, it contains the details of the interchanges, i.e.,
*>          the row and column k of A were interchanged with the
*>          row and column IPIV(k).
*> \endverbatim
*>
*> \param[out] IPIV2
*> \verbatim
*>          IPIV2 is INTEGER array, dimension (N)
*>          On exit, it contains the details of the interchanges, i.e.,
*>          the row and column k of T were interchanged with the
*>          row and column IPIV(k).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 workspace of size (MAX(1,LWORK))
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The size of WORK. LWORK >= MAX(1,N), internally used to
*>          select NB such that LWORK >= N*NB.
*>
*>          If LWORK = -1, then a workspace query is assumed; the
*>          routine only calculates the optimal size of the WORK array,
*>          returns this value as the first entry of the WORK array, and
*>          no error message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*>          > 0:  if INFO = i, band LU factorization failed on i-th column
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hetrf_aa_2stage
*
*  =====================================================================
      SUBROUTINE zhetrf_aa_2stage( UPLO, N, A, LDA, TB, LTB, IPIV,
     $                             IPIV2, WORK, LWORK, 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..--
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            N, LDA, LTB, LWORK, INFO
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IPIV2( * )
      COMPLEX*16         A( LDA, * ), TB( * ), WORK( * )
*     ..
*
*  =====================================================================
*     .. Parameters ..
      COMPLEX*16         ZERO, ONE
      parameter( zero = ( 0.0e+0, 0.0e+0 ),
     $                     one  = ( 1.0e+0, 0.0e+0 ) )
*
*     .. Local Scalars ..
      LOGICAL            UPPER, TQUERY, WQUERY
      INTEGER            I, J, K, I1, I2, TD
      INTEGER            LWKOPT, LDTB, NB, KB, JB, NT, IINFO
      COMPLEX*16         PIV
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           lsame, ilaenv
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zcopy, zlacgv, zlacpy,
     $                   zlaset, zgbtrf, zgemm,  zgetrf, 
     $                   zhegst, zswap, ztrsm 
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dconjg, min, max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      upper = lsame( uplo, 'U' )
      wquery = ( lwork.EQ.-1 )
      tquery = ( ltb.EQ.-1 )
      IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -4
      ELSE IF( ltb.LT.max( 1, 4*n ) .AND. .NOT.tquery ) THEN
         info = -6
      ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.wquery ) THEN
         info = -10
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZHETRF_AA_2STAGE', -info )
         RETURN
      END IF
*
*     Answer the query
*
      nb = ilaenv( 1, 'ZHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
      IF( info.EQ.0 ) THEN
         IF( tquery ) THEN
            tb( 1 ) = max( 1, (3*nb+1)*n )
         END IF
         IF( wquery ) THEN
            work( 1 ) = max( 1, n*nb )
         END IF
      END IF
      IF( tquery .OR. wquery ) THEN
         RETURN
      END IF
*
*     Quick return
*
      IF( n.EQ.0 ) THEN
         RETURN
      ENDIF
*
*     Determine the number of the block size
*
      ldtb = ltb/n
      IF( ldtb .LT. 3*nb+1 ) THEN
         nb = (ldtb-1)/3
      END IF
      IF( lwork .LT. nb*n ) THEN
         nb = lwork/n
      END IF
*
*     Determine the number of the block columns
*
      nt = (n+nb-1)/nb
      td = 2*nb
      kb = min(nb, n)
*
*     Initialize vectors/matrices
*
      DO j = 1, kb
         ipiv( j ) = j
      END DO
*
*     Save NB
*
      tb( 1 ) = nb
*
      IF( upper ) THEN
*
*        .....................................................
*        Factorize A as U**H*D*U using the upper triangle of A
*        .....................................................
*
         DO j = 0, nt-1
*         
*           Generate Jth column of W and H
*
            kb = min(nb, n-j*nb)
            DO i = 1, j-1
               IF( i.EQ.1 ) THEN
*                  H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
                  IF( i .EQ. (j-1) ) THEN
                     jb = nb+kb
                  ELSE
                     jb = 2*nb
                  END IF
                  CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                    nb, kb, jb,
     $                    one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
     $                         a( (i-1)*nb+1, j*nb+1 ), lda,
     $                    zero, work( i*nb+1 ), n )
               ELSE
*                 H(I,J) = T(I,I-1)*U(I-1,J) + T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
                  IF( i .EQ. (j-1) ) THEN
                     jb = 2*nb+kb
                  ELSE
                     jb = 3*nb
                  END IF
                  CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                    nb, kb, jb,
     $                    one,  tb( td+nb+1 + ((i-1)*nb)*ldtb ),
     $                       ldtb-1,
     $                          a( (i-2)*nb+1, j*nb+1 ), lda,
     $                    zero, work( i*nb+1 ), n )
               END IF
            END DO
*         
*           Compute T(J,J)
*     
            CALL zlacpy( 'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
     $                   tb( td+1 + (j*nb)*ldtb ), ldtb-1 ) 
            IF( j.GT.1 ) THEN
*              T(J,J) = U(1:J,J)'*H(1:J)             
               CALL zgemm( 'Conjugate transpose', 'NoTranspose',
     $                 kb, kb, (j-1)*nb,
     $                -one, a( 1, j*nb+1 ), lda,
     $                      work( nb+1 ), n,
     $                 one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
*              T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
               CALL zgemm( 'Conjugate transpose', 'NoTranspose',
     $                 kb, nb, kb,
     $                 one,  a( (j-1)*nb+1, j*nb+1 ), lda,
     $                       tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
     $                 zero, work( 1 ), n )
               CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, kb, nb,
     $                -one, work( 1 ), n,
     $                      a( (j-2)*nb+1, j*nb+1 ), lda,
     $                 one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
            END IF
            IF( j.GT.0 ) THEN 
               CALL zhegst( 1, 'Upper', kb, 
     $                      tb( td+1 + (j*nb)*ldtb ), ldtb-1, 
     $                      a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
            END IF
*
*           Expand T(J,J) into full format
*
            DO i = 1, kb
               tb( td+1 + (j*nb+i-1)*ldtb )
     $            = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
               DO k = i+1, kb
                  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
     $               = dconjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
               END DO
            END DO
*
            IF( j.LT.nt-1 ) THEN
               IF( j.GT.0 ) THEN
*
*                 Compute H(J,J)
*
                  IF( j.EQ.1 ) THEN
                     CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                       kb, kb, kb,
     $                       one,  tb( td+1 + (j*nb)*ldtb ), ldtb-1,
     $                             a( (j-1)*nb+1, j*nb+1 ), lda,
     $                       zero, work( j*nb+1 ), n )
                  ELSE
                     CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                      kb, kb, nb+kb,
     $                      one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
     $                         ldtb-1,
     $                            a( (j-2)*nb+1, j*nb+1 ), lda,
     $                      zero, work( j*nb+1 ), n )
                  END IF
*
*                 Update with the previous column
*
                  CALL zgemm( 'Conjugate transpose', 'NoTranspose',
     $                    nb, n-(j+1)*nb, j*nb,
     $                    -one, work( nb+1 ), n,
     $                          a( 1, (j+1)*nb+1 ), lda,
     $                     one, a( j*nb+1, (j+1)*nb+1 ), lda )
               END IF
*
*              Copy panel to workspace to call ZGETRF
*
               DO k = 1, nb
                   CALL zcopy( n-(j+1)*nb,
     $                         a( j*nb+k, (j+1)*nb+1 ), lda,
     $                         work( 1+(k-1)*n ), 1 )
               END DO
*
*              Factorize panel
*
               CALL zgetrf( n-(j+1)*nb, nb, 
     $                      work, n,
     $                      ipiv( (j+1)*nb+1 ), iinfo )
c               IF( IINFO.NE.0 .AND. INFO.EQ.0 ) THEN
c                  INFO = IINFO+(J+1)*NB
c               END IF
*
*              Copy panel back
*
               DO k = 1, nb
*
*                  Copy only L-factor
*
                   CALL zcopy( n-k-(j+1)*nb,
     $                         work( k+1+(k-1)*n ), 1,
     $                         a( j*nb+k, (j+1)*nb+k+1 ), lda )
*
*                  Transpose U-factor to be copied back into T(J+1, J)
*
                   CALL zlacgv( k, work( 1+(k-1)*n ), 1 )
               END DO
*         
*              Compute T(J+1, J), zero out for GEMM update
*     
               kb = min(nb, n-(j+1)*nb)
               CALL zlaset( 'Full', kb, nb, zero, zero, 
     $                      tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
               CALL zlacpy( 'Upper', kb, nb,
     $                      work, n,
     $                      tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
               IF( j.GT.0 ) THEN 
                  CALL ztrsm( 'R', 'U', 'N', 'U', kb, nb, one,
     $                        a( (j-1)*nb+1, j*nb+1 ), lda,
     $                        tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
               END IF
*
*              Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
*              updates
*
               DO k = 1, nb
                  DO i = 1, kb
                     tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
     $                  = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
                  END DO
               END DO
               CALL zlaset( 'Lower', kb, nb, zero, one, 
     $                      a( j*nb+1, (j+1)*nb+1), lda )
*              
*              Apply pivots to trailing submatrix of A
*     
               DO k = 1, kb
*                 > Adjust ipiv
                  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
*                  
                  i1 = (j+1)*nb+k
                  i2 = ipiv( (j+1)*nb+k )
                  IF( i1.NE.i2 ) THEN 
*                    > Apply pivots to previous columns of L
                     CALL zswap( k-1, a( (j+1)*nb+1, i1 ), 1, 
     $                                a( (j+1)*nb+1, i2 ), 1 )
*                    > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
                     IF( i2.GT.(i1+1) ) THEN
                        CALL zswap( i2-i1-1, a( i1, i1+1 ), lda,
     $                                       a( i1+1, i2 ), 1 )
                        CALL zlacgv( i2-i1-1, a( i1+1, i2 ), 1 )
                     END IF
                     CALL zlacgv( i2-i1, a( i1, i1+1 ), lda )
*                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
                     IF( i2.LT.n )
     $                  CALL zswap( n-i2, a( i1, i2+1 ), lda,
     $                                    a( i2, i2+1 ), lda ) 
*                    > Swap A(I1, I1) with A(I2, I2)
                     piv = a( i1, i1 )
                     a( i1, i1 ) = a( i2, i2 )
                     a( i2, i2 ) = piv
*                    > Apply pivots to previous columns of L
                     IF( j.GT.0 ) THEN
                        CALL zswap( j*nb, a( 1, i1 ), 1,
     $                                    a( 1, i2 ), 1 )
                     END IF
                  ENDIF   
               END DO   
            END IF
         END DO
      ELSE
*
*        .....................................................
*        Factorize A as L*D*L**H using the lower triangle of A
*        .....................................................
*
         DO j = 0, nt-1
*         
*           Generate Jth column of W and H
*
            kb = min(nb, n-j*nb)
            DO i = 1, j-1
               IF( i.EQ.1 ) THEN
*                  H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
                  IF( i .EQ. (j-1) ) THEN
                     jb = nb+kb
                  ELSE
                     jb = 2*nb
                  END IF
                  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
     $                    nb, kb, jb,
     $                    one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
     $                         a( j*nb+1, (i-1)*nb+1 ), lda,
     $                    zero, work( i*nb+1 ), n )
               ELSE
*                 H(I,J) = T(I,I-1)*L(J,I-1)' + T(I,I)*L(J,I)' + T(I,I+1)*L(J,I+1)'
                  IF( i .EQ. (j-1) ) THEN
                     jb = 2*nb+kb
                  ELSE
                     jb = 3*nb
                  END IF
                  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
     $                    nb, kb, jb,
     $                    one,  tb( td+nb+1 + ((i-1)*nb)*ldtb ),
     $                       ldtb-1,
     $                          a( j*nb+1, (i-2)*nb+1 ), lda,
     $                    zero, work( i*nb+1 ), n )
               END IF
            END DO
*         
*           Compute T(J,J)
*     
            CALL zlacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
     $                   tb( td+1 + (j*nb)*ldtb ), ldtb-1 ) 
            IF( j.GT.1 ) THEN
*              T(J,J) = L(J,1:J)*H(1:J)             
               CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, kb, (j-1)*nb,
     $                -one, a( j*nb+1, 1 ), lda,
     $                      work( nb+1 ), n,
     $                 one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
*              T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
               CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, nb, kb,
     $                 one,  a( j*nb+1, (j-1)*nb+1 ), lda,
     $                       tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
     $                 zero, work( 1 ), n )
               CALL zgemm( 'NoTranspose', 'Conjugate transpose',
     $                 kb, kb, nb,
     $                -one, work( 1 ), n,
     $                      a( j*nb+1, (j-2)*nb+1 ), lda,
     $                 one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
            END IF
            IF( j.GT.0 ) THEN 
               CALL zhegst( 1, 'Lower', kb, 
     $                      tb( td+1 + (j*nb)*ldtb ), ldtb-1,
     $                      a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
            END IF
*
*           Expand T(J,J) into full format
*
            DO i = 1, kb
               tb( td+1 + (j*nb+i-1)*ldtb ) 
     $            = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
               DO k = i+1, kb
                  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
     $               = dconjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
               END DO
            END DO
*
            IF( j.LT.nt-1 ) THEN
               IF( j.GT.0 ) THEN
*
*                 Compute H(J,J)
*
                  IF( j.EQ.1 ) THEN
                     CALL zgemm( 'NoTranspose',
     $                           'Conjugate transpose',
     $                       kb, kb, kb,
     $                       one,  tb( td+1 + (j*nb)*ldtb ), ldtb-1,
     $                             a( j*nb+1, (j-1)*nb+1 ), lda,
     $                       zero, work( j*nb+1 ), n )
                  ELSE
                     CALL zgemm( 'NoTranspose',
     $                           'Conjugate transpose',
     $                      kb, kb, nb+kb,
     $                      one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
     $                         ldtb-1,
     $                            a( j*nb+1, (j-2)*nb+1 ), lda,
     $                      zero, work( j*nb+1 ), n )
                  END IF
*
*                 Update with the previous column
*
                  CALL zgemm( 'NoTranspose', 'NoTranspose',
     $                    n-(j+1)*nb, nb, j*nb,
     $                    -one, a( (j+1)*nb+1, 1 ), lda,
     $                          work( nb+1 ), n,
     $                     one, a( (j+1)*nb+1, j*nb+1 ), lda )
               END IF
*
*              Factorize panel
*
               CALL zgetrf( n-(j+1)*nb, nb, 
     $                      a( (j+1)*nb+1, j*nb+1 ), lda,
     $                      ipiv( (j+1)*nb+1 ), iinfo )
c               IF( IINFO.NE.0 .AND. INFO.EQ.0 ) THEN
c                  INFO = IINFO+(J+1)*NB
c               END IF
*         
*              Compute T(J+1, J), zero out for GEMM update
*     
               kb = min(nb, n-(j+1)*nb)
               CALL zlaset( 'Full', kb, nb, zero, zero, 
     $                      tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
               CALL zlacpy( 'Upper', kb, nb,
     $                      a( (j+1)*nb+1, j*nb+1 ), lda,
     $                      tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
               IF( j.GT.0 ) THEN 
                  CALL ztrsm( 'R', 'L', 'C', 'U', kb, nb, one,
     $                        a( j*nb+1, (j-1)*nb+1 ), lda,
     $                        tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
               END IF
*
*              Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
*              updates
*
               DO k = 1, nb
                  DO i = 1, kb
                     tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
     $                  = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
                  END DO
               END DO
               CALL zlaset( 'Upper', kb, nb, zero, one, 
     $                      a( (j+1)*nb+1, j*nb+1), lda )
*              
*              Apply pivots to trailing submatrix of A
*     
               DO k = 1, kb
*                 > Adjust ipiv               
                  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
*                  
                  i1 = (j+1)*nb+k
                  i2 = ipiv( (j+1)*nb+k )
                  IF( i1.NE.i2 ) THEN 
*                    > Apply pivots to previous columns of L
                     CALL zswap( k-1, a( i1, (j+1)*nb+1 ), lda, 
     $                                a( i2, (j+1)*nb+1 ), lda )
*                    > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
                     IF( i2.GT.(i1+1) ) THEN
                        CALL zswap( i2-i1-1, a( i1+1, i1 ), 1,
     $                                       a( i2, i1+1 ), lda )
                        CALL zlacgv( i2-i1-1, a( i2, i1+1 ), lda )
                     END IF
                     CALL zlacgv( i2-i1, a( i1+1, i1 ), 1 )
*                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
                     IF( i2.LT.n )
     $                  CALL zswap( n-i2, a( i2+1, i1 ), 1,
     $                                    a( i2+1, i2 ), 1 ) 
*                    > Swap A(I1, I1) with A(I2, I2)
                     piv = a( i1, i1 )
                     a( i1, i1 ) = a( i2, i2 )
                     a( i2, i2 ) = piv
*                    > Apply pivots to previous columns of L
                     IF( j.GT.0 ) THEN
                        CALL zswap( j*nb, a( i1, 1 ), lda,
     $                                    a( i2, 1 ), lda )
                     END IF
                  ENDIF   
               END DO   
*         
*              Apply pivots to previous columns of L
*         
c               CALL ZLASWP( J*NB, A( 1, 1 ), LDA, 
c     $                     (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
            END IF
         END DO
      END IF
*
*     Factor the band matrix
      CALL zgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
*
      RETURN
*
*     End of ZHETRF_AA_2STAGE
*
      END