zlahef_aa()

subroutine zlahef_aa	(	character	uplo,
		integer	j1,
		integer	m,
		integer	nb,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		integer, dimension( * )	ipiv,
		complex*16, dimension( ldh, * )	h,
		integer	ldh,
		complex*16, dimension( * )	work
	)

ZLAHEF_AA

Download ZLAHEF_AA + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 DLAHEF_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.

Parameters

[in]	UPLO	UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.
[in]	J1	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 ZHETRF_AA, for the first panel, J1 is 1, while for the remaining panels, J1 is 2.
[in]	M	M is INTEGER The dimension of the submatrix. M >= 0.
[in]	NB	NB is INTEGER The dimension of the panel to be facotorized.
[in,out]	A	A is COMPLEX*16 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.
[in]	LDA	LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).
[out]	IPIV	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).
[in,out]	H	H is COMPLEX*16 workspace, dimension (LDH,NB).
[in]	LDH	LDH is INTEGER The leading dimension of the workspace H. LDH >= max(1,M).
[out]	WORK	WORK is COMPLEX*16 workspace, dimension (M).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 142 of file zlahef_aa.f.

*
*  -- 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*16   A( LDA, * ), H( LDH, * ), WORK( * )
*     ..
*
*  =====================================================================
*     .. Parameters ..
      COMPLEX*16   ZERO, ONE
      parameter( zero = (0.0d+0, 0.0d+0), one = (1.0d+0, 0.0d+0) )
*
*     .. Local Scalars ..
      INTEGER      J, K, K1, I1, I2, MJ
      COMPLEX*16   PIV, ALPHA
*     ..
*     .. External Functions ..
      LOGICAL      LSAME
      INTEGER      IZAMAX, ILAENV
      EXTERNAL     lsame, ilaenv, izamax
*     ..
*     .. External Subroutines ..
      EXTERNAL     zgemm, zgemv, zaxpy, zlacgv, zcopy, zscal, zswap,
     $             zlaset, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC    dble, dconjg, 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 ZHETRF_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 zlacgv( j-k1, a( 1, j ), 1 )
            CALL zgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( 1, j ), 1,
     $                  one, h( j, j ), 1 )
            CALL zlacgv( j-k1, a( 1, j ), 1 )
         END IF
*
*        Copy H(i:n, i) into WORK
*
         CALL zcopy( 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 = -dconjg( a( k-1, j ) )
            CALL zaxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
         END IF
*
*        Set A(J, J) = T(J, J)
*
         a( k, j ) = dble( 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 zaxpy( m-j, alpha, a( k-1, j+1 ), lda,
     $                                 work( 2 ), 1 )
            ENDIF
*
*           Find max(|WORK(2:n)|)
*
            i2 = izamax( 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 zswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
     $                              a( j1+i1, i2 ), 1 )
               CALL zlacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
               CALL zlacgv( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
                  CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
               ELSE
                  CALL zlaset( '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 ZHETRF_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 zlacgv( j-k1, a( j, 1 ), lda )
            CALL zgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( j, 1 ), lda,
     $                  one, h( j, j ), 1 )
            CALL zlacgv( j-k1, a( j, 1 ), lda )
         END IF
*
*        Copy H(J:N, J) into WORK
*
         CALL zcopy( 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 = -dconjg( a( j, k-1 ) )
            CALL zaxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
         END IF
*
*        Set A(J, J) = T(J, J)
*
         a( j, k ) = dble( 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 zaxpy( m-j, alpha, a( j+1, k-1 ), 1,
     $                                 work( 2 ), 1 )
            ENDIF
*
*           Find max(|WORK(2:n)|)
*
            i2 = izamax( 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 zswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
     $                              a( i2, j1+i1 ), lda )
               CALL zlacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
               CALL zlacgv( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
                  CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
               ELSE
                  CALL zlaset( '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 ZLAHEF_AA
*

Here is the call graph for this function:

Here is the caller graph for this function: