1 SUBROUTINE STRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
2 * .. Scalar Arguments ..
3 INTEGER INCX,LDA,N
4 CHARACTER DIAG,TRANS,UPLO
5 * ..
6 * .. Array Arguments ..
7 REAL A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * STRSV solves one of the systems of equations
14 *
15 * A*x = b, or A**T*x = b,
16 *
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.
19 *
20 * No test for singularity or near-singularity is included in this
21 * routine. Such tests must be performed before calling this routine.
22 *
23 * Arguments
24 * ==========
25 *
26 * UPLO - CHARACTER*1.
27 * On entry, UPLO specifies whether the matrix is an upper or
28 * lower triangular matrix as follows:
29 *
30 * UPLO = 'U' or 'u' A is an upper triangular matrix.
31 *
32 * UPLO = 'L' or 'l' A is a lower triangular matrix.
33 *
34 * Unchanged on exit.
35 *
36 * TRANS - CHARACTER*1.
37 * On entry, TRANS specifies the equations to be solved as
38 * follows:
39 *
40 * TRANS = 'N' or 'n' A*x = b.
41 *
42 * TRANS = 'T' or 't' A**T*x = b.
43 *
44 * TRANS = 'C' or 'c' A**T*x = b.
45 *
46 * Unchanged on exit.
47 *
48 * DIAG - CHARACTER*1.
49 * On entry, DIAG specifies whether or not A is unit
50 * triangular as follows:
51 *
52 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
53 *
54 * DIAG = 'N' or 'n' A is not assumed to be unit
55 * triangular.
56 *
57 * Unchanged on exit.
58 *
59 * N - INTEGER.
60 * On entry, N specifies the order of the matrix A.
61 * N must be at least zero.
62 * Unchanged on exit.
63 *
64 * A - REAL array of DIMENSION ( LDA, n ).
65 * Before entry with UPLO = 'U' or 'u', the leading n by n
66 * upper triangular part of the array A must contain the upper
67 * triangular matrix and the strictly lower triangular part of
68 * A is not referenced.
69 * Before entry with UPLO = 'L' or 'l', the leading n by n
70 * lower triangular part of the array A must contain the lower
71 * triangular matrix and the strictly upper triangular part of
72 * A is not referenced.
73 * Note that when DIAG = 'U' or 'u', the diagonal elements of
74 * A are not referenced either, but are assumed to be unity.
75 * Unchanged on exit.
76 *
77 * LDA - INTEGER.
78 * On entry, LDA specifies the first dimension of A as declared
79 * in the calling (sub) program. LDA must be at least
80 * max( 1, n ).
81 * Unchanged on exit.
82 *
83 * X - REAL array of dimension at least
84 * ( 1 + ( n - 1 )*abs( INCX ) ).
85 * Before entry, the incremented array X must contain the n
86 * element right-hand side vector b. On exit, X is overwritten
87 * with the solution vector x.
88 *
89 * INCX - INTEGER.
90 * On entry, INCX specifies the increment for the elements of
91 * X. INCX must not be zero.
92 * Unchanged on exit.
93 *
94 * Further Details
95 * ===============
96 *
97 * Level 2 Blas routine.
98 *
99 * -- Written on 22-October-1986.
100 * Jack Dongarra, Argonne National Lab.
101 * Jeremy Du Croz, Nag Central Office.
102 * Sven Hammarling, Nag Central Office.
103 * Richard Hanson, Sandia National Labs.
104 *
105 * =====================================================================
106 *
107 * .. Parameters ..
108 REAL ZERO
109 PARAMETER (ZERO=0.0E+0)
110 * ..
111 * .. Local Scalars ..
112 REAL TEMP
113 INTEGER I,INFO,IX,J,JX,KX
114 LOGICAL NOUNIT
115 * ..
116 * .. External Functions ..
117 LOGICAL LSAME
118 EXTERNAL LSAME
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL XERBLA
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC MAX
125 * ..
126 *
127 * Test the input parameters.
128 *
129 INFO = 0
130 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
131 INFO = 1
132 ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
133 + .NOT.LSAME(TRANS,'C')) THEN
134 INFO = 2
135 ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
136 INFO = 3
137 ELSE IF (N.LT.0) THEN
138 INFO = 4
139 ELSE IF (LDA.LT.MAX(1,N)) THEN
140 INFO = 6
141 ELSE IF (INCX.EQ.0) THEN
142 INFO = 8
143 END IF
144 IF (INFO.NE.0) THEN
145 CALL XERBLA('STRSV ',INFO)
146 RETURN
147 END IF
148 *
149 * Quick return if possible.
150 *
151 IF (N.EQ.0) RETURN
152 *
153 NOUNIT = LSAME(DIAG,'N')
154 *
155 * Set up the start point in X if the increment is not unity. This
156 * will be ( N - 1 )*INCX too small for descending loops.
157 *
158 IF (INCX.LE.0) THEN
159 KX = 1 - (N-1)*INCX
160 ELSE IF (INCX.NE.1) THEN
161 KX = 1
162 END IF
163 *
164 * Start the operations. In this version the elements of A are
165 * accessed sequentially with one pass through A.
166 *
167 IF (LSAME(TRANS,'N')) THEN
168 *
169 * Form x := inv( A )*x.
170 *
171 IF (LSAME(UPLO,'U')) THEN
172 IF (INCX.EQ.1) THEN
173 DO 20 J = N,1,-1
174 IF (X(J).NE.ZERO) THEN
175 IF (NOUNIT) X(J) = X(J)/A(J,J)
176 TEMP = X(J)
177 DO 10 I = J - 1,1,-1
178 X(I) = X(I) - TEMP*A(I,J)
179 10 CONTINUE
180 END IF
181 20 CONTINUE
182 ELSE
183 JX = KX + (N-1)*INCX
184 DO 40 J = N,1,-1
185 IF (X(JX).NE.ZERO) THEN
186 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
187 TEMP = X(JX)
188 IX = JX
189 DO 30 I = J - 1,1,-1
190 IX = IX - INCX
191 X(IX) = X(IX) - TEMP*A(I,J)
192 30 CONTINUE
193 END IF
194 JX = JX - INCX
195 40 CONTINUE
196 END IF
197 ELSE
198 IF (INCX.EQ.1) THEN
199 DO 60 J = 1,N
200 IF (X(J).NE.ZERO) THEN
201 IF (NOUNIT) X(J) = X(J)/A(J,J)
202 TEMP = X(J)
203 DO 50 I = J + 1,N
204 X(I) = X(I) - TEMP*A(I,J)
205 50 CONTINUE
206 END IF
207 60 CONTINUE
208 ELSE
209 JX = KX
210 DO 80 J = 1,N
211 IF (X(JX).NE.ZERO) THEN
212 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
213 TEMP = X(JX)
214 IX = JX
215 DO 70 I = J + 1,N
216 IX = IX + INCX
217 X(IX) = X(IX) - TEMP*A(I,J)
218 70 CONTINUE
219 END IF
220 JX = JX + INCX
221 80 CONTINUE
222 END IF
223 END IF
224 ELSE
225 *
226 * Form x := inv( A**T )*x.
227 *
228 IF (LSAME(UPLO,'U')) THEN
229 IF (INCX.EQ.1) THEN
230 DO 100 J = 1,N
231 TEMP = X(J)
232 DO 90 I = 1,J - 1
233 TEMP = TEMP - A(I,J)*X(I)
234 90 CONTINUE
235 IF (NOUNIT) TEMP = TEMP/A(J,J)
236 X(J) = TEMP
237 100 CONTINUE
238 ELSE
239 JX = KX
240 DO 120 J = 1,N
241 TEMP = X(JX)
242 IX = KX
243 DO 110 I = 1,J - 1
244 TEMP = TEMP - A(I,J)*X(IX)
245 IX = IX + INCX
246 110 CONTINUE
247 IF (NOUNIT) TEMP = TEMP/A(J,J)
248 X(JX) = TEMP
249 JX = JX + INCX
250 120 CONTINUE
251 END IF
252 ELSE
253 IF (INCX.EQ.1) THEN
254 DO 140 J = N,1,-1
255 TEMP = X(J)
256 DO 130 I = N,J + 1,-1
257 TEMP = TEMP - A(I,J)*X(I)
258 130 CONTINUE
259 IF (NOUNIT) TEMP = TEMP/A(J,J)
260 X(J) = TEMP
261 140 CONTINUE
262 ELSE
263 KX = KX + (N-1)*INCX
264 JX = KX
265 DO 160 J = N,1,-1
266 TEMP = X(JX)
267 IX = KX
268 DO 150 I = N,J + 1,-1
269 TEMP = TEMP - A(I,J)*X(IX)
270 IX = IX - INCX
271 150 CONTINUE
272 IF (NOUNIT) TEMP = TEMP/A(J,J)
273 X(JX) = TEMP
274 JX = JX - INCX
275 160 CONTINUE
276 END IF
277 END IF
278 END IF
279 *
280 RETURN
281 *
282 * End of STRSV .
283 *
284 END
2 * .. Scalar Arguments ..
3 INTEGER INCX,LDA,N
4 CHARACTER DIAG,TRANS,UPLO
5 * ..
6 * .. Array Arguments ..
7 REAL A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * STRSV solves one of the systems of equations
14 *
15 * A*x = b, or A**T*x = b,
16 *
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.
19 *
20 * No test for singularity or near-singularity is included in this
21 * routine. Such tests must be performed before calling this routine.
22 *
23 * Arguments
24 * ==========
25 *
26 * UPLO - CHARACTER*1.
27 * On entry, UPLO specifies whether the matrix is an upper or
28 * lower triangular matrix as follows:
29 *
30 * UPLO = 'U' or 'u' A is an upper triangular matrix.
31 *
32 * UPLO = 'L' or 'l' A is a lower triangular matrix.
33 *
34 * Unchanged on exit.
35 *
36 * TRANS - CHARACTER*1.
37 * On entry, TRANS specifies the equations to be solved as
38 * follows:
39 *
40 * TRANS = 'N' or 'n' A*x = b.
41 *
42 * TRANS = 'T' or 't' A**T*x = b.
43 *
44 * TRANS = 'C' or 'c' A**T*x = b.
45 *
46 * Unchanged on exit.
47 *
48 * DIAG - CHARACTER*1.
49 * On entry, DIAG specifies whether or not A is unit
50 * triangular as follows:
51 *
52 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
53 *
54 * DIAG = 'N' or 'n' A is not assumed to be unit
55 * triangular.
56 *
57 * Unchanged on exit.
58 *
59 * N - INTEGER.
60 * On entry, N specifies the order of the matrix A.
61 * N must be at least zero.
62 * Unchanged on exit.
63 *
64 * A - REAL array of DIMENSION ( LDA, n ).
65 * Before entry with UPLO = 'U' or 'u', the leading n by n
66 * upper triangular part of the array A must contain the upper
67 * triangular matrix and the strictly lower triangular part of
68 * A is not referenced.
69 * Before entry with UPLO = 'L' or 'l', the leading n by n
70 * lower triangular part of the array A must contain the lower
71 * triangular matrix and the strictly upper triangular part of
72 * A is not referenced.
73 * Note that when DIAG = 'U' or 'u', the diagonal elements of
74 * A are not referenced either, but are assumed to be unity.
75 * Unchanged on exit.
76 *
77 * LDA - INTEGER.
78 * On entry, LDA specifies the first dimension of A as declared
79 * in the calling (sub) program. LDA must be at least
80 * max( 1, n ).
81 * Unchanged on exit.
82 *
83 * X - REAL array of dimension at least
84 * ( 1 + ( n - 1 )*abs( INCX ) ).
85 * Before entry, the incremented array X must contain the n
86 * element right-hand side vector b. On exit, X is overwritten
87 * with the solution vector x.
88 *
89 * INCX - INTEGER.
90 * On entry, INCX specifies the increment for the elements of
91 * X. INCX must not be zero.
92 * Unchanged on exit.
93 *
94 * Further Details
95 * ===============
96 *
97 * Level 2 Blas routine.
98 *
99 * -- Written on 22-October-1986.
100 * Jack Dongarra, Argonne National Lab.
101 * Jeremy Du Croz, Nag Central Office.
102 * Sven Hammarling, Nag Central Office.
103 * Richard Hanson, Sandia National Labs.
104 *
105 * =====================================================================
106 *
107 * .. Parameters ..
108 REAL ZERO
109 PARAMETER (ZERO=0.0E+0)
110 * ..
111 * .. Local Scalars ..
112 REAL TEMP
113 INTEGER I,INFO,IX,J,JX,KX
114 LOGICAL NOUNIT
115 * ..
116 * .. External Functions ..
117 LOGICAL LSAME
118 EXTERNAL LSAME
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL XERBLA
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC MAX
125 * ..
126 *
127 * Test the input parameters.
128 *
129 INFO = 0
130 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
131 INFO = 1
132 ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
133 + .NOT.LSAME(TRANS,'C')) THEN
134 INFO = 2
135 ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
136 INFO = 3
137 ELSE IF (N.LT.0) THEN
138 INFO = 4
139 ELSE IF (LDA.LT.MAX(1,N)) THEN
140 INFO = 6
141 ELSE IF (INCX.EQ.0) THEN
142 INFO = 8
143 END IF
144 IF (INFO.NE.0) THEN
145 CALL XERBLA('STRSV ',INFO)
146 RETURN
147 END IF
148 *
149 * Quick return if possible.
150 *
151 IF (N.EQ.0) RETURN
152 *
153 NOUNIT = LSAME(DIAG,'N')
154 *
155 * Set up the start point in X if the increment is not unity. This
156 * will be ( N - 1 )*INCX too small for descending loops.
157 *
158 IF (INCX.LE.0) THEN
159 KX = 1 - (N-1)*INCX
160 ELSE IF (INCX.NE.1) THEN
161 KX = 1
162 END IF
163 *
164 * Start the operations. In this version the elements of A are
165 * accessed sequentially with one pass through A.
166 *
167 IF (LSAME(TRANS,'N')) THEN
168 *
169 * Form x := inv( A )*x.
170 *
171 IF (LSAME(UPLO,'U')) THEN
172 IF (INCX.EQ.1) THEN
173 DO 20 J = N,1,-1
174 IF (X(J).NE.ZERO) THEN
175 IF (NOUNIT) X(J) = X(J)/A(J,J)
176 TEMP = X(J)
177 DO 10 I = J - 1,1,-1
178 X(I) = X(I) - TEMP*A(I,J)
179 10 CONTINUE
180 END IF
181 20 CONTINUE
182 ELSE
183 JX = KX + (N-1)*INCX
184 DO 40 J = N,1,-1
185 IF (X(JX).NE.ZERO) THEN
186 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
187 TEMP = X(JX)
188 IX = JX
189 DO 30 I = J - 1,1,-1
190 IX = IX - INCX
191 X(IX) = X(IX) - TEMP*A(I,J)
192 30 CONTINUE
193 END IF
194 JX = JX - INCX
195 40 CONTINUE
196 END IF
197 ELSE
198 IF (INCX.EQ.1) THEN
199 DO 60 J = 1,N
200 IF (X(J).NE.ZERO) THEN
201 IF (NOUNIT) X(J) = X(J)/A(J,J)
202 TEMP = X(J)
203 DO 50 I = J + 1,N
204 X(I) = X(I) - TEMP*A(I,J)
205 50 CONTINUE
206 END IF
207 60 CONTINUE
208 ELSE
209 JX = KX
210 DO 80 J = 1,N
211 IF (X(JX).NE.ZERO) THEN
212 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
213 TEMP = X(JX)
214 IX = JX
215 DO 70 I = J + 1,N
216 IX = IX + INCX
217 X(IX) = X(IX) - TEMP*A(I,J)
218 70 CONTINUE
219 END IF
220 JX = JX + INCX
221 80 CONTINUE
222 END IF
223 END IF
224 ELSE
225 *
226 * Form x := inv( A**T )*x.
227 *
228 IF (LSAME(UPLO,'U')) THEN
229 IF (INCX.EQ.1) THEN
230 DO 100 J = 1,N
231 TEMP = X(J)
232 DO 90 I = 1,J - 1
233 TEMP = TEMP - A(I,J)*X(I)
234 90 CONTINUE
235 IF (NOUNIT) TEMP = TEMP/A(J,J)
236 X(J) = TEMP
237 100 CONTINUE
238 ELSE
239 JX = KX
240 DO 120 J = 1,N
241 TEMP = X(JX)
242 IX = KX
243 DO 110 I = 1,J - 1
244 TEMP = TEMP - A(I,J)*X(IX)
245 IX = IX + INCX
246 110 CONTINUE
247 IF (NOUNIT) TEMP = TEMP/A(J,J)
248 X(JX) = TEMP
249 JX = JX + INCX
250 120 CONTINUE
251 END IF
252 ELSE
253 IF (INCX.EQ.1) THEN
254 DO 140 J = N,1,-1
255 TEMP = X(J)
256 DO 130 I = N,J + 1,-1
257 TEMP = TEMP - A(I,J)*X(I)
258 130 CONTINUE
259 IF (NOUNIT) TEMP = TEMP/A(J,J)
260 X(J) = TEMP
261 140 CONTINUE
262 ELSE
263 KX = KX + (N-1)*INCX
264 JX = KX
265 DO 160 J = N,1,-1
266 TEMP = X(JX)
267 IX = KX
268 DO 150 I = N,J + 1,-1
269 TEMP = TEMP - A(I,J)*X(IX)
270 IX = IX - INCX
271 150 CONTINUE
272 IF (NOUNIT) TEMP = TEMP/A(J,J)
273 X(JX) = TEMP
274 JX = JX - INCX
275 160 CONTINUE
276 END IF
277 END IF
278 END IF
279 *
280 RETURN
281 *
282 * End of STRSV .
283 *
284 END