LAPACK: SRC/dtftri.f Source File

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