csytrf_aa_2stage.f

Go to the documentation of this file.

*> \brief \b CSYTRF_AA_2STAGE
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CSYTRF_AA_2STAGE + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE CSYTRF_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            A( LDA, * ), TB( * ), WORK( * )
*       ..
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CSYTRF_AA_2STAGE computes the factorization of a complex symmetric matrix A
*> using the Aasen's algorithm.  The form of the factorization is
*>
*>    A = U**T*T*U  or  A = L*T*L**T
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a complex symmetric 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 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 subdiaonal 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 array, dimension (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 >= 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 workspace of size LWORK
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The size of WORK. LWORK >= 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 complexSYcomputational
*
*  =====================================================================
      SUBROUTINE csytrf_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            A( LDA, * ), TB( * ), WORK( * )
*     ..
*
*  =====================================================================
*     .. Parameters ..
      COMPLEX            CZERO, CONE
      parameter( czero = ( 0.0e+0, 0.0e+0 ),
     $                     cone  = ( 1.0e+0, 0.0e+0 ) )
*
*     .. Local Scalars ..
      LOGICAL            UPPER, TQUERY, WQUERY
      INTEGER            I, J, K, I1, I2, TD
      INTEGER            LDTB, NB, KB, JB, NT, IINFO
      COMPLEX            PIV
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           lsame, ilaenv
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, cgbtrf, cgemm, cgetrf, clacpy,  
     $                   claset, ctrsm, cswap, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          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. 4*n .AND. .NOT.tquery ) THEN
         info = -6
      ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
         info = -10
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CSYTRF_AA_2STAGE', -info )
         RETURN
      END IF
*
*     Answer the query
*
      nb = ilaenv( 1, 'CSYTRF_AA_2STAGE', uplo, n, -1, -1, -1 )
      IF( info.EQ.0 ) THEN
         IF( tquery ) THEN
            tb( 1 ) = (3*nb+1)*n
         END IF
         IF( wquery ) THEN
            work( 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**T*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 cgemm( 'NoTranspose', 'NoTranspose',
     $                    nb, kb, jb,
     $                    cone,  tb( td+1 + (i*nb)*ldtb ), ldtb-1,
     $                           a( (i-1)*nb+1, j*nb+1 ), lda,
     $                    czero, 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 cgemm( 'NoTranspose', 'NoTranspose',
     $                    nb, kb, jb,
     $                    cone,  tb( td+nb+1 + ((i-1)*nb)*ldtb ),
     $                       ldtb-1,
     $                           a( (i-2)*nb+1, j*nb+1 ), lda,
     $                    czero, work( i*nb+1 ), n )
               END IF
            END DO
*         
*           Compute T(J,J)
*     
            CALL clacpy( '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 cgemm( 'Transpose', 'NoTranspose',
     $                 kb, kb, (j-1)*nb,
     $                -cone, a( 1, j*nb+1 ), lda,
     $                       work( nb+1 ), n,
     $                 cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
*              T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
               CALL cgemm( 'Transpose', 'NoTranspose',
     $                 kb, nb, kb,
     $                 cone,  a( (j-1)*nb+1, j*nb+1 ), lda,
     $                        tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
     $                 czero, work( 1 ), n )
               CALL cgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, kb, nb,
     $                -cone, work( 1 ), n,
     $                       a( (j-2)*nb+1, j*nb+1 ), lda,
     $                 cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
            END IF
*
*           Expand T(J,J) into full format
*
            DO i = 1, kb
               DO k = i+1, kb
                  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
     $                = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
               END DO
            END DO
            IF( j.GT.0 ) THEN 
c               CALL CHEGST( 1, 'Upper', KB, 
c     $                      TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 
c     $                      A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
               CALL ctrsm( 'L', 'U', 'T', 'N', kb, kb, cone,
     $                     a( (j-1)*nb+1, j*nb+1 ), lda, 
     $                     tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
               CALL ctrsm( 'R', 'U', 'N', 'n', KB, KB, CONE,
     $                     A( (J-1)*NB+1, J*NB+1 ), LDA, 
     $                     TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
            END IF
*
.LT.            IF( JNT-1 ) THEN
.GT.               IF( J0 ) THEN
*
*                 Compute H(J,J)
*
.EQ.                  IF( J1 ) THEN
                     CALL CGEMM( 'notranspose', 'notranspose',
     $                       KB, KB, KB,
     $                       CONE,  TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
     $                              A( (J-1)*NB+1, J*NB+1 ), LDA,
     $                       CZERO, WORK( J*NB+1 ), N )
                  ELSE
                     CALL CGEMM( 'notranspose', 'notranspose',
     $                      KB, KB, NB+KB,
     $                      CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
     $                         LDTB-1,
     $                             A( (J-2)*NB+1, J*NB+1 ), LDA,
     $                      CZERO, WORK( J*NB+1 ), N )
                  END IF
*
*                 Update with the previous column
*
                  CALL CGEMM( 'transpose', 'notranspose',
     $                    NB, N-(J+1)*NB, J*NB,
     $                    -CONE, WORK( NB+1 ), N,
     $                           A( 1, (J+1)*NB+1 ), LDA,
     $                     CONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
               END IF
*
*              Copy panel to workspace to call CGETRF
*
               DO K = 1, NB
                   CALL CCOPY( N-(J+1)*NB,
     $                         A( J*NB+K, (J+1)*NB+1 ), LDA,
     $                         WORK( 1+(K-1)*N ), 1 )
               END DO
*
*              Factorize panel
*
               CALL CGETRF( 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
                   CALL CCOPY( N-(J+1)*NB,
     $                         WORK( 1+(K-1)*N ), 1,
     $                         A( J*NB+K, (J+1)*NB+1 ), LDA )
               END DO
*         
*              Compute T(J+1, J), zero out for GEMM update
*     
               KB = MIN(NB, N-(J+1)*NB)
               CALL CLASET( 'full', KB, NB, CZERO, CZERO, 
     $                      TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
               CALL CLACPY( 'upper', KB, NB,
     $                      WORK, N,
     $                      TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
.GT.               IF( J0 ) THEN 
                  CALL CTRSM( 'r', 'u', 'n', 'u', KB, NB, CONE,
     $                        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 )
     $                  = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
                  END DO
               END DO
               CALL CLASET( 'lower', KB, NB, CZERO, CONE, 
     $                      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 )
.NE.                  IF( I1I2 ) THEN 
*                    > Apply pivots to previous columns of L
                     CALL CSWAP( 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)
.GT.                     IF( I2(I1+1) )
     $                  CALL CSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
     $                                       A( I1+1, I2 ), 1 )
*                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
.LT.                     IF( I2N )
     $                  CALL CSWAP( 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
.GT.                     IF( J0 ) THEN
                        CALL CSWAP( 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**T 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
.EQ.               IF( I1 ) THEN
*                  H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
.EQ.                  IF( I  (J-1) ) THEN
                     JB = NB+KB
                  ELSE
                     JB = 2*NB
                  END IF
                  CALL CGEMM( 'notranspose', 'transpose',
     $                    NB, KB, JB,
     $                    CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
     $                          A( J*NB+1, (I-1)*NB+1 ), LDA,
     $                    CZERO, 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)'
.EQ.                  IF( I  (J-1) ) THEN
                     JB = 2*NB+KB
                  ELSE
                     JB = 3*NB
                  END IF
                  CALL CGEMM( 'notranspose', 'Transpose',
     $                    nb, kb, jb,
     $                    cone,  tb( td+nb+1 + ((i-1)*nb)*ldtb ),
     $                       ldtb-1,
     $                           a( j*nb+1, (i-2)*nb+1 ), lda,
     $                    czero, work( i*nb+1 ), n )
               END IF
            END DO
*         
*           Compute T(J,J)
*     
            CALL clacpy( '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 cgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, kb, (j-1)*nb,
     $                -cone, a( j*nb+1, 1 ), lda,
     $                       work( nb+1 ), n,
     $                 cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
*              T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
               CALL cgemm( 'NoTranspose', 'NoTranspose',
     $                 kb, nb, kb,
     $                 cone,  a( j*nb+1, (j-1)*nb+1 ), lda,
     $                        tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
     $                 czero, work( 1 ), n )
               CALL cgemm( 'NoTranspose', 'Transpose',
     $                 kb, kb, nb,
     $                -cone, work( 1 ), n,
     $                       a( j*nb+1, (j-2)*nb+1 ), lda,
     $                 cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
            END IF
*
*           Expand T(J,J) into full format
*
            DO i = 1, kb
               DO k = i+1, kb
                  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
     $                = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
               END DO
            END DO
            IF( j.GT.0 ) THEN 
c               CALL CHEGST( 1, 'Lower', KB, 
c     $                      TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
c     $                      A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
               CALL ctrsm( 'L', 'L', 'N', 'N', kb, kb, cone,
     $                     a( j*nb+1, (j-1)*nb+1 ), lda,
     $                     tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
               CALL ctrsm( 'R', 'L', 'T', 'N', kb, kb, cone,
     $                     a( j*nb+1, (j-1)*nb+1 ), lda,
     $                     tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
            END IF
*
*           Symmetrize T(J,J)
*
            DO i = 1, kb
               DO k = i+1, kb
                  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
     $                = 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 cgemm( 'NoTranspose', 'Transpose',
     $                       kb, kb, kb,
     $                       cone,  tb( td+1 + (j*nb)*ldtb ), ldtb-1,
     $                              a( j*nb+1, (j-1)*nb+1 ), lda,
     $                       czero, work( j*nb+1 ), n )
                  ELSE
                     CALL cgemm( 'NoTranspose', 'Transpose',
     $                      kb, kb, nb+kb,
     $                      cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
     $                         ldtb-1,
     $                             a( j*nb+1, (j-2)*nb+1 ), lda,
     $                      czero, work( j*nb+1 ), n )
                  END IF
*
*                 Update with the previous column
*
                  CALL cgemm( 'NoTranspose', 'NoTranspose',
     $                    n-(j+1)*nb, nb, j*nb,
     $                    -cone, a( (j+1)*nb+1, 1 ), lda,
     $                           work( nb+1 ), n,
     $                     cone, a( (j+1)*nb+1, j*nb+1 ), lda )
               END IF
*
*              Factorize panel
*
               CALL cgetrf( 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 claset( 'Full', kb, nb, czero, czero, 
     $                      tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
               CALL clacpy( '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 ctrsm( 'R', 'L', 'T', 'U', kb, nb, cone,
     $                        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 ) =
     $                  tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
                  END DO
               END DO
               CALL claset( 'Upper', kb, nb, czero, cone, 
     $                      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 cswap( 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) )
     $                  CALL cswap( i2-i1-1, a( i1+1, i1 ), 1,
     $                                       a( i2, i1+1 ), lda )
*                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
                     IF( i2.LT.n )
     $                  CALL cswap( 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 cswap( j*nb, a( i1, 1 ), lda,
     $                                    a( i2, 1 ), lda )
                     END IF
                  ENDIF   
               END DO   
*         
*              Apply pivots to previous columns of L
*         
c               CALL CLASWP( 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 cgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
*
      RETURN
*
*     End of CSYTRF_AA_2STAGE
*
      END