zpftri()

subroutine zpftri	(	character	transr,
		character	uplo,
		integer	n,
		complex*16, dimension( 0: * )	a,
		integer	info )

ZPFTRI

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

Purpose:

!>
!> ZPFTRI computes the inverse of a complex Hermitian positive definite
!> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
!> computed by ZPFTRF.
!> 

Parameters

[in]	TRANSR	!> TRANSR is CHARACTER*1 !> = 'N': The Normal TRANSR of RFP A is stored; !> = 'C': The Conjugate-transpose TRANSR of RFP A is stored. !>
[in]	UPLO	!> UPLO is CHARACTER*1 !> = 'U': Upper triangle of A is stored; !> = 'L': Lower triangle of A is stored. !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in,out]	A	!> A is COMPLEX*16 array, dimension ( N*(N+1)/2 ); !> On entry, the Hermitian matrix A in RFP format. RFP format is !> described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' !> then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is !> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is !> the Conjugate-transpose of RFP A as defined when !> TRANSR = 'N'. The contents of RFP A are defined by UPLO as !> follows: If UPLO = 'U' the RFP A contains the nt elements of !> upper packed A. If UPLO = 'L' the RFP A contains the elements !> of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = !> 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N !> is odd. See the Note below for more details. !> !> On exit, the Hermitian inverse of the original matrix, in the !> same storage format. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !> > 0: if INFO = i, the (i,i) element of the factor U or L is !> zero, and the inverse could not be computed. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  We first consider Standard Packed Format when N is even.
!>  We give an example where N = 6.
!>
!>      AP is Upper             AP is Lower
!>
!>   00 01 02 03 04 05       00
!>      11 12 13 14 15       10 11
!>         22 23 24 25       20 21 22
!>            33 34 35       30 31 32 33
!>               44 45       40 41 42 43 44
!>                  55       50 51 52 53 54 55
!>
!>
!>  Let TRANSR = 'N'. RFP holds AP as follows:
!>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
!>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
!>  conjugate-transpose of the first three columns of AP upper.
!>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
!>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
!>  conjugate-transpose of the last three columns of AP lower.
!>  To denote conjugate we place -- above the element. This covers the
!>  case N even and TRANSR = 'N'.
!>
!>         RFP A                   RFP A
!>
!>                                -- -- --
!>        03 04 05                33 43 53
!>                                   -- --
!>        13 14 15                00 44 54
!>                                      --
!>        23 24 25                10 11 55
!>
!>        33 34 35                20 21 22
!>        --
!>        00 44 45                30 31 32
!>        -- --
!>        01 11 55                40 41 42
!>        -- -- --
!>        02 12 22                50 51 52
!>
!>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
!>  transpose of RFP A above. One therefore gets:
!>
!>
!>           RFP A                   RFP A
!>
!>     -- -- -- --                -- -- -- -- -- --
!>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
!>     -- -- -- -- --                -- -- -- -- --
!>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
!>     -- -- -- -- -- --                -- -- -- --
!>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
!>
!>
!>  We next  consider Standard Packed Format when N is odd.
!>  We give an example where N = 5.
!>
!>     AP is Upper                 AP is Lower
!>
!>   00 01 02 03 04              00
!>      11 12 13 14              10 11
!>         22 23 24              20 21 22
!>            33 34              30 31 32 33
!>               44              40 41 42 43 44
!>
!>
!>  Let TRANSR = 'N'. RFP holds AP as follows:
!>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
!>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
!>  conjugate-transpose of the first two   columns of AP upper.
!>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
!>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
!>  conjugate-transpose of the last two   columns of AP lower.
!>  To denote conjugate we place -- above the element. This covers the
!>  case N odd  and TRANSR = 'N'.
!>
!>         RFP A                   RFP A
!>
!>                                   -- --
!>        02 03 04                00 33 43
!>                                      --
!>        12 13 14                10 11 44
!>
!>        22 23 24                20 21 22
!>        --
!>        00 33 34                30 31 32
!>        -- --
!>        01 11 44                40 41 42
!>
!>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
!>  transpose of RFP A above. One therefore gets:
!>
!>
!>           RFP A                   RFP A
!>
!>     -- -- --                   -- -- -- -- -- --
!>     02 12 22 00 01             00 10 20 30 40 50
!>     -- -- -- --                   -- -- -- -- --
!>     03 13 23 33 11             33 11 21 31 41 51
!>     -- -- -- -- --                   -- -- -- --
!>     04 14 24 34 44             43 44 22 32 42 52
!> 

Definition at line 209 of file zpftri.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..--
*
*     .. Scalar Arguments ..
      CHARACTER          TRANSR, UPLO
      INTEGER            INFO, N
