clahef_aa.f

Go to the documentation of this file.

*> \brief \b CLAHEF_AA
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CLAHEF_AA + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahef_aa.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahef_aa.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahef_aa.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
*                             H, LDH, WORK )
*
*       .. Scalar Arguments ..
*       CHARACTER    UPLO
*       INTEGER      J1, M, NB, LDA, LDH
*       ..
*       .. Array Arguments ..
*       INTEGER      IPIV( * )
*       COMPLEX      A( LDA, * ), H( LDH, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CLAHEF_AA factorizes a panel of a complex hermitian matrix A using
*> the Aasen's algorithm. The panel consists of a set of NB rows of A
*> when UPLO is U, or a set of NB columns when UPLO is L.
*>
*> In order to factorize the panel, the Aasen's algorithm requires the
*> last row, or column, of the previous panel. The first row, or column,
*> of A is set to be the first row, or column, of an identity matrix,
*> which is used to factorize the first panel.
*>
*> The resulting J-th row of U, or J-th column of L, is stored in the
*> (J-1)-th row, or column, of A (without the unit diagonals), while
*> the diagonal and subdiagonal of A are overwritten by those of T.
*>
*> \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] J1
*> \verbatim
*>          J1 is INTEGER
*>          The location of the first row, or column, of the panel
*>          within the submatrix of A, passed to this routine, e.g.,
*>          when called by CHETRF_AA, for the first panel, J1 is 1,
*>          while for the remaining panels, J1 is 2.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The dimension of the submatrix. M >= 0.
*> \endverbatim
*>
*> \param[in] NB
*> \verbatim
*>          NB is INTEGER
*>          The dimension of the panel to be facotorized.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX array, dimension (LDA,M) for
*>          the first panel, while dimension (LDA,M+1) for the
*>          remaining panels.
*>
*>          On entry, A contains the last row, or column, of
*>          the previous panel, and the trailing submatrix of A
*>          to be factorized, except for the first panel, only
*>          the panel is passed.
*>
*>          On exit, the leading panel is factorized.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension (N)
*>          Details of the row and column interchanges,
*>          the row and column k were interchanged with the row and
*>          column IPIV(k).
*> \endverbatim
*>
*> \param[in,out] H
*> \verbatim
*>          H is COMPLEX workspace, dimension (LDH,NB).
*>
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*>          LDH is INTEGER
*>          The leading dimension of the workspace H. LDH >= max(1,M).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX workspace, dimension (M).
*> \endverbatim
*>
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complexSYcomputational
*
*  =====================================================================
      SUBROUTINE clahef_aa( UPLO, J1, M, NB, A, LDA, IPIV,
     $                      H, LDH, WORK )
*
*  -- 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      M, NB, J1, LDA, LDH
*     ..
*     .. Array Arguments ..
      INTEGER      IPIV( * )
      COMPLEX      A( LDA, * ), H( LDH, * ), WORK( * )
*     ..
*
*  =====================================================================
*     .. Parameters ..
      COMPLEX      ZERO, ONE
      parameter( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
*
*     .. Local Scalars ..
      INTEGER      J, K, K1, I1, I2, MJ
      COMPLEX      PIV, ALPHA
*     ..
*     .. External Functions ..
      LOGICAL      LSAME
      INTEGER      ICAMAX, ILAENV
      EXTERNAL     lsame, ilaenv, icamax
*     ..
*     .. External Subroutines ..
      EXTERNAL     clacgv, cgemv, cscal, caxpy, ccopy, cswap, claset,
     $             xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC    real, conjg, max
*     ..
*     .. Executable Statements ..
*
      j = 1
*
*     K1 is the first column of the panel to be factorized
*     i.e.,  K1 is 2 for the first block column, and 1 for the rest of the blocks
*
      k1 = (2-j1)+1
*
      IF( lsame( uplo, 'U' ) ) THEN
*
*        .....................................................
*        Factorize A as U**T*D*U using the upper triangle of A
*        .....................................................
*
 10      CONTINUE
         IF ( j.GT.min(m, nb) )
     $      GO TO 20
*
*        K is the column to be factorized
*         when being called from CHETRF_AA,
*         > for the first block column, J1 is 1, hence J1+J-1 is J,
*         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
*
         k = j1+j-1
         IF( j.EQ.m ) THEN
*
*            Only need to compute T(J, J)
*
             mj = 1
         ELSE
             mj = m-j+1
         END IF
*
*        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
*         where H(J:N, J) has been initialized to be A(J, J:N)
*
         IF( k.GT.2 ) THEN
*
*        K is the column to be factorized
*         > for the first block column, K is J, skipping the first two
*           columns
*         > for the rest of the columns, K is J+1, skipping only the
*           first column
*
            CALL clacgv( j-k1, a( 1, j ), 1 )
            CALL cgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( 1, j ), 1,
     $                  one, h( j, j ), 1 )
            CALL clacgv( j-k1, a( 1, j ), 1 )
         END IF
*
*        Copy H(i:n, i) into WORK
*
         CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
*
         IF( j.GT.k1 ) THEN
*
*           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
*
            alpha = -conjg( a( k-1, j ) )
            CALL caxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
         END IF
*
*        Set A(J, J) = T(J, J)
*
         a( k, j ) = real( work( 1 ) )
*
         IF( j.LT.m ) THEN
*
*           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
*            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
*
            IF( k.GT.1 ) THEN
               alpha = -a( k, j )
               CALL caxpy( m-j, alpha, a( k-1, j+1 ), lda,
     $                                 work( 2 ), 1 )
            ENDIF
*
*           Find max(|WORK(2:n)|)
*
            i2 = icamax( m-j, work( 2 ), 1 ) + 1
            piv = work( i2 )
*
*           Apply hermitian pivot
*
            IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
*
*              Swap WORK(I1) and WORK(I2)
*
               i1 = 2
               work( i2 ) = work( i1 )
               work( i1 ) = piv
*
*              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
*
               i1 = i1+j-1
               i2 = i2+j-1
               CALL cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
     $                              a( j1+i1, i2 ), 1 )
               CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
               CALL clacgv( i2-i1-1, a( j1+i1, i2 ), 1 )
*
*              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
*
               IF( i2.LT.m )
     $            CALL cswap( m-i2, a( j1+i1-1, i2+1 ), lda,
     $                              a( j1+i2-1, i2+1 ), lda )
*
*              Swap A(I1, I1) with A(I2,I2)
*
               piv = a( i1+j1-1, i1 )
               a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
               a( j1+i2-1, i2 ) = piv
*
*              Swap H(I1, 1:J1) with H(I2, 1:J1)
*
               CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
               ipiv( i1 ) = i2
*
               IF( i1.GT.(k1-1) ) THEN
*
*                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
*                  skipping the first column
*
                  CALL cswap( i1-k1+1, a( 1, i1 ), 1,
     $                                 a( 1, i2 ), 1 )
               END IF
            ELSE
               ipiv( j+1 ) = j+1
            ENDIF
*
*           Set A(J, J+1) = T(J, J+1)
*
            a( k, j+1 ) = work( 2 )
*
            IF( j.LT.nb ) THEN
*
*              Copy A(J+1:N, J+1) into H(J:N, J),
*
               CALL ccopy( m-j, a( k+1, j+1 ), lda,
     $                          h( j+1, j+1 ), 1 )
            END IF
*
*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
*
            IF( j.LT.(m-1) ) THEN
               IF( a( k, j+1 ).NE.zero ) THEN
                  alpha = one / a( k, j+1 )
                  CALL ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
                  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
               ELSE
                  CALL claset( 'Full', 1, m-j-1, zero, zero,
     $                         a( k, j+2 ), lda)
               END IF
            END IF
         END IF
         j = j + 1
         GO TO 10
 20      CONTINUE
*
      ELSE
*
*        .....................................................
*        Factorize A as L*D*L**T using the lower triangle of A
*        .....................................................
*
 30      CONTINUE
         IF( j.GT.min( m, nb ) )
     $      GO TO 40
*
*        K is the column to be factorized
*         when being called from CHETRF_AA,
*         > for the first block column, J1 is 1, hence J1+J-1 is J,
*         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
*
         k = j1+j-1
         IF( j.EQ.m ) THEN
*
*            Only need to compute T(J, J)
*
             mj = 1
         ELSE
             mj = m-j+1
         END IF
*
*        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
*         where H(J:N, J) has been initialized to be A(J:N, J)
*
         IF( k.GT.2 ) THEN
*
*        K is the column to be factorized
*         > for the first block column, K is J, skipping the first two
*           columns
*         > for the rest of the columns, K is J+1, skipping only the
*           first column
*
            CALL clacgv( j-k1, a( j, 1 ), lda )
            CALL cgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( j, 1 ), lda,
     $                  one, h( j, j ), 1 )
            CALL clacgv( j-k1, a( j, 1 ), lda )
         END IF
*
*        Copy H(J:N, J) into WORK
*
         CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
*
         IF( j.GT.k1 ) THEN
*
*           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
*            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
*
            alpha = -conjg( a( j, k-1 ) )
            CALL caxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
         END IF
*
*        Set A(J, J) = T(J, J)
*
         a( j, k ) = real( work( 1 ) )
*
         IF( j.LT.m ) THEN
*
*           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
*            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
*
            IF( k.GT.1 ) THEN
               alpha = -a( j, k )
               CALL caxpy( m-j, alpha, a( j+1, k-1 ), 1,
     $                                 work( 2 ), 1 )
            ENDIF
*
*           Find max(|WORK(2:n)|)
*
            i2 = icamax( m-j, work( 2 ), 1 ) + 1
            piv = work( i2 )
*
*           Apply hermitian pivot
*
            IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
*
*              Swap WORK(I1) and WORK(I2)
*
               i1 = 2
               work( i2 ) = work( i1 )
               work( i1 ) = piv
*
*              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
*
               i1 = i1+j-1
               i2 = i2+j-1
               CALL cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
     $                              a( i2, j1+i1 ), lda )
               CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
               CALL clacgv( i2-i1-1, a( i2, j1+i1 ), lda )
*
*              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
*
               IF( i2.LT.m )
     $            CALL cswap( m-i2, a( i2+1, j1+i1-1 ), 1,
     $                              a( i2+1, j1+i2-1 ), 1 )
*
*              Swap A(I1, I1) with A(I2, I2)
*
               piv = a( i1, j1+i1-1 )
               a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
               a( i2, j1+i2-1 ) = piv
*
*              Swap H(I1, I1:J1) with H(I2, I2:J1)
*
               CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
               ipiv( i1 ) = i2
*
               IF( i1.GT.(k1-1) ) THEN
*
*                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
*                  skipping the first column
*
                  CALL cswap( i1-k1+1, a( i1, 1 ), lda,
     $                                 a( i2, 1 ), lda )
               END IF
            ELSE
               ipiv( j+1 ) = j+1
            ENDIF
*
*           Set A(J+1, J) = T(J+1, J)
*
            a( j+1, k ) = work( 2 )
*
            IF( j.LT.nb ) THEN
*
*              Copy A(J+1:N, J+1) into H(J+1:N, J),
*
               CALL ccopy( m-j, a( j+1, k+1 ), 1,
     $                          h( j+1, j+1 ), 1 )
            END IF
*
*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
*
            IF( j.LT.(m-1) ) THEN
               IF( a( j+1, k ).NE.zero ) THEN
                  alpha = one / a( j+1, k )
                  CALL ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
                  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
               ELSE
                  CALL claset( 'Full', m-j-1, 1, zero, zero,
     $                         a( j+2, k ), lda )
               END IF
            END IF
         END IF
         j = j + 1
         GO TO 30
 40      CONTINUE
      END IF
      RETURN
*
*     End of CLAHEF_AA
*
      END