clahef_aa()

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

CLAHEF_AA

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

Purpose:

 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.

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 CHETRF_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 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 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 workspace, dimension (M).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 142 of file clahef_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      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
*

Here is the call graph for this function:

Here is the caller graph for this function: