1 SUBROUTINE DTRSV(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 DOUBLE PRECISION A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * DTRSV 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 *
95 * Level 2 Blas routine.
96 *
97 * -- Written on 22-October-1986.
98 * Jack Dongarra, Argonne National Lab.
99 * Jeremy Du Croz, Nag Central Office.
100 * Sven Hammarling, Nag Central Office.
101 * Richard Hanson, Sandia National Labs.
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106 DOUBLE PRECISION ZERO
107 PARAMETER (ZERO=0.0D+0)
108 * ..
109 * .. Local Scalars ..
110 DOUBLE PRECISION TEMP
111 INTEGER I,INFO,IX,J,JX,KX
112 LOGICAL NOUNIT
113 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 EXTERNAL LSAME
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL XERBLA
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC MAX
123 * ..
124 *
125 * Test the input parameters.
126 *
127 INFO = 0
128 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
129 INFO = 1
130 ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
131 + .NOT.LSAME(TRANS,'C')) THEN
132 INFO = 2
133 ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
134 INFO = 3
135 ELSE IF (N.LT.0) THEN
136 INFO = 4
137 ELSE IF (LDA.LT.MAX(1,N)) THEN
138 INFO = 6
139 ELSE IF (INCX.EQ.0) THEN
140 INFO = 8
141 END IF
142 IF (INFO.NE.0) THEN
143 CALL XERBLA('DTRSV ',INFO)
144 RETURN
145 END IF
146 *
147 * Quick return if possible.
148 *
149 IF (N.EQ.0) RETURN
150 *
151 NOUNIT = LSAME(DIAG,'N')
152 *
153 * Set up the start point in X if the increment is not unity. This
154 * will be ( N - 1 )*INCX too small for descending loops.
155 *
156 IF (INCX.LE.0) THEN
157 KX = 1 - (N-1)*INCX
158 ELSE IF (INCX.NE.1) THEN
159 KX = 1
160 END IF
161 *
162 * Start the operations. In this version the elements of A are
163 * accessed sequentially with one pass through A.
164 *
165 IF (LSAME(TRANS,'N')) THEN
166 *
167 * Form x := inv( A )*x.
168 *
169 IF (LSAME(UPLO,'U')) THEN
170 IF (INCX.EQ.1) THEN
171 DO 20 J = N,1,-1
172 IF (X(J).NE.ZERO) THEN
173 IF (NOUNIT) X(J) = X(J)/A(J,J)
174 TEMP = X(J)
175 DO 10 I = J - 1,1,-1
176 X(I) = X(I) - TEMP*A(I,J)
177 10 CONTINUE
178 END IF
179 20 CONTINUE
180 ELSE
181 JX = KX + (N-1)*INCX
182 DO 40 J = N,1,-1
183 IF (X(JX).NE.ZERO) THEN
184 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
185 TEMP = X(JX)
186 IX = JX
187 DO 30 I = J - 1,1,-1
188 IX = IX - INCX
189 X(IX) = X(IX) - TEMP*A(I,J)
190 30 CONTINUE
191 END IF
192 JX = JX - INCX
193 40 CONTINUE
194 END IF
195 ELSE
196 IF (INCX.EQ.1) THEN
197 DO 60 J = 1,N
198 IF (X(J).NE.ZERO) THEN
199 IF (NOUNIT) X(J) = X(J)/A(J,J)
200 TEMP = X(J)
201 DO 50 I = J + 1,N
202 X(I) = X(I) - TEMP*A(I,J)
203 50 CONTINUE
204 END IF
205 60 CONTINUE
206 ELSE
207 JX = KX
208 DO 80 J = 1,N
209 IF (X(JX).NE.ZERO) THEN
210 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
211 TEMP = X(JX)
212 IX = JX
213 DO 70 I = J + 1,N
214 IX = IX + INCX
215 X(IX) = X(IX) - TEMP*A(I,J)
216 70 CONTINUE
217 END IF
218 JX = JX + INCX
219 80 CONTINUE
220 END IF
221 END IF
222 ELSE
223 *
224 * Form x := inv( A**T )*x.
225 *
226 IF (LSAME(UPLO,'U')) THEN
227 IF (INCX.EQ.1) THEN
228 DO 100 J = 1,N
229 TEMP = X(J)
230 DO 90 I = 1,J - 1
231 TEMP = TEMP - A(I,J)*X(I)
232 90 CONTINUE
233 IF (NOUNIT) TEMP = TEMP/A(J,J)
234 X(J) = TEMP
235 100 CONTINUE
236 ELSE
237 JX = KX
238 DO 120 J = 1,N
239 TEMP = X(JX)
240 IX = KX
241 DO 110 I = 1,J - 1
242 TEMP = TEMP - A(I,J)*X(IX)
243 IX = IX + INCX
244 110 CONTINUE
245 IF (NOUNIT) TEMP = TEMP/A(J,J)
246 X(JX) = TEMP
247 JX = JX + INCX
248 120 CONTINUE
249 END IF
250 ELSE
251 IF (INCX.EQ.1) THEN
252 DO 140 J = N,1,-1
253 TEMP = X(J)
254 DO 130 I = N,J + 1,-1
255 TEMP = TEMP - A(I,J)*X(I)
256 130 CONTINUE
257 IF (NOUNIT) TEMP = TEMP/A(J,J)
258 X(J) = TEMP
259 140 CONTINUE
260 ELSE
261 KX = KX + (N-1)*INCX
262 JX = KX
263 DO 160 J = N,1,-1
264 TEMP = X(JX)
265 IX = KX
266 DO 150 I = N,J + 1,-1
267 TEMP = TEMP - A(I,J)*X(IX)
268 IX = IX - INCX
269 150 CONTINUE
270 IF (NOUNIT) TEMP = TEMP/A(J,J)
271 X(JX) = TEMP
272 JX = JX - INCX
273 160 CONTINUE
274 END IF
275 END IF
276 END IF
277 *
278 RETURN
279 *
280 * End of DTRSV .
281 *
282 END
2 * .. Scalar Arguments ..
3 INTEGER INCX,LDA,N
4 CHARACTER DIAG,TRANS,UPLO
5 * ..
6 * .. Array Arguments ..
7 DOUBLE PRECISION A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * DTRSV 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 *
95 * Level 2 Blas routine.
96 *
97 * -- Written on 22-October-1986.
98 * Jack Dongarra, Argonne National Lab.
99 * Jeremy Du Croz, Nag Central Office.
100 * Sven Hammarling, Nag Central Office.
101 * Richard Hanson, Sandia National Labs.
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106 DOUBLE PRECISION ZERO
107 PARAMETER (ZERO=0.0D+0)
108 * ..
109 * .. Local Scalars ..
110 DOUBLE PRECISION TEMP
111 INTEGER I,INFO,IX,J,JX,KX
112 LOGICAL NOUNIT
113 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 EXTERNAL LSAME
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL XERBLA
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC MAX
123 * ..
124 *
125 * Test the input parameters.
126 *
127 INFO = 0
128 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
129 INFO = 1
130 ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
131 + .NOT.LSAME(TRANS,'C')) THEN
132 INFO = 2
133 ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
134 INFO = 3
135 ELSE IF (N.LT.0) THEN
136 INFO = 4
137 ELSE IF (LDA.LT.MAX(1,N)) THEN
138 INFO = 6
139 ELSE IF (INCX.EQ.0) THEN
140 INFO = 8
141 END IF
142 IF (INFO.NE.0) THEN
143 CALL XERBLA('DTRSV ',INFO)
144 RETURN
145 END IF
146 *
147 * Quick return if possible.
148 *
149 IF (N.EQ.0) RETURN
150 *
151 NOUNIT = LSAME(DIAG,'N')
152 *
153 * Set up the start point in X if the increment is not unity. This
154 * will be ( N - 1 )*INCX too small for descending loops.
155 *
156 IF (INCX.LE.0) THEN
157 KX = 1 - (N-1)*INCX
158 ELSE IF (INCX.NE.1) THEN
159 KX = 1
160 END IF
161 *
162 * Start the operations. In this version the elements of A are
163 * accessed sequentially with one pass through A.
164 *
165 IF (LSAME(TRANS,'N')) THEN
166 *
167 * Form x := inv( A )*x.
168 *
169 IF (LSAME(UPLO,'U')) THEN
170 IF (INCX.EQ.1) THEN
171 DO 20 J = N,1,-1
172 IF (X(J).NE.ZERO) THEN
173 IF (NOUNIT) X(J) = X(J)/A(J,J)
174 TEMP = X(J)
175 DO 10 I = J - 1,1,-1
176 X(I) = X(I) - TEMP*A(I,J)
177 10 CONTINUE
178 END IF
179 20 CONTINUE
180 ELSE
181 JX = KX + (N-1)*INCX
182 DO 40 J = N,1,-1
183 IF (X(JX).NE.ZERO) THEN
184 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
185 TEMP = X(JX)
186 IX = JX
187 DO 30 I = J - 1,1,-1
188 IX = IX - INCX
189 X(IX) = X(IX) - TEMP*A(I,J)
190 30 CONTINUE
191 END IF
192 JX = JX - INCX
193 40 CONTINUE
194 END IF
195 ELSE
196 IF (INCX.EQ.1) THEN
197 DO 60 J = 1,N
198 IF (X(J).NE.ZERO) THEN
199 IF (NOUNIT) X(J) = X(J)/A(J,J)
200 TEMP = X(J)
201 DO 50 I = J + 1,N
202 X(I) = X(I) - TEMP*A(I,J)
203 50 CONTINUE
204 END IF
205 60 CONTINUE
206 ELSE
207 JX = KX
208 DO 80 J = 1,N
209 IF (X(JX).NE.ZERO) THEN
210 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
211 TEMP = X(JX)
212 IX = JX
213 DO 70 I = J + 1,N
214 IX = IX + INCX
215 X(IX) = X(IX) - TEMP*A(I,J)
216 70 CONTINUE
217 END IF
218 JX = JX + INCX
219 80 CONTINUE
220 END IF
221 END IF
222 ELSE
223 *
224 * Form x := inv( A**T )*x.
225 *
226 IF (LSAME(UPLO,'U')) THEN
227 IF (INCX.EQ.1) THEN
228 DO 100 J = 1,N
229 TEMP = X(J)
230 DO 90 I = 1,J - 1
231 TEMP = TEMP - A(I,J)*X(I)
232 90 CONTINUE
233 IF (NOUNIT) TEMP = TEMP/A(J,J)
234 X(J) = TEMP
235 100 CONTINUE
236 ELSE
237 JX = KX
238 DO 120 J = 1,N
239 TEMP = X(JX)
240 IX = KX
241 DO 110 I = 1,J - 1
242 TEMP = TEMP - A(I,J)*X(IX)
243 IX = IX + INCX
244 110 CONTINUE
245 IF (NOUNIT) TEMP = TEMP/A(J,J)
246 X(JX) = TEMP
247 JX = JX + INCX
248 120 CONTINUE
249 END IF
250 ELSE
251 IF (INCX.EQ.1) THEN
252 DO 140 J = N,1,-1
253 TEMP = X(J)
254 DO 130 I = N,J + 1,-1
255 TEMP = TEMP - A(I,J)*X(I)
256 130 CONTINUE
257 IF (NOUNIT) TEMP = TEMP/A(J,J)
258 X(J) = TEMP
259 140 CONTINUE
260 ELSE
261 KX = KX + (N-1)*INCX
262 JX = KX
263 DO 160 J = N,1,-1
264 TEMP = X(JX)
265 IX = KX
266 DO 150 I = N,J + 1,-1
267 TEMP = TEMP - A(I,J)*X(IX)
268 IX = IX - INCX
269 150 CONTINUE
270 IF (NOUNIT) TEMP = TEMP/A(J,J)
271 X(JX) = TEMP
272 JX = JX - INCX
273 160 CONTINUE
274 END IF
275 END IF
276 END IF
277 *
278 RETURN
279 *
280 * End of DTRSV .
281 *
282 END