LAPACK 3.3.0

zpftri.f

Go to the documentation of this file.
00001       SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.0)                                    --
00004 *
00005 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
00006 *     November 2010
00007 *
00008 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00009 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          TRANSR, UPLO
00013       INTEGER            INFO, N
00014 *     .. Array Arguments ..
00015       COMPLEX*16         A( 0: * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZPFTRI computes the inverse of a complex Hermitian positive definite
00022 *  matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
00023 *  computed by ZPFTRF.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  TRANSR    (input) CHARACTER*1
00029 *          = 'N':  The Normal TRANSR of RFP A is stored;
00030 *          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
00031 *
00032 *  UPLO    (input) CHARACTER*1
00033 *          = 'U':  Upper triangle of A is stored;
00034 *          = 'L':  Lower triangle of A is stored.
00035 *
00036 *  N       (input) INTEGER
00037 *          The order of the matrix A.  N >= 0.
00038 *
00039 *  A       (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 );
00040 *          On entry, the Hermitian matrix A in RFP format. RFP format is
00041 *          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
00042 *          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
00043 *          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
00044 *          the Conjugate-transpose of RFP A as defined when
00045 *          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
00046 *          follows: If UPLO = 'U' the RFP A contains the nt elements of
00047 *          upper packed A. If UPLO = 'L' the RFP A contains the elements
00048 *          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
00049 *          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
00050 *          is odd. See the Note below for more details.
00051 *
00052 *          On exit, the Hermitian inverse of the original matrix, in the
00053 *          same storage format.
00054 *
00055 *  INFO    (output) INTEGER
00056 *          = 0:  successful exit
00057 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00058 *          > 0:  if INFO = i, the (i,i) element of the factor U or L is
00059 *                zero, and the inverse could not be computed.
00060 *
00061 *  Further Details
00062 *  ===============
00063 *
00064 *  We first consider Standard Packed Format when N is even.
00065 *  We give an example where N = 6.
00066 *
00067 *      AP is Upper             AP is Lower
00068 *
00069 *   00 01 02 03 04 05       00
00070 *      11 12 13 14 15       10 11
00071 *         22 23 24 25       20 21 22
00072 *            33 34 35       30 31 32 33
00073 *               44 45       40 41 42 43 44
00074 *                  55       50 51 52 53 54 55
00075 *
00076 *
00077 *  Let TRANSR = 'N'. RFP holds AP as follows:
00078 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
00079 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
00080 *  conjugate-transpose of the first three columns of AP upper.
00081 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
00082 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
00083 *  conjugate-transpose of the last three columns of AP lower.
00084 *  To denote conjugate we place -- above the element. This covers the
00085 *  case N even and TRANSR = 'N'.
00086 *
00087 *         RFP A                   RFP A
00088 *
00089 *                                -- -- --
00090 *        03 04 05                33 43 53
00091 *                                   -- --
00092 *        13 14 15                00 44 54
00093 *                                      --
00094 *        23 24 25                10 11 55
00095 *
00096 *        33 34 35                20 21 22
00097 *        --
00098 *        00 44 45                30 31 32
00099 *        -- --
00100 *        01 11 55                40 41 42
00101 *        -- -- --
00102 *        02 12 22                50 51 52
00103 *
00104 *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
00105 *  transpose of RFP A above. One therefore gets:
00106 *
00107 *
00108 *           RFP A                   RFP A
00109 *
00110 *     -- -- -- --                -- -- -- -- -- --
00111 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
00112 *     -- -- -- -- --                -- -- -- -- --
00113 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
00114 *     -- -- -- -- -- --                -- -- -- --
00115 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
00116 *
00117 *
00118 *  We next  consider Standard Packed Format when N is odd.
00119 *  We give an example where N = 5.
00120 *
00121 *     AP is Upper                 AP is Lower
00122 *
00123 *   00 01 02 03 04              00
00124 *      11 12 13 14              10 11
00125 *         22 23 24              20 21 22
00126 *            33 34              30 31 32 33
00127 *               44              40 41 42 43 44
00128 *
00129 *
00130 *  Let TRANSR = 'N'. RFP holds AP as follows:
00131 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
00132 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
00133 *  conjugate-transpose of the first two   columns of AP upper.
00134 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
00135 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
00136 *  conjugate-transpose of the last two   columns of AP lower.
00137 *  To denote conjugate we place -- above the element. This covers the
00138 *  case N odd  and TRANSR = 'N'.
00139 *
00140 *         RFP A                   RFP A
00141 *
00142 *                                   -- --
00143 *        02 03 04                00 33 43
00144 *                                      --
00145 *        12 13 14                10 11 44
00146 *
00147 *        22 23 24                20 21 22
00148 *        --
00149 *        00 33 34                30 31 32
00150 *        -- --
00151 *        01 11 44                40 41 42
00152 *
00153 *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
00154 *  transpose of RFP A above. One therefore gets:
00155 *
00156 *
00157 *           RFP A                   RFP A
00158 *
00159 *     -- -- --                   -- -- -- -- -- --
00160 *     02 12 22 00 01             00 10 20 30 40 50
00161 *     -- -- -- --                   -- -- -- -- --
00162 *     03 13 23 33 11             33 11 21 31 41 51
00163 *     -- -- -- -- --                   -- -- -- --
00164 *     04 14 24 34 44             43 44 22 32 42 52
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ONE
00170       COMPLEX*16         CONE
00171       PARAMETER          ( ONE = 1.D0, CONE = ( 1.D0, 0.D0 ) )
00172 *     ..
00173 *     .. Local Scalars ..
00174       LOGICAL            LOWER, NISODD, NORMALTRANSR
00175       INTEGER            N1, N2, K
00176 *     ..
00177 *     .. External Functions ..
00178       LOGICAL            LSAME
00179       EXTERNAL           LSAME
00180 *     ..
00181 *     .. External Subroutines ..
00182       EXTERNAL           XERBLA, ZTFTRI, ZLAUUM, ZTRMM, ZHERK
00183 *     ..
00184 *     .. Intrinsic Functions ..
00185       INTRINSIC          MOD
00186 *     ..
00187 *     .. Executable Statements ..
00188 *
00189 *     Test the input parameters.
00190 *
00191       INFO = 0
00192       NORMALTRANSR = LSAME( TRANSR, 'N' )
00193       LOWER = LSAME( UPLO, 'L' )
00194       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
00195          INFO = -1
00196       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00197          INFO = -2
00198       ELSE IF( N.LT.0 ) THEN
00199          INFO = -3
00200       END IF
00201       IF( INFO.NE.0 ) THEN
00202          CALL XERBLA( 'ZPFTRI', -INFO )
00203          RETURN
00204       END IF
00205 *
00206 *     Quick return if possible
00207 *
00208       IF( N.EQ.0 )
00209      +   RETURN
00210 *
00211 *     Invert the triangular Cholesky factor U or L.
00212 *
00213       CALL ZTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
00214       IF( INFO.GT.0 )
00215      +   RETURN
00216 *
00217 *     If N is odd, set NISODD = .TRUE.
00218 *     If N is even, set K = N/2 and NISODD = .FALSE.
00219 *
00220       IF( MOD( N, 2 ).EQ.0 ) THEN
00221          K = N / 2
00222          NISODD = .FALSE.
00223       ELSE
00224          NISODD = .TRUE.
00225       END IF
00226 *
00227 *     Set N1 and N2 depending on LOWER
00228 *
00229       IF( LOWER ) THEN
00230          N2 = N / 2
00231          N1 = N - N2
00232       ELSE
00233          N1 = N / 2
00234          N2 = N - N1
00235       END IF
00236 *
00237 *     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
00238 *     inv(L)^C*inv(L). There are eight cases.
00239 *
00240       IF( NISODD ) THEN
00241 *
00242 *        N is odd
00243 *
00244          IF( NORMALTRANSR ) THEN
00245 *
00246 *           N is odd and TRANSR = 'N'
00247 *
00248             IF( LOWER ) THEN
00249 *
00250 *              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
00251 *              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
00252 *              T1 -> a(0), T2 -> a(n), S -> a(N1)
00253 *
00254                CALL ZLAUUM( 'L', N1, A( 0 ), N, INFO )
00255                CALL ZHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE,
00256      +                     A( 0 ), N )
00257                CALL ZTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N,
00258      +                     A( N1 ), N )
00259                CALL ZLAUUM( 'U', N2, A( N ), N, INFO )
00260 *
00261             ELSE
00262 *
00263 *              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
00264 *              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
00265 *              T1 -> a(N2), T2 -> a(N1), S -> a(0)
00266 *
00267                CALL ZLAUUM( 'L', N1, A( N2 ), N, INFO )
00268                CALL ZHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
00269      +                     A( N2 ), N )
00270                CALL ZTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N,
00271      +                     A( 0 ), N )
00272                CALL ZLAUUM( 'U', N2, A( N1 ), N, INFO )
00273 *
00274             END IF
00275 *
00276          ELSE
00277 *
00278 *           N is odd and TRANSR = 'C'
00279 *
00280             IF( LOWER ) THEN
00281 *
00282 *              SRPA for LOWER, TRANSPOSE, and N is odd
00283 *              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
00284 *
00285                CALL ZLAUUM( 'U', N1, A( 0 ), N1, INFO )
00286                CALL ZHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
00287      +                     A( 0 ), N1 )
00288                CALL ZTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1,
00289      +                     A( N1*N1 ), N1 )
00290                CALL ZLAUUM( 'L', N2, A( 1 ), N1, INFO )
00291 *
00292             ELSE
00293 *
00294 *              SRPA for UPPER, TRANSPOSE, and N is odd
00295 *              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
00296 *
00297                CALL ZLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
00298                CALL ZHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE,
00299      +                     A( N2*N2 ), N2 )
00300                CALL ZTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ),
00301      +                     N2, A( 0 ), N2 )
00302                CALL ZLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
00303 *
00304             END IF
00305 *
00306          END IF
00307 *
00308       ELSE
00309 *
00310 *        N is even
00311 *
00312          IF( NORMALTRANSR ) THEN
00313 *
00314 *           N is even and TRANSR = 'N'
00315 *
00316             IF( LOWER ) THEN
00317 *
00318 *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
00319 *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
00320 *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
00321 *
00322                CALL ZLAUUM( 'L', K, A( 1 ), N+1, INFO )
00323                CALL ZHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE,
00324      +                     A( 1 ), N+1 )
00325                CALL ZTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1,
00326      +                     A( K+1 ), N+1 )
00327                CALL ZLAUUM( 'U', K, A( 0 ), N+1, INFO )
00328 *
00329             ELSE
00330 *
00331 *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
00332 *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
00333 *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
00334 *
00335                CALL ZLAUUM( 'L', K, A( K+1 ), N+1, INFO )
00336                CALL ZHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
00337      +                     A( K+1 ), N+1 )
00338                CALL ZTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1,
00339      +                     A( 0 ), N+1 )
00340                CALL ZLAUUM( 'U', K, A( K ), N+1, INFO )
00341 *
00342             END IF
00343 *
00344          ELSE
00345 *
00346 *           N is even and TRANSR = 'C'
00347 *
00348             IF( LOWER ) THEN
00349 *
00350 *              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
00351 *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
00352 *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
00353 *
00354                CALL ZLAUUM( 'U', K, A( K ), K, INFO )
00355                CALL ZHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
00356      +                     A( K ), K )
00357                CALL ZTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K,
00358      +                     A( K*( K+1 ) ), K )
00359                CALL ZLAUUM( 'L', K, A( 0 ), K, INFO )
00360 *
00361             ELSE
00362 *
00363 *              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
00364 *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
00365 *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
00366 *
00367                CALL ZLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
00368                CALL ZHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE,
00369      +                     A( K*( K+1 ) ), K )
00370                CALL ZTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K,
00371      +                     A( 0 ), K )
00372                CALL ZLAUUM( 'L', K, A( K*K ), K, INFO )
00373 *
00374             END IF
00375 *
00376          END IF
00377 *
00378       END IF
00379 *
00380       RETURN
00381 *
00382 *     End of ZPFTRI
00383 *
00384       END
 All Files Functions