SRC/dtftri.f

   1 *> \brief \b DTFTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTFTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtftri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtftri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtftri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          TRANSR, UPLO, DIAG
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   A( 0: * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DTFTRI computes the inverse of a triangular matrix A stored in RFP
  38 *> format.
  39 *>
  40 *> This is a Level 3 BLAS version of the algorithm.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] TRANSR
  47 *> \verbatim
  48 *>          TRANSR is CHARACTER*1
  49 *>          = 'N':  The Normal TRANSR of RFP A is stored;
  50 *>          = 'T':  The Transpose TRANSR of RFP A is stored.
  51 *> \endverbatim
  52 *>
  53 *> \param[in] UPLO
  54 *> \verbatim
  55 *>          UPLO is CHARACTER*1
  56 *>          = 'U':  A is upper triangular;
  57 *>          = 'L':  A is lower triangular.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] DIAG
  61 *> \verbatim
  62 *>          DIAG is CHARACTER*1
  63 *>          = 'N':  A is non-unit triangular;
  64 *>          = 'U':  A is unit triangular.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The order of the matrix A.  N >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in,out] A
  74 *> \verbatim
  75 *>          A is DOUBLE PRECISION array, dimension (0:nt-1);
  76 *>          nt=N*(N+1)/2. On entry, the triangular factor of a Hermitian
  77 *>          Positive Definite matrix A in RFP format. RFP format is
  78 *>          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
  79 *>          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
  80 *>          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
  81 *>          the transpose of RFP A as defined when
  82 *>          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
  83 *>          follows: If UPLO = 'U' the RFP A contains the nt elements of
  84 *>          upper packed A; If UPLO = 'L' the RFP A contains the nt
  85 *>          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
  86 *>          TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
  87 *>          even and N is odd. See the Note below for more details.
  88 *>
  89 *>          On exit, the (triangular) inverse of the original matrix, in
  90 *>          the same storage format.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] INFO
  94 *> \verbatim
  95 *>          INFO is INTEGER
  96 *>          = 0: successful exit
  97 *>          < 0: if INFO = -i, the i-th argument had an illegal value
  98 *>          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
  99 *>               matrix is singular and its inverse can not be computed.
 100 *> \endverbatim
 101 *
 102 *  Authors:
 103 *  ========
 104 *
 105 *> \author Univ. of Tennessee
 106 *> \author Univ. of California Berkeley
 107 *> \author Univ. of Colorado Denver
 108 *> \author NAG Ltd.
 109 *
 110 *> \date November 2011
 111 *
 112 *> \ingroup doubleOTHERcomputational
 113 *
 114 *> \par Further Details:
 115 *  =====================
 116 *>
 117 *> \verbatim
 118 *>
 119 *>  We first consider Rectangular Full Packed (RFP) Format when N is
 120 *>  even. We give an example where N = 6.
 121 *>
 122 *>      AP is Upper             AP is Lower
 123 *>
 124 *>   00 01 02 03 04 05       00
 125 *>      11 12 13 14 15       10 11
 126 *>         22 23 24 25       20 21 22
 127 *>            33 34 35       30 31 32 33
 128 *>               44 45       40 41 42 43 44
 129 *>                  55       50 51 52 53 54 55
 130 *>
 131 *>
 132 *>  Let TRANSR = 'N'. RFP holds AP as follows:
 133 *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
 134 *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
 135 *>  the transpose of the first three columns of AP upper.
 136 *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
 137 *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
 138 *>  the transpose of the last three columns of AP lower.
 139 *>  This covers the case N even and TRANSR = 'N'.
 140 *>
 141 *>         RFP A                   RFP A
 142 *>
 143 *>        03 04 05                33 43 53
 144 *>        13 14 15                00 44 54
 145 *>        23 24 25                10 11 55
 146 *>        33 34 35                20 21 22
 147 *>        00 44 45                30 31 32
 148 *>        01 11 55                40 41 42
 149 *>        02 12 22                50 51 52
 150 *>
 151 *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
 152 *>  transpose of RFP A above. One therefore gets:
 153 *>
 154 *>
 155 *>           RFP A                   RFP A
 156 *>
 157 *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
 158 *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
 159 *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
 160 *>
 161 *>
 162 *>  We then consider Rectangular Full Packed (RFP) Format when N is
 163 *>  odd. We give an example where N = 5.
 164 *>
 165 *>     AP is Upper                 AP is Lower
 166 *>
 167 *>   00 01 02 03 04              00
 168 *>      11 12 13 14              10 11
 169 *>         22 23 24              20 21 22
 170 *>            33 34              30 31 32 33
 171 *>               44              40 41 42 43 44
 172 *>
 173 *>
 174 *>  Let TRANSR = 'N'. RFP holds AP as follows:
 175 *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
 176 *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
 177 *>  the transpose of the first two columns of AP upper.
 178 *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
 179 *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
 180 *>  the transpose of the last two columns of AP lower.
 181 *>  This covers the case N odd and TRANSR = 'N'.
 182 *>
 183 *>         RFP A                   RFP A
 184 *>
 185 *>        02 03 04                00 33 43
 186 *>        12 13 14                10 11 44
 187 *>        22 23 24                20 21 22
 188 *>        00 33 34                30 31 32
 189 *>        01 11 44                40 41 42
 190 *>
 191 *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
 192 *>  transpose of RFP A above. One therefore gets:
 193 *>
 194 *>           RFP A                   RFP A
 195 *>
 196 *>     02 12 22 00 01             00 10 20 30 40 50
 197 *>     03 13 23 33 11             33 11 21 31 41 51
 198 *>     04 14 24 34 44             43 44 22 32 42 52
 199 *> \endverbatim
 200 *>
 201 *  =====================================================================
 202       SUBROUTINE DTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
 203 *
 204 *  -- LAPACK computational routine (version 3.4.0) --
 205 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 206 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 207 *     November 2011
 208 *
 209 *     .. Scalar Arguments ..
 210       CHARACTER          TRANSR, UPLO, DIAG
 211       INTEGER            INFO, N
 212 *     ..
 213 *     .. Array Arguments ..
 214       DOUBLE PRECISION   A( 0: * )
 215 *     ..
 216 *
 217 *  =====================================================================
 218 *
 219 *     .. Parameters ..
 220       DOUBLE PRECISION   ONE
 221       PARAMETER          ( ONE = 1.0D+0 )
 222 *     ..
 223 *     .. Local Scalars ..
 224       LOGICAL            LOWER, NISODD, NORMALTRANSR
 225       INTEGER            N1, N2, K
 226 *     ..
 227 *     .. External Functions ..
 228       LOGICAL            LSAME
 229       EXTERNAL           LSAME
 230 *     ..
 231 *     .. External Subroutines ..
 232       EXTERNAL           XERBLA, DTRMM, DTRTRI
 233 *     ..
 234 *     .. Intrinsic Functions ..
 235       INTRINSIC          MOD
 236 *     ..
 237 *     .. Executable Statements ..
 238 *
 239 *     Test the input parameters.
 240 *
 241       INFO = 0
 242       NORMALTRANSR = LSAME( TRANSR, 'N' )
 243       LOWER = LSAME( UPLO, 'L' )
 244       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
 245          INFO = -1
 246       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
 247          INFO = -2
 248       ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
 249      $         THEN
 250          INFO = -3
 251       ELSE IF( N.LT.0 ) THEN
 252          INFO = -4
 253       END IF
 254       IF( INFO.NE.0 ) THEN
 255          CALL XERBLA( 'DTFTRI', -INFO )
 256          RETURN
 257       END IF
 258 *
 259 *     Quick return if possible
 260 *
 261       IF( N.EQ.0 )
 262      $   RETURN
 263 *
 264 *     If N is odd, set NISODD = .TRUE.
 265 *     If N is even, set K = N/2 and NISODD = .FALSE.
 266 *
 267       IF( MOD( N, 2 ).EQ.0 ) THEN
 268          K = N / 2
 269          NISODD = .FALSE.
 270       ELSE
 271          NISODD = .TRUE.
 272       END IF
 273 *
 274 *     Set N1 and N2 depending on LOWER
 275 *
 276       IF( LOWER ) THEN
 277          N2 = N / 2
 278          N1 = N - N2
 279       ELSE
 280          N1 = N / 2
 281          N2 = N - N1
 282       END IF
 283 *
 284 *
 285 *     start execution: there are eight cases
 286 *
 287       IF( NISODD ) THEN
 288 *
 289 *        N is odd
 290 *
 291          IF( NORMALTRANSR ) THEN
 292 *
 293 *           N is odd and TRANSR = 'N'
 294 *
 295             IF( LOWER ) THEN
 296 *
 297 *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
 298 *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
 299 *             T1 -> a(0), T2 -> a(n), S -> a(n1)
 300 *
 301                CALL DTRTRI( 'L', DIAG, N1, A( 0 ), N, INFO )
 302                IF( INFO.GT.0 )
 303      $            RETURN
 304                CALL DTRMM( 'R', 'L', 'N', DIAG, N2, N1, -ONE, A( 0 ),
 305      $                     N, A( N1 ), N )
 306                CALL DTRTRI( 'U', DIAG, N2, A( N ), N, INFO )
 307                IF( INFO.GT.0 )
 308      $            INFO = INFO + N1
 309                IF( INFO.GT.0 )
 310      $            RETURN
 311                CALL DTRMM( 'L', 'U', 'T', DIAG, N2, N1, ONE, A( N ), N,
 312      $                     A( N1 ), N )
 313 *
 314             ELSE
 315 *
 316 *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
 317 *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
 318 *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
 319 *
 320                CALL DTRTRI( 'L', DIAG, N1, A( N2 ), N, INFO )
 321                IF( INFO.GT.0 )
 322      $            RETURN
 323                CALL DTRMM( 'L', 'L', 'T', DIAG, N1, N2, -ONE, A( N2 ),
 324      $                     N, A( 0 ), N )
 325                CALL DTRTRI( 'U', DIAG, N2, A( N1 ), N, INFO )
 326                IF( INFO.GT.0 )
 327      $            INFO = INFO + N1
 328                IF( INFO.GT.0 )
 329      $            RETURN
 330                CALL DTRMM( 'R', 'U', 'N', DIAG, N1, N2, ONE, A( N1 ),
 331      $                     N, A( 0 ), N )
 332 *
 333             END IF
 334 *
 335          ELSE
 336 *
 337 *           N is odd and TRANSR = 'T'
 338 *
 339             IF( LOWER ) THEN
 340 *
 341 *              SRPA for LOWER, TRANSPOSE and N is odd
 342 *              T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
 343 *
 344                CALL DTRTRI( 'U', DIAG, N1, A( 0 ), N1, INFO )
 345                IF( INFO.GT.0 )
 346      $            RETURN
 347                CALL DTRMM( 'L', 'U', 'N', DIAG, N1, N2, -ONE, A( 0 ),
 348      $                     N1, A( N1*N1 ), N1 )
 349                CALL DTRTRI( 'L', DIAG, N2, A( 1 ), N1, INFO )
 350                IF( INFO.GT.0 )
 351      $            INFO = INFO + N1
 352                IF( INFO.GT.0 )
 353      $            RETURN
 354                CALL DTRMM( 'R', 'L', 'T', DIAG, N1, N2, ONE, A( 1 ),
 355      $                     N1, A( N1*N1 ), N1 )
 356 *
 357             ELSE
 358 *
 359 *              SRPA for UPPER, TRANSPOSE and N is odd
 360 *              T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
 361 *
 362                CALL DTRTRI( 'U', DIAG, N1, A( N2*N2 ), N2, INFO )
 363                IF( INFO.GT.0 )
 364      $            RETURN
 365                CALL DTRMM( 'R', 'U', 'T', DIAG, N2, N1, -ONE,
 366      $                     A( N2*N2 ), N2, A( 0 ), N2 )
 367                CALL DTRTRI( 'L', DIAG, N2, A( N1*N2 ), N2, INFO )
 368                IF( INFO.GT.0 )
 369      $            INFO = INFO + N1
 370                IF( INFO.GT.0 )
 371      $            RETURN
 372                CALL DTRMM( 'L', 'L', 'N', DIAG, N2, N1, ONE,
 373      $                     A( N1*N2 ), N2, A( 0 ), N2 )
 374             END IF
 375 *
 376          END IF
 377 *
 378       ELSE
 379 *
 380 *        N is even
 381 *
 382          IF( NORMALTRANSR ) THEN
 383 *
 384 *           N is even and TRANSR = 'N'
 385 *
 386             IF( LOWER ) THEN
 387 *
 388 *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
 389 *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
 390 *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
 391 *
 392                CALL DTRTRI( 'L', DIAG, K, A( 1 ), N+1, INFO )
 393                IF( INFO.GT.0 )
 394      $            RETURN
 395                CALL DTRMM( 'R', 'L', 'N', DIAG, K, K, -ONE, A( 1 ),
 396      $                     N+1, A( K+1 ), N+1 )
 397                CALL DTRTRI( 'U', DIAG, K, A( 0 ), N+1, INFO )
 398                IF( INFO.GT.0 )
 399      $            INFO = INFO + K
 400                IF( INFO.GT.0 )
 401      $            RETURN
 402                CALL DTRMM( 'L', 'U', 'T', DIAG, K, K, ONE, A( 0 ), N+1,
 403      $                     A( K+1 ), N+1 )
 404 *
 405             ELSE
 406 *
 407 *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
 408 *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
 409 *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
 410 *
 411                CALL DTRTRI( 'L', DIAG, K, A( K+1 ), N+1, INFO )
 412                IF( INFO.GT.0 )
 413      $            RETURN
 414                CALL DTRMM( 'L', 'L', 'T', DIAG, K, K, -ONE, A( K+1 ),
 415      $                     N+1, A( 0 ), N+1 )
 416                CALL DTRTRI( 'U', DIAG, K, A( K ), N+1, INFO )
 417                IF( INFO.GT.0 )
 418      $            INFO = INFO + K
 419                IF( INFO.GT.0 )
 420      $            RETURN
 421                CALL DTRMM( 'R', 'U', 'N', DIAG, K, K, ONE, A( K ), N+1,
 422      $                     A( 0 ), N+1 )
 423             END IF
 424          ELSE
 425 *
 426 *           N is even and TRANSR = 'T'
 427 *
 428             IF( LOWER ) THEN
 429 *
 430 *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
 431 *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
 432 *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
 433 *
 434                CALL DTRTRI( 'U', DIAG, K, A( K ), K, INFO )
 435                IF( INFO.GT.0 )
 436      $            RETURN
 437                CALL DTRMM( 'L', 'U', 'N', DIAG, K, K, -ONE, A( K ), K,
 438      $                     A( K*( K+1 ) ), K )
 439                CALL DTRTRI( 'L', DIAG, K, A( 0 ), K, INFO )
 440                IF( INFO.GT.0 )
 441      $            INFO = INFO + K
 442                IF( INFO.GT.0 )
 443      $            RETURN
 444                CALL DTRMM( 'R', 'L', 'T', DIAG, K, K, ONE, A( 0 ), K,
 445      $                     A( K*( K+1 ) ), K )
 446             ELSE
 447 *
 448 *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
 449 *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
 450 *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
 451 *
 452                CALL DTRTRI( 'U', DIAG, K, A( K*( K+1 ) ), K, INFO )
 453                IF( INFO.GT.0 )
 454      $            RETURN
 455                CALL DTRMM( 'R', 'U', 'T', DIAG, K, K, -ONE,
 456      $                     A( K*( K+1 ) ), K, A( 0 ), K )
 457                CALL DTRTRI( 'L', DIAG, K, A( K*K ), K, INFO )
 458                IF( INFO.GT.0 )
 459      $            INFO = INFO + K
 460                IF( INFO.GT.0 )
 461      $            RETURN
 462                CALL DTRMM( 'L', 'L', 'N', DIAG, K, K, ONE, A( K*K ), K,
 463      $                     A( 0 ), K )
 464             END IF
 465          END IF
 466       END IF
 467 *
 468       RETURN
 469 *
 470 *     End of DTFTRI
 471 *
 472       END