1 SUBROUTINE CTRSV(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 COMPLEX A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * CTRSV solves one of the systems of equations
14 *
15 * A*x = b, or A**T*x = b, or A**H*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**H*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 - COMPLEX 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 - COMPLEX 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 COMPLEX ZERO
109 PARAMETER (ZERO= (0.0E+0,0.0E+0))
110 * ..
111 * .. Local Scalars ..
112 COMPLEX TEMP
113 INTEGER I,INFO,IX,J,JX,KX
114 LOGICAL NOCONJ,NOUNIT
115 * ..
116 * .. External Functions ..
117 LOGICAL LSAME
118 EXTERNAL LSAME
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL XERBLA
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC CONJG,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('CTRSV ',INFO)
146 RETURN
147 END IF
148 *
149 * Quick return if possible.
150 *
151 IF (N.EQ.0) RETURN
152 *
153 NOCONJ = LSAME(TRANS,'T')
154 NOUNIT = LSAME(DIAG,'N')
155 *
156 * Set up the start point in X if the increment is not unity. This
157 * will be ( N - 1 )*INCX too small for descending loops.
158 *
159 IF (INCX.LE.0) THEN
160 KX = 1 - (N-1)*INCX
161 ELSE IF (INCX.NE.1) THEN
162 KX = 1
163 END IF
164 *
165 * Start the operations. In this version the elements of A are
166 * accessed sequentially with one pass through A.
167 *
168 IF (LSAME(TRANS,'N')) THEN
169 *
170 * Form x := inv( A )*x.
171 *
172 IF (LSAME(UPLO,'U')) THEN
173 IF (INCX.EQ.1) THEN
174 DO 20 J = N,1,-1
175 IF (X(J).NE.ZERO) THEN
176 IF (NOUNIT) X(J) = X(J)/A(J,J)
177 TEMP = X(J)
178 DO 10 I = J - 1,1,-1
179 X(I) = X(I) - TEMP*A(I,J)
180 10 CONTINUE
181 END IF
182 20 CONTINUE
183 ELSE
184 JX = KX + (N-1)*INCX
185 DO 40 J = N,1,-1
186 IF (X(JX).NE.ZERO) THEN
187 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
188 TEMP = X(JX)
189 IX = JX
190 DO 30 I = J - 1,1,-1
191 IX = IX - INCX
192 X(IX) = X(IX) - TEMP*A(I,J)
193 30 CONTINUE
194 END IF
195 JX = JX - INCX
196 40 CONTINUE
197 END IF
198 ELSE
199 IF (INCX.EQ.1) THEN
200 DO 60 J = 1,N
201 IF (X(J).NE.ZERO) THEN
202 IF (NOUNIT) X(J) = X(J)/A(J,J)
203 TEMP = X(J)
204 DO 50 I = J + 1,N
205 X(I) = X(I) - TEMP*A(I,J)
206 50 CONTINUE
207 END IF
208 60 CONTINUE
209 ELSE
210 JX = KX
211 DO 80 J = 1,N
212 IF (X(JX).NE.ZERO) THEN
213 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
214 TEMP = X(JX)
215 IX = JX
216 DO 70 I = J + 1,N
217 IX = IX + INCX
218 X(IX) = X(IX) - TEMP*A(I,J)
219 70 CONTINUE
220 END IF
221 JX = JX + INCX
222 80 CONTINUE
223 END IF
224 END IF
225 ELSE
226 *
227 * Form x := inv( A**T )*x or x := inv( A**H )*x.
228 *
229 IF (LSAME(UPLO,'U')) THEN
230 IF (INCX.EQ.1) THEN
231 DO 110 J = 1,N
232 TEMP = X(J)
233 IF (NOCONJ) THEN
234 DO 90 I = 1,J - 1
235 TEMP = TEMP - A(I,J)*X(I)
236 90 CONTINUE
237 IF (NOUNIT) TEMP = TEMP/A(J,J)
238 ELSE
239 DO 100 I = 1,J - 1
240 TEMP = TEMP - CONJG(A(I,J))*X(I)
241 100 CONTINUE
242 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
243 END IF
244 X(J) = TEMP
245 110 CONTINUE
246 ELSE
247 JX = KX
248 DO 140 J = 1,N
249 IX = KX
250 TEMP = X(JX)
251 IF (NOCONJ) THEN
252 DO 120 I = 1,J - 1
253 TEMP = TEMP - A(I,J)*X(IX)
254 IX = IX + INCX
255 120 CONTINUE
256 IF (NOUNIT) TEMP = TEMP/A(J,J)
257 ELSE
258 DO 130 I = 1,J - 1
259 TEMP = TEMP - CONJG(A(I,J))*X(IX)
260 IX = IX + INCX
261 130 CONTINUE
262 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
263 END IF
264 X(JX) = TEMP
265 JX = JX + INCX
266 140 CONTINUE
267 END IF
268 ELSE
269 IF (INCX.EQ.1) THEN
270 DO 170 J = N,1,-1
271 TEMP = X(J)
272 IF (NOCONJ) THEN
273 DO 150 I = N,J + 1,-1
274 TEMP = TEMP - A(I,J)*X(I)
275 150 CONTINUE
276 IF (NOUNIT) TEMP = TEMP/A(J,J)
277 ELSE
278 DO 160 I = N,J + 1,-1
279 TEMP = TEMP - CONJG(A(I,J))*X(I)
280 160 CONTINUE
281 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
282 END IF
283 X(J) = TEMP
284 170 CONTINUE
285 ELSE
286 KX = KX + (N-1)*INCX
287 JX = KX
288 DO 200 J = N,1,-1
289 IX = KX
290 TEMP = X(JX)
291 IF (NOCONJ) THEN
292 DO 180 I = N,J + 1,-1
293 TEMP = TEMP - A(I,J)*X(IX)
294 IX = IX - INCX
295 180 CONTINUE
296 IF (NOUNIT) TEMP = TEMP/A(J,J)
297 ELSE
298 DO 190 I = N,J + 1,-1
299 TEMP = TEMP - CONJG(A(I,J))*X(IX)
300 IX = IX - INCX
301 190 CONTINUE
302 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
303 END IF
304 X(JX) = TEMP
305 JX = JX - INCX
306 200 CONTINUE
307 END IF
308 END IF
309 END IF
310 *
311 RETURN
312 *
313 * End of CTRSV .
314 *
315 END
2 * .. Scalar Arguments ..
3 INTEGER INCX,LDA,N
4 CHARACTER DIAG,TRANS,UPLO
5 * ..
6 * .. Array Arguments ..
7 COMPLEX A(LDA,*),X(*)
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * CTRSV solves one of the systems of equations
14 *
15 * A*x = b, or A**T*x = b, or A**H*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**H*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 - COMPLEX 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 - COMPLEX 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 COMPLEX ZERO
109 PARAMETER (ZERO= (0.0E+0,0.0E+0))
110 * ..
111 * .. Local Scalars ..
112 COMPLEX TEMP
113 INTEGER I,INFO,IX,J,JX,KX
114 LOGICAL NOCONJ,NOUNIT
115 * ..
116 * .. External Functions ..
117 LOGICAL LSAME
118 EXTERNAL LSAME
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL XERBLA
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC CONJG,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('CTRSV ',INFO)
146 RETURN
147 END IF
148 *
149 * Quick return if possible.
150 *
151 IF (N.EQ.0) RETURN
152 *
153 NOCONJ = LSAME(TRANS,'T')
154 NOUNIT = LSAME(DIAG,'N')
155 *
156 * Set up the start point in X if the increment is not unity. This
157 * will be ( N - 1 )*INCX too small for descending loops.
158 *
159 IF (INCX.LE.0) THEN
160 KX = 1 - (N-1)*INCX
161 ELSE IF (INCX.NE.1) THEN
162 KX = 1
163 END IF
164 *
165 * Start the operations. In this version the elements of A are
166 * accessed sequentially with one pass through A.
167 *
168 IF (LSAME(TRANS,'N')) THEN
169 *
170 * Form x := inv( A )*x.
171 *
172 IF (LSAME(UPLO,'U')) THEN
173 IF (INCX.EQ.1) THEN
174 DO 20 J = N,1,-1
175 IF (X(J).NE.ZERO) THEN
176 IF (NOUNIT) X(J) = X(J)/A(J,J)
177 TEMP = X(J)
178 DO 10 I = J - 1,1,-1
179 X(I) = X(I) - TEMP*A(I,J)
180 10 CONTINUE
181 END IF
182 20 CONTINUE
183 ELSE
184 JX = KX + (N-1)*INCX
185 DO 40 J = N,1,-1
186 IF (X(JX).NE.ZERO) THEN
187 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
188 TEMP = X(JX)
189 IX = JX
190 DO 30 I = J - 1,1,-1
191 IX = IX - INCX
192 X(IX) = X(IX) - TEMP*A(I,J)
193 30 CONTINUE
194 END IF
195 JX = JX - INCX
196 40 CONTINUE
197 END IF
198 ELSE
199 IF (INCX.EQ.1) THEN
200 DO 60 J = 1,N
201 IF (X(J).NE.ZERO) THEN
202 IF (NOUNIT) X(J) = X(J)/A(J,J)
203 TEMP = X(J)
204 DO 50 I = J + 1,N
205 X(I) = X(I) - TEMP*A(I,J)
206 50 CONTINUE
207 END IF
208 60 CONTINUE
209 ELSE
210 JX = KX
211 DO 80 J = 1,N
212 IF (X(JX).NE.ZERO) THEN
213 IF (NOUNIT) X(JX) = X(JX)/A(J,J)
214 TEMP = X(JX)
215 IX = JX
216 DO 70 I = J + 1,N
217 IX = IX + INCX
218 X(IX) = X(IX) - TEMP*A(I,J)
219 70 CONTINUE
220 END IF
221 JX = JX + INCX
222 80 CONTINUE
223 END IF
224 END IF
225 ELSE
226 *
227 * Form x := inv( A**T )*x or x := inv( A**H )*x.
228 *
229 IF (LSAME(UPLO,'U')) THEN
230 IF (INCX.EQ.1) THEN
231 DO 110 J = 1,N
232 TEMP = X(J)
233 IF (NOCONJ) THEN
234 DO 90 I = 1,J - 1
235 TEMP = TEMP - A(I,J)*X(I)
236 90 CONTINUE
237 IF (NOUNIT) TEMP = TEMP/A(J,J)
238 ELSE
239 DO 100 I = 1,J - 1
240 TEMP = TEMP - CONJG(A(I,J))*X(I)
241 100 CONTINUE
242 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
243 END IF
244 X(J) = TEMP
245 110 CONTINUE
246 ELSE
247 JX = KX
248 DO 140 J = 1,N
249 IX = KX
250 TEMP = X(JX)
251 IF (NOCONJ) THEN
252 DO 120 I = 1,J - 1
253 TEMP = TEMP - A(I,J)*X(IX)
254 IX = IX + INCX
255 120 CONTINUE
256 IF (NOUNIT) TEMP = TEMP/A(J,J)
257 ELSE
258 DO 130 I = 1,J - 1
259 TEMP = TEMP - CONJG(A(I,J))*X(IX)
260 IX = IX + INCX
261 130 CONTINUE
262 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
263 END IF
264 X(JX) = TEMP
265 JX = JX + INCX
266 140 CONTINUE
267 END IF
268 ELSE
269 IF (INCX.EQ.1) THEN
270 DO 170 J = N,1,-1
271 TEMP = X(J)
272 IF (NOCONJ) THEN
273 DO 150 I = N,J + 1,-1
274 TEMP = TEMP - A(I,J)*X(I)
275 150 CONTINUE
276 IF (NOUNIT) TEMP = TEMP/A(J,J)
277 ELSE
278 DO 160 I = N,J + 1,-1
279 TEMP = TEMP - CONJG(A(I,J))*X(I)
280 160 CONTINUE
281 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
282 END IF
283 X(J) = TEMP
284 170 CONTINUE
285 ELSE
286 KX = KX + (N-1)*INCX
287 JX = KX
288 DO 200 J = N,1,-1
289 IX = KX
290 TEMP = X(JX)
291 IF (NOCONJ) THEN
292 DO 180 I = N,J + 1,-1
293 TEMP = TEMP - A(I,J)*X(IX)
294 IX = IX - INCX
295 180 CONTINUE
296 IF (NOUNIT) TEMP = TEMP/A(J,J)
297 ELSE
298 DO 190 I = N,J + 1,-1
299 TEMP = TEMP - CONJG(A(I,J))*X(IX)
300 IX = IX - INCX
301 190 CONTINUE
302 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
303 END IF
304 X(JX) = TEMP
305 JX = JX - INCX
306 200 CONTINUE
307 END IF
308 END IF
309 END IF
310 *
311 RETURN
312 *
313 * End of CTRSV .
314 *
315 END