1
SUBROUTINE STPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
2
* .. Scalar Arguments ..
4
CHARACTER DIAG,TRANS,UPLO
6
* .. Array Arguments ..
13
* STPSV solves one of the systems of equations
15
* A*x = b, or A'*x = b,
17
* where b and x are n element vectors and A is an n by n unit, or
18
* non-unit, upper or lower triangular matrix, supplied in packed form.
20
* No test for singularity or near-singularity is included in this
21
* routine. Such tests must be performed before calling this routine.
27
* On entry, UPLO specifies whether the matrix is an upper or
28
* lower triangular matrix as follows:
30
* UPLO = 'U' or 'u' A is an upper triangular matrix.
32
* UPLO = 'L' or 'l' A is a lower triangular matrix.
36
* TRANS - CHARACTER*1.
37
* On entry, TRANS specifies the equations to be solved as
40
* TRANS = 'N' or 'n' A*x = b.
42
* TRANS = 'T' or 't' A'*x = b.
44
* TRANS = 'C' or 'c' A'*x = b.
49
* On entry, DIAG specifies whether or not A is unit
50
* triangular as follows:
52
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
54
* DIAG = 'N' or 'n' A is not assumed to be unit
60
* On entry, N specifies the order of the matrix A.
61
* N must be at least zero.
64
* AP - REAL array of DIMENSION at least
65
* ( ( n*( n + 1 ) )/2 ).
66
* Before entry with UPLO = 'U' or 'u', the array AP must
67
* contain the upper triangular matrix packed sequentially,
68
* column by column, so that AP( 1 ) contains a( 1, 1 ),
69
* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
70
* respectively, and so on.
71
* Before entry with UPLO = 'L' or 'l', the array AP must
72
* contain the lower triangular matrix packed sequentially,
73
* column by column, so that AP( 1 ) contains a( 1, 1 ),
74
* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
75
* respectively, and so on.
76
* Note that when DIAG = 'U' or 'u', the diagonal elements of
77
* A are not referenced, but are assumed to be unity.
80
* X - REAL array of dimension at least
81
* ( 1 + ( n - 1 )*abs( INCX ) ).
82
* Before entry, the incremented array X must contain the n
83
* element right-hand side vector b. On exit, X is overwritten
84
* with the solution vector x.
87
* On entry, INCX specifies the increment for the elements of
88
* X. INCX must not be zero.
92
* Level 2 Blas routine.
94
* -- Written on 22-October-1986.
95
* Jack Dongarra, Argonne National Lab.
96
* Jeremy Du Croz, Nag Central Office.
97
* Sven Hammarling, Nag Central Office.
98
* Richard Hanson, Sandia National Labs.
103
PARAMETER (ZERO=0.0E+0)
105
* .. Local Scalars ..
107
INTEGER I,INFO,IX,J,JX,K,KK,KX
110
* .. External Functions ..
114
* .. External Subroutines ..
118
* Test the input parameters.
121
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
123
ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
124
+ .NOT.LSAME(TRANS,'C')) THEN
126
ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
128
ELSE IF (N.LT.0) THEN
130
ELSE IF (INCX.EQ.0) THEN
134
CALL XERBLA('STPSV ',INFO)
138
* Quick return if possible.
142
NOUNIT = LSAME(DIAG,'N')
144
* Set up the start point in X if the increment is not unity. This
145
* will be ( N - 1 )*INCX too small for descending loops.
149
ELSE IF (INCX.NE.1) THEN
153
* Start the operations. In this version the elements of AP are
154
* accessed sequentially with one pass through AP.
156
IF (LSAME(TRANS,'N')) THEN
158
* Form x := inv( A )*x.
160
IF (LSAME(UPLO,'U')) THEN
164
IF (X(J).NE.ZERO) THEN
165
IF (NOUNIT) X(J) = X(J)/AP(KK)
169
X(I) = X(I) - TEMP*AP(K)
178
IF (X(JX).NE.ZERO) THEN
179
IF (NOUNIT) X(JX) = X(JX)/AP(KK)
182
DO 30 K = KK - 1,KK - J + 1,-1
184
X(IX) = X(IX) - TEMP*AP(K)
195
IF (X(J).NE.ZERO) THEN
196
IF (NOUNIT) X(J) = X(J)/AP(KK)
200
X(I) = X(I) - TEMP*AP(K)
209
IF (X(JX).NE.ZERO) THEN
210
IF (NOUNIT) X(JX) = X(JX)/AP(KK)
213
DO 70 K = KK + 1,KK + N - J
215
X(IX) = X(IX) - TEMP*AP(K)
225
* Form x := inv( A' )*x.
227
IF (LSAME(UPLO,'U')) THEN
234
TEMP = TEMP - AP(K)*X(I)
237
IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
246
DO 110 K = KK,KK + J - 2
247
TEMP = TEMP - AP(K)*X(IX)
250
IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
262
DO 130 I = N,J + 1,-1
263
TEMP = TEMP - AP(K)*X(I)
266
IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
276
DO 150 K = KK,KK - (N- (J+1)),-1
277
TEMP = TEMP - AP(K)*X(IX)
280
IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)