LAPACK 3.3.1
Linear Algebra PACKage

ctptri.f

Go to the documentation of this file.
00001       SUBROUTINE CTPTRI( 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       COMPLEX            AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CTPTRI computes the inverse of a complex 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) COMPLEX 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       COMPLEX            ONE, ZERO
00072       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
00073      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00074 *     ..
00075 *     .. Local Scalars ..
00076       LOGICAL            NOUNIT, UPPER
00077       INTEGER            J, JC, JCLAST, JJ
00078       COMPLEX            AJJ
00079 *     ..
00080 *     .. External Functions ..
00081       LOGICAL            LSAME
00082       EXTERNAL           LSAME
00083 *     ..
00084 *     .. External Subroutines ..
00085       EXTERNAL           CSCAL, CTPMV, XERBLA
00086 *     ..
00087 *     .. Executable Statements ..
00088 *
00089 *     Test the input parameters.
00090 *
00091       INFO = 0
00092       UPPER = LSAME( UPLO, 'U' )
00093       NOUNIT = LSAME( DIAG, 'N' )
00094       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00095          INFO = -1
00096       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00097          INFO = -2
00098       ELSE IF( N.LT.0 ) THEN
00099          INFO = -3
00100       END IF
00101       IF( INFO.NE.0 ) THEN
00102          CALL XERBLA( 'CTPTRI', -INFO )
00103          RETURN
00104       END IF
00105 *
00106 *     Check for singularity if non-unit.
00107 *
00108       IF( NOUNIT ) THEN
00109          IF( UPPER ) THEN
00110             JJ = 0
00111             DO 10 INFO = 1, N
00112                JJ = JJ + INFO
00113                IF( AP( JJ ).EQ.ZERO )
00114      $            RETURN
00115    10       CONTINUE
00116          ELSE
00117             JJ = 1
00118             DO 20 INFO = 1, N
00119                IF( AP( JJ ).EQ.ZERO )
00120      $            RETURN
00121                JJ = JJ + N - INFO + 1
00122    20       CONTINUE
00123          END IF
00124          INFO = 0
00125       END IF
00126 *
00127       IF( UPPER ) THEN
00128 *
00129 *        Compute inverse of upper triangular matrix.
00130 *
00131          JC = 1
00132          DO 30 J = 1, N
00133             IF( NOUNIT ) THEN
00134                AP( JC+J-1 ) = ONE / AP( JC+J-1 )
00135                AJJ = -AP( JC+J-1 )
00136             ELSE
00137                AJJ = -ONE
00138             END IF
00139 *
00140 *           Compute elements 1:j-1 of j-th column.
00141 *
00142             CALL CTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
00143      $                  AP( JC ), 1 )
00144             CALL CSCAL( J-1, AJJ, AP( JC ), 1 )
00145             JC = JC + J
00146    30    CONTINUE
00147 *
00148       ELSE
00149 *
00150 *        Compute inverse of lower triangular matrix.
00151 *
00152          JC = N*( N+1 ) / 2
00153          DO 40 J = N, 1, -1
00154             IF( NOUNIT ) THEN
00155                AP( JC ) = ONE / AP( JC )
00156                AJJ = -AP( JC )
00157             ELSE
00158                AJJ = -ONE
00159             END IF
00160             IF( J.LT.N ) THEN
00161 *
00162 *              Compute elements j+1:n of j-th column.
00163 *
00164                CALL CTPMV( 'Lower', 'No transpose', DIAG, N-J,
00165      $                     AP( JCLAST ), AP( JC+1 ), 1 )
00166                CALL CSCAL( N-J, AJJ, AP( JC+1 ), 1 )
00167             END IF
00168             JCLAST = JC
00169             JC = JC - N + J - 2
00170    40    CONTINUE
00171       END IF
00172 *
00173       RETURN
00174 *
00175 *     End of CTPTRI
00176 *
00177       END
 All Files Functions