*     .. Array Arguments ..
      COMPLEX*16         A( 0: * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE
      COMPLEX*16         CONE
      parameter( one = 1.d0, cone = ( 1.d0, 0.d0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            LOWER, NISODD, NORMALTRANSR
      INTEGER            N1, N2, K
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, ztftri, zlauum, ztrmm,
     $                   zherk
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          mod
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      normaltransr = lsame( transr, 'N' )
      lower = lsame( uplo, 'L' )
      IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
         info = -1
      ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZPFTRI', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Invert the triangular Cholesky factor U or L.
*
      CALL ztftri( transr, uplo, 'N', n, a, info )
      IF( info.GT.0 )
     $   RETURN
*
*     If N is odd, set NISODD = .TRUE.
*     If N is even, set K = N/2 and NISODD = .FALSE.
*
      IF( mod( n, 2 ).EQ.0 ) THEN
         k = n / 2
         nisodd = .false.
      ELSE
         nisodd = .true.
      END IF
*
*     Set N1 and N2 depending on LOWER
*
      IF( lower ) THEN
         n2 = n / 2
         n1 = n - n2
      ELSE
         n1 = n / 2
         n2 = n - n1
      END IF
*
*     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
*     inv(L)^C*inv(L). There are eight cases.
*
      IF( nisodd ) THEN
*
*        N is odd
*
         IF( normaltransr ) THEN
*
*           N is odd and TRANSR = 'N'
*
            IF( lower ) THEN
*
*              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
*              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
*              T1 -> a(0), T2 -> a(n), S -> a(N1)
*
               CALL zlauum( 'L', n1, a( 0 ), n, info )
               CALL zherk( 'L', 'C', n1, n2, one, a( n1 ), n, one,
     $                     a( 0 ), n )
               CALL ztrmm( 'L', 'U', 'N', 'N', n2, n1, cone, a( n ),
     $                     n,
     $                     a( n1 ), n )
               CALL zlauum( 'U', n2, a( n ), n, info )
*
            ELSE
*
*              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
*              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
*              T1 -> a(N2), T2 -> a(N1), S -> a(0)
*
               CALL zlauum( 'L', n1, a( n2 ), n, info )
               CALL zherk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
     $                     a( n2 ), n )
               CALL ztrmm( 'R', 'U', 'C', 'N', n1, n2, cone, a( n1 ),
     $                     n,
     $                     a( 0 ), n )
               CALL zlauum( 'U', n2, a( n1 ), n, info )
*
            END IF
*
         ELSE
*
*           N is odd and TRANSR = 'C'
*
            IF( lower ) THEN
*
*              SRPA for LOWER, TRANSPOSE, and N is odd
*              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
*
               CALL zlauum( 'U', n1, a( 0 ), n1, info )
               CALL zherk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1,
     $                     one,
     $                     a( 0 ), n1 )
               CALL ztrmm( 'R', 'L', 'N', 'N', n1, n2, cone, a( 1 ),
     $                     n1,
     $                     a( n1*n1 ), n1 )
               CALL zlauum( 'L', n2, a( 1 ), n1, info )
*
            ELSE
*
*              SRPA for UPPER, TRANSPOSE, and N is odd
*              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
*
               CALL zlauum( 'U', n1, a( n2*n2 ), n2, info )
               CALL zherk( 'U', 'C', n1, n2, one, a( 0 ), n2, one,
     $                     a( n2*n2 ), n2 )
               CALL ztrmm( 'L', 'L', 'C', 'N', n2, n1, cone,
     $                     a( n1*n2 ),
     $                     n2, a( 0 ), n2 )
               CALL zlauum( 'L', n2, a( n1*n2 ), n2, info )
*
            END IF
*
         END IF
*
      ELSE
*
*        N is even
*
         IF( normaltransr ) THEN
*
*           N is even and TRANSR = 'N'
*
            IF( lower ) THEN
*
*              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
*              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
*              T1 -> a(1), T2 -> a(0), S -> a(k+1)
*
               CALL zlauum( 'L', k, a( 1 ), n+1, info )
               CALL zherk( 'L', 'C', k, k, one, a( k+1 ), n+1, one,
     $                     a( 1 ), n+1 )
               CALL ztrmm( 'L', 'U', 'N', 'N', k, k, cone, a( 0 ),
     $                     n+1,
     $                     a( k+1 ), n+1 )
               CALL zlauum( 'U', k, a( 0 ), n+1, info )
*
            ELSE
*
*              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
*              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
*              T1 -> a(k+1), T2 -> a(k), S -> a(0)
*
               CALL zlauum( 'L', k, a( k+1 ), n+1, info )
               CALL zherk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
     $                     a( k+1 ), n+1 )
               CALL ztrmm( 'R', 'U', 'C', 'N', k, k, cone, a( k ),
     $                     n+1,
     $                     a( 0 ), n+1 )
               CALL zlauum( 'U', k, a( k ), n+1, info )
*
            END IF
*
         ELSE
*
*           N is even and TRANSR = 'C'
*
            IF( lower ) THEN
*
*              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
*              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
*              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
*
               CALL zlauum( 'U', k, a( k ), k, info )
               CALL zherk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k,
     $                     one,
     $                     a( k ), k )
               CALL ztrmm( 'R', 'L', 'N', 'N', k, k, cone, a( 0 ), k,
     $                     a( k*( k+1 ) ), k )
               CALL zlauum( 'L', k, a( 0 ), k, info )
*
            ELSE
*
*              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
*              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
*              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
*
               CALL zlauum( 'U', k, a( k*( k+1 ) ), k, info )
               CALL zherk( 'U', 'C', k, k, one, a( 0 ), k, one,
     $                     a( k*( k+1 ) ), k )
               CALL ztrmm( 'L', 'L', 'C', 'N', k, k, cone, a( k*k ),
     $                     k,
     $                     a( 0 ), k )
               CALL zlauum( 'L', k, a( k*k ), k, info )
*
            END IF
*
         END IF
*
      END IF
*
      RETURN
*
*     End of ZPFTRI
*

Here is the call graph for this function:

Here is the caller graph for this function: