|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 |
SUBROUTINE STPTRI( UPLO, DIAG, N, AP, INFO )
* * -- LAPACK routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER INFO, N * .. * .. Array Arguments .. REAL AP( * ) * .. * * Purpose * ======= * * STPTRI computes the inverse of a real upper or lower triangular * matrix A stored in packed format. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * = 'U': A is upper triangular; * = 'L': A is lower triangular. * * DIAG (input) CHARACTER*1 * = 'N': A is non-unit triangular; * = 'U': A is unit triangular. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * AP (input/output) REAL array, dimension (N*(N+1)/2) * On entry, the upper or lower triangular matrix A, stored * columnwise in a linear array. The j-th column of A is stored * in the array AP as follows: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; * if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n. * See below for further details. * On exit, the (triangular) inverse of the original matrix, in * the same packed storage format. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, A(i,i) is exactly zero. The triangular * matrix is singular and its inverse can not be computed. * * Further Details * =============== * * A triangular matrix A can be transferred to packed storage using one * of the following program segments: * * UPLO = 'U': UPLO = 'L': * * JC = 1 JC = 1 * DO 2 J = 1, N DO 2 J = 1, N * DO 1 I = 1, J DO 1 I = J, N * AP(JC+I-1) = A(I,J) AP(JC+I-J) = A(I,J) * 1 CONTINUE 1 CONTINUE * JC = JC + J JC = JC + N - J + 1 * 2 CONTINUE 2 CONTINUE * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER J, JC, JCLAST, JJ REAL AJJ * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL SSCAL, STPMV, XERBLA * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) NOUNIT = LSAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'STPTRI', -INFO ) RETURN END IF * * Check for singularity if non-unit. * IF( NOUNIT ) THEN IF( UPPER ) THEN JJ = 0 DO 10 INFO = 1, N JJ = JJ + INFO IF( AP( JJ ).EQ.ZERO ) $ RETURN 10 CONTINUE ELSE JJ = 1 DO 20 INFO = 1, N IF( AP( JJ ).EQ.ZERO ) $ RETURN JJ = JJ + N - INFO + 1 20 CONTINUE END IF INFO = 0 END IF * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix. * JC = 1 DO 30 J = 1, N IF( NOUNIT ) THEN AP( JC+J-1 ) = ONE / AP( JC+J-1 ) AJJ = -AP( JC+J-1 ) ELSE AJJ = -ONE END IF * * Compute elements 1:j-1 of j-th column. * CALL STPMV( 'Upper', 'No transpose', DIAG, J-1, AP, $ AP( JC ), 1 ) CALL SSCAL( J-1, AJJ, AP( JC ), 1 ) JC = JC + J 30 CONTINUE * ELSE * * Compute inverse of lower triangular matrix. * JC = N*( N+1 ) / 2 DO 40 J = N, 1, -1 IF( NOUNIT ) THEN AP( JC ) = ONE / AP( JC ) AJJ = -AP( JC ) ELSE AJJ = -ONE END IF IF( J.LT.N ) THEN * * Compute elements j+1:n of j-th column. * CALL STPMV( 'Lower', 'No transpose', DIAG, N-J, $ AP( JCLAST ), AP( JC+1 ), 1 ) CALL SSCAL( N-J, AJJ, AP( JC+1 ), 1 ) END IF JCLAST = JC JC = JC - N + J - 2 40 CONTINUE END IF * RETURN * * End of STPTRI * END |