LAPACK 3.3.1 Linear Algebra PACKage

# cpftrf.f

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