Annotation of rpl/lapack/lapack/dtptri.f, revision 1.17

1.8       bertrand    1: *> \brief \b DTPTRI
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.14      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.14      bertrand    9: *> Download DTPTRI + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtptri.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtptri.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtptri.f">
1.8       bertrand   15: *> [TXT]</a>
1.14      bertrand   16: *> \endhtmlonly
1.8       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DTPTRI( UPLO, DIAG, N, AP, INFO )
1.14      bertrand   22: *
1.8       bertrand   23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          DIAG, UPLO
                     25: *       INTEGER            INFO, N
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       DOUBLE PRECISION   AP( * )
                     29: *       ..
1.14      bertrand   30: *
1.8       bertrand   31: *
                     32: *> \par Purpose:
                     33: *  =============
                     34: *>
                     35: *> \verbatim
                     36: *>
                     37: *> DTPTRI computes the inverse of a real upper or lower triangular
                     38: *> matrix A stored in packed format.
                     39: *> \endverbatim
                     40: *
                     41: *  Arguments:
                     42: *  ==========
                     43: *
                     44: *> \param[in] UPLO
                     45: *> \verbatim
                     46: *>          UPLO is CHARACTER*1
                     47: *>          = 'U':  A is upper triangular;
                     48: *>          = 'L':  A is lower triangular.
                     49: *> \endverbatim
                     50: *>
                     51: *> \param[in] DIAG
                     52: *> \verbatim
                     53: *>          DIAG is CHARACTER*1
                     54: *>          = 'N':  A is non-unit triangular;
                     55: *>          = 'U':  A is unit triangular.
                     56: *> \endverbatim
                     57: *>
                     58: *> \param[in] N
                     59: *> \verbatim
                     60: *>          N is INTEGER
                     61: *>          The order of the matrix A.  N >= 0.
                     62: *> \endverbatim
                     63: *>
                     64: *> \param[in,out] AP
                     65: *> \verbatim
                     66: *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     67: *>          On entry, the upper or lower triangular matrix A, stored
                     68: *>          columnwise in a linear array.  The j-th column of A is stored
                     69: *>          in the array AP as follows:
                     70: *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
                     71: *>          if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
                     72: *>          See below for further details.
                     73: *>          On exit, the (triangular) inverse of the original matrix, in
                     74: *>          the same packed storage format.
                     75: *> \endverbatim
                     76: *>
                     77: *> \param[out] INFO
                     78: *> \verbatim
                     79: *>          INFO is INTEGER
                     80: *>          = 0:  successful exit
                     81: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                     82: *>          > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular
                     83: *>                matrix is singular and its inverse can not be computed.
                     84: *> \endverbatim
                     85: *
                     86: *  Authors:
                     87: *  ========
                     88: *
1.14      bertrand   89: *> \author Univ. of Tennessee
                     90: *> \author Univ. of California Berkeley
                     91: *> \author Univ. of Colorado Denver
                     92: *> \author NAG Ltd.
1.8       bertrand   93: *
                     94: *> \ingroup doubleOTHERcomputational
                     95: *
                     96: *> \par Further Details:
                     97: *  =====================
                     98: *>
                     99: *> \verbatim
                    100: *>
                    101: *>  A triangular matrix A can be transferred to packed storage using one
                    102: *>  of the following program segments:
                    103: *>
                    104: *>  UPLO = 'U':                      UPLO = 'L':
                    105: *>
                    106: *>        JC = 1                           JC = 1
                    107: *>        DO 2 J = 1, N                    DO 2 J = 1, N
                    108: *>           DO 1 I = 1, J                    DO 1 I = J, N
                    109: *>              AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)
                    110: *>      1    CONTINUE                    1    CONTINUE
                    111: *>           JC = JC + J                      JC = JC + N - J + 1
                    112: *>      2 CONTINUE                       2 CONTINUE
                    113: *> \endverbatim
                    114: *>
                    115: *  =====================================================================
1.1       bertrand  116:       SUBROUTINE DTPTRI( UPLO, DIAG, N, AP, INFO )
                    117: *
1.17    ! bertrand  118: *  -- LAPACK computational routine --
1.1       bertrand  119: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    120: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    121: *
                    122: *     .. Scalar Arguments ..
                    123:       CHARACTER          DIAG, UPLO
                    124:       INTEGER            INFO, N
                    125: *     ..
                    126: *     .. Array Arguments ..
                    127:       DOUBLE PRECISION   AP( * )
                    128: *     ..
                    129: *
                    130: *  =====================================================================
                    131: *
                    132: *     .. Parameters ..
                    133:       DOUBLE PRECISION   ONE, ZERO
                    134:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                    135: *     ..
                    136: *     .. Local Scalars ..
                    137:       LOGICAL            NOUNIT, UPPER
                    138:       INTEGER            J, JC, JCLAST, JJ
                    139:       DOUBLE PRECISION   AJJ
                    140: *     ..
                    141: *     .. External Functions ..
                    142:       LOGICAL            LSAME
                    143:       EXTERNAL           LSAME
                    144: *     ..
                    145: *     .. External Subroutines ..
                    146:       EXTERNAL           DSCAL, DTPMV, XERBLA
                    147: *     ..
                    148: *     .. Executable Statements ..
                    149: *
                    150: *     Test the input parameters.
                    151: *
                    152:       INFO = 0
                    153:       UPPER = LSAME( UPLO, 'U' )
                    154:       NOUNIT = LSAME( DIAG, 'N' )
                    155:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                    156:          INFO = -1
                    157:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
                    158:          INFO = -2
                    159:       ELSE IF( N.LT.0 ) THEN
                    160:          INFO = -3
                    161:       END IF
                    162:       IF( INFO.NE.0 ) THEN
                    163:          CALL XERBLA( 'DTPTRI', -INFO )
                    164:          RETURN
                    165:       END IF
                    166: *
                    167: *     Check for singularity if non-unit.
                    168: *
                    169:       IF( NOUNIT ) THEN
                    170:          IF( UPPER ) THEN
                    171:             JJ = 0
                    172:             DO 10 INFO = 1, N
                    173:                JJ = JJ + INFO
                    174:                IF( AP( JJ ).EQ.ZERO )
                    175:      $            RETURN
                    176:    10       CONTINUE
                    177:          ELSE
                    178:             JJ = 1
                    179:             DO 20 INFO = 1, N
                    180:                IF( AP( JJ ).EQ.ZERO )
                    181:      $            RETURN
                    182:                JJ = JJ + N - INFO + 1
                    183:    20       CONTINUE
                    184:          END IF
                    185:          INFO = 0
                    186:       END IF
                    187: *
                    188:       IF( UPPER ) THEN
                    189: *
                    190: *        Compute inverse of upper triangular matrix.
                    191: *
                    192:          JC = 1
                    193:          DO 30 J = 1, N
                    194:             IF( NOUNIT ) THEN
                    195:                AP( JC+J-1 ) = ONE / AP( JC+J-1 )
                    196:                AJJ = -AP( JC+J-1 )
                    197:             ELSE
                    198:                AJJ = -ONE
                    199:             END IF
                    200: *
                    201: *           Compute elements 1:j-1 of j-th column.
                    202: *
                    203:             CALL DTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
                    204:      $                  AP( JC ), 1 )
                    205:             CALL DSCAL( J-1, AJJ, AP( JC ), 1 )
                    206:             JC = JC + J
                    207:    30    CONTINUE
                    208: *
                    209:       ELSE
                    210: *
                    211: *        Compute inverse of lower triangular matrix.
                    212: *
                    213:          JC = N*( N+1 ) / 2
                    214:          DO 40 J = N, 1, -1
                    215:             IF( NOUNIT ) THEN
                    216:                AP( JC ) = ONE / AP( JC )
                    217:                AJJ = -AP( JC )
                    218:             ELSE
                    219:                AJJ = -ONE
                    220:             END IF
                    221:             IF( J.LT.N ) THEN
                    222: *
                    223: *              Compute elements j+1:n of j-th column.
                    224: *
                    225:                CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J,
                    226:      $                     AP( JCLAST ), AP( JC+1 ), 1 )
                    227:                CALL DSCAL( N-J, AJJ, AP( JC+1 ), 1 )
                    228:             END IF
                    229:             JCLAST = JC
                    230:             JC = JC - N + J - 2
                    231:    40    CONTINUE
                    232:       END IF
                    233: *
                    234:       RETURN
                    235: *
                    236: *     End of DTPTRI
                    237: *
                    238:       END

CVSweb interface <joel.bertrand@systella.fr>