Mercurial > octave-nkf

comparison libcruft/blas/ctrsv.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
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'*x = b, or conjg( A' )*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'*x = b.
43 *
44 * TRANS = 'C' or 'c' conjg( A' )*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 *
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 * .. Parameters ..
105 COMPLEX ZERO
106 PARAMETER (ZERO= (0.0E+0,0.0E+0))
107 * ..
108 * .. Local Scalars ..
109 COMPLEX TEMP
110 INTEGER I,INFO,IX,J,JX,KX
111 LOGICAL NOCONJ,NOUNIT
112 * ..
113 * .. External Functions ..
114 LOGICAL LSAME
115 EXTERNAL LSAME
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL XERBLA
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC CONJG,MAX
122 * ..
123 *
124 * Test the input parameters.
125 *
126 INFO = 0
127 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
128 INFO = 1
129 ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
130 + .NOT.LSAME(TRANS,'C')) THEN
131 INFO = 2
132 ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
133 INFO = 3
134 ELSE IF (N.LT.0) THEN
135 INFO = 4
136 ELSE IF (LDA.LT.MAX(1,N)) THEN
137 INFO = 6
138 ELSE IF (INCX.EQ.0) THEN
139 INFO = 8
140 END IF
141 IF (INFO.NE.0) THEN
142 CALL XERBLA('CTRSV ',INFO)
143 RETURN
144 END IF
145 *
146 * Quick return if possible.
147 *
148 IF (N.EQ.0) RETURN
149 *
150 NOCONJ = LSAME(TRANS,'T')
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' )*x or x := inv( conjg( A' ) )*x.
225 *
226 IF (LSAME(UPLO,'U')) THEN
227 IF (INCX.EQ.1) THEN
228 DO 110 J = 1,N
229 TEMP = X(J)
230 IF (NOCONJ) THEN
231 DO 90 I = 1,J - 1
232 TEMP = TEMP - A(I,J)*X(I)
233 90 CONTINUE
234 IF (NOUNIT) TEMP = TEMP/A(J,J)
235 ELSE
236 DO 100 I = 1,J - 1
237 TEMP = TEMP - CONJG(A(I,J))*X(I)
238 100 CONTINUE
239 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
240 END IF
241 X(J) = TEMP
242 110 CONTINUE
243 ELSE
244 JX = KX
245 DO 140 J = 1,N
246 IX = KX
247 TEMP = X(JX)
248 IF (NOCONJ) THEN
249 DO 120 I = 1,J - 1
250 TEMP = TEMP - A(I,J)*X(IX)
251 IX = IX + INCX
252 120 CONTINUE
253 IF (NOUNIT) TEMP = TEMP/A(J,J)
254 ELSE
255 DO 130 I = 1,J - 1
256 TEMP = TEMP - CONJG(A(I,J))*X(IX)
257 IX = IX + INCX
258 130 CONTINUE
259 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
260 END IF
261 X(JX) = TEMP
262 JX = JX + INCX
263 140 CONTINUE
264 END IF
265 ELSE
266 IF (INCX.EQ.1) THEN
267 DO 170 J = N,1,-1
268 TEMP = X(J)
269 IF (NOCONJ) THEN
270 DO 150 I = N,J + 1,-1
271 TEMP = TEMP - A(I,J)*X(I)
272 150 CONTINUE
273 IF (NOUNIT) TEMP = TEMP/A(J,J)
274 ELSE
275 DO 160 I = N,J + 1,-1
276 TEMP = TEMP - CONJG(A(I,J))*X(I)
277 160 CONTINUE
278 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
279 END IF
280 X(J) = TEMP
281 170 CONTINUE
282 ELSE
283 KX = KX + (N-1)*INCX
284 JX = KX
285 DO 200 J = N,1,-1
286 IX = KX
287 TEMP = X(JX)
288 IF (NOCONJ) THEN
289 DO 180 I = N,J + 1,-1
290 TEMP = TEMP - A(I,J)*X(IX)
291 IX = IX - INCX
292 180 CONTINUE
293 IF (NOUNIT) TEMP = TEMP/A(J,J)
294 ELSE
295 DO 190 I = N,J + 1,-1
296 TEMP = TEMP - CONJG(A(I,J))*X(IX)
297 IX = IX - INCX
298 190 CONTINUE
299 IF (NOUNIT) TEMP = TEMP/CONJG(A(J,J))
300 END IF
301 X(JX) = TEMP
302 JX = JX - INCX
303 200 CONTINUE
304 END IF
305 END IF
306 END IF
307 *
308 RETURN
309 *
310 * End of CTRSV .
311 *
312 END