clasyf_aa()

subroutine clasyf_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
	)

CLASYF_AA

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

Purpose:

 DLATRF_AA factorizes a panel of a complex symmetric 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 CSYTRF_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,M).
[out]	IPIV	IPIV is INTEGER array, dimension (M) 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 clasyf_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, one = 1.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           caxpy, cgemv, cscal, ccopy, cswap, claset,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          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 CSYTRF_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:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
*         where H(J:M, J) has been initialized to be A(J, J:M)
*
         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 cgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( 1, j ), 1,
     $                  one, h( j, j ), 1 )
         END IF
*
*        Copy H(i:M, 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:M) * T(J-1,J),
*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
*
            alpha = -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 ) = work( 1 )
*
         IF( j.LT.m ) THEN
*
*           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
*            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
*
            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:M)|)
*
            i2 = icamax( m-j, work( 2 ), 1 ) + 1
            piv = work( i2 )
*
*           Apply symmetric 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:M) with A(I1+1:M, 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 )
*
*              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
*
               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:M, J+1) into H(J:M, 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:M ) / T(J, J+1),
*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, 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 CSYTRF_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:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
*         where H(J:M, J) has been initialized to be A(J:M, 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 cgemv( 'No transpose', mj, j-k1,
     $                 -one, h( j, k1 ), ldh,
     $                       a( j, 1 ), lda,
     $                  one, h( j, j ), 1 )
         END IF
*
*        Copy H(J:M, J) into WORK
*
         CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
*
         IF( j.GT.k1 ) THEN
*
*           Compute WORK := WORK - L(J:M, 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 = -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 ) = work( 1 )
*
         IF( j.LT.m ) THEN
*
*           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
*            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, 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:M)|)
*
            i2 = icamax( m-j, work( 2 ), 1 ) + 1
            piv = work( i2 )
*
*           Apply symmetric 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:M, I1) with A(I2, I1+1:M)
*
               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 )
*
*              Swap A(I2+1:M, I1) with A(I2+1:M, 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:M, J+1) into H(J+1:M, 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:M ) / T(J, J+1),
*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, 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 CLASYF_AA
*

Here is the call graph for this function:

Here is the caller graph for this function: