LAPACK 3.3.0

dtptri.f

Go to the documentation of this file.
00001       SUBROUTINE DTPTRI( UPLO, DIAG, N, AP, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          DIAG, UPLO
00010       INTEGER            INFO, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DTPTRI computes the inverse of a real upper or lower triangular
00020 *  matrix A stored in packed format.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  UPLO    (input) CHARACTER*1
00026 *          = 'U':  A is upper triangular;
00027 *          = 'L':  A is lower triangular.
00028 *
00029 *  DIAG    (input) CHARACTER*1
00030 *          = 'N':  A is non-unit triangular;
00031 *          = 'U':  A is unit triangular.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00037 *          On entry, the upper or lower triangular matrix A, stored
00038 *          columnwise in a linear array.  The j-th column of A is stored
00039 *          in the array AP as follows:
00040 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00041 *          if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
00042 *          See below for further details.
00043 *          On exit, the (triangular) inverse of the original matrix, in
00044 *          the same packed storage format.
00045 *
00046 *  INFO    (output) INTEGER
00047 *          = 0:  successful exit
00048 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00049 *          > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular
00050 *                matrix is singular and its inverse can not be computed.
00051 *
00052 *  Further Details
00053 *  ===============
00054 *
00055 *  A triangular matrix A can be transferred to packed storage using one
00056 *  of the following program segments:
00057 *
00058 *  UPLO = 'U':                      UPLO = 'L':
00059 *
00060 *        JC = 1                           JC = 1
00061 *        DO 2 J = 1, N                    DO 2 J = 1, N
00062 *           DO 1 I = 1, J                    DO 1 I = J, N
00063 *              AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)
00064 *      1    CONTINUE                    1    CONTINUE
00065 *           JC = JC + J                      JC = JC + N - J + 1
00066 *      2 CONTINUE                       2 CONTINUE
00067 *
00068 *  =====================================================================
00069 *
00070 *     .. Parameters ..
00071       DOUBLE PRECISION   ONE, ZERO
00072       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00073 *     ..
00074 *     .. Local Scalars ..
00075       LOGICAL            NOUNIT, UPPER
00076       INTEGER            J, JC, JCLAST, JJ
00077       DOUBLE PRECISION   AJJ
00078 *     ..
00079 *     .. External Functions ..
00080       LOGICAL            LSAME
00081       EXTERNAL           LSAME
00082 *     ..
00083 *     .. External Subroutines ..
00084       EXTERNAL           DSCAL, DTPMV, XERBLA
00085 *     ..
00086 *     .. Executable Statements ..
00087 *
00088 *     Test the input parameters.
00089 *
00090       INFO = 0
00091       UPPER = LSAME( UPLO, 'U' )
00092       NOUNIT = LSAME( DIAG, 'N' )
00093       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00094          INFO = -1
00095       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00096          INFO = -2
00097       ELSE IF( N.LT.0 ) THEN
00098          INFO = -3
00099       END IF
00100       IF( INFO.NE.0 ) THEN
00101          CALL XERBLA( 'DTPTRI', -INFO )
00102          RETURN
00103       END IF
00104 *
00105 *     Check for singularity if non-unit.
00106 *
00107       IF( NOUNIT ) THEN
00108          IF( UPPER ) THEN
00109             JJ = 0
00110             DO 10 INFO = 1, N
00111                JJ = JJ + INFO
00112                IF( AP( JJ ).EQ.ZERO )
00113      $            RETURN
00114    10       CONTINUE
00115          ELSE
00116             JJ = 1
00117             DO 20 INFO = 1, N
00118                IF( AP( JJ ).EQ.ZERO )
00119      $            RETURN
00120                JJ = JJ + N - INFO + 1
00121    20       CONTINUE
00122          END IF
00123          INFO = 0
00124       END IF
00125 *
00126       IF( UPPER ) THEN
00127 *
00128 *        Compute inverse of upper triangular matrix.
00129 *
00130          JC = 1
00131          DO 30 J = 1, N
00132             IF( NOUNIT ) THEN
00133                AP( JC+J-1 ) = ONE / AP( JC+J-1 )
00134                AJJ = -AP( JC+J-1 )
00135             ELSE
00136                AJJ = -ONE
00137             END IF
00138 *
00139 *           Compute elements 1:j-1 of j-th column.
00140 *
00141             CALL DTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
00142      $                  AP( JC ), 1 )
00143             CALL DSCAL( J-1, AJJ, AP( JC ), 1 )
00144             JC = JC + J
00145    30    CONTINUE
00146 *
00147       ELSE
00148 *
00149 *        Compute inverse of lower triangular matrix.
00150 *
00151          JC = N*( N+1 ) / 2
00152          DO 40 J = N, 1, -1
00153             IF( NOUNIT ) THEN
00154                AP( JC ) = ONE / AP( JC )
00155                AJJ = -AP( JC )
00156             ELSE
00157                AJJ = -ONE
00158             END IF
00159             IF( J.LT.N ) THEN
00160 *
00161 *              Compute elements j+1:n of j-th column.
00162 *
00163                CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J,
00164      $                     AP( JCLAST ), AP( JC+1 ), 1 )
00165                CALL DSCAL( N-J, AJJ, AP( JC+1 ), 1 )
00166             END IF
00167             JCLAST = JC
00168             JC = JC - N + J - 2
00169    40    CONTINUE
00170       END IF
00171 *
00172       RETURN
00173 *
00174 *     End of DTPTRI
00175 *
00176       END
 All Files Functions