/[ascend]/trunk/blas/dtrsv.f
ViewVC logotype

Contents of /trunk/blas/dtrsv.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (14 years ago) by aw0a
File size: 7983 byte(s)
Setting up web subdirectory in repository
1 SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
2 *
3 ************************************************************************
4 *
5 * .. Scalar Arguments ..
6 INTEGER INCX, LDA, N
7 CHARACTER*1 DIAG, TRANS, UPLO
8 * .. Array Arguments ..
9 DOUBLE PRECISION A( LDA, N ), X( * )
10 * ..
11 *
12 * Purpose
13 * =======
14 *
15 * DTRSV solves one of the systems of equations
16 *
17 * A*x = b, or A'*x = b,
18 *
19 * where b and x are n element vectors and A is an n by n unit, or
20 * non-unit, upper or lower triangular matrix.
21 *
22 * No test for singularity or near-singularity is included in this
23 * routine. Such tests must be performed before calling this routine.
24 *
25 * Parameters
26 * ==========
27 *
28 * UPLO - CHARACTER*1.
29 * On entry, UPLO specifies whether the matrix is an upper or
30 * lower triangular matrix as follows:
31 *
32 * UPLO = 'U' or 'u' A is an upper triangular matrix.
33 *
34 * UPLO = 'L' or 'l' A is a lower triangular matrix.
35 *
36 * Unchanged on exit.
37 *
38 * TRANS - CHARACTER*1.
39 * On entry, TRANS specifies the equations to be solved as
40 * follows:
41 *
42 * TRANS = 'N' or 'n' A*x = b.
43 *
44 * TRANS = 'T' or 't' A'*x = b.
45 *
46 * TRANS = 'C' or 'c' A'*x = b.
47 *
48 * Unchanged on exit.
49 *
50 * DIAG - CHARACTER*1.
51 * On entry, DIAG specifies whether or not A is unit
52 * triangular as follows:
53 *
54 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
55 *
56 * DIAG = 'N' or 'n' A is not assumed to be unit
57 * triangular.
58 *
59 * Unchanged on exit.
60 *
61 * N - INTEGER.
62 * On entry, N specifies the order of the matrix A.
63 * N must be at least zero.
64 * Unchanged on exit.
65 *
66 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
67 * Before entry with UPLO = 'U' or 'u', the leading n by n
68 * upper triangular part of the array A must contain the upper
69 * triangular matrix and the strictly lower triangular part of
70 * A is not referenced.
71 * Before entry with UPLO = 'L' or 'l', the leading n by n
72 * lower triangular part of the array A must contain the lower
73 * triangular matrix and the strictly upper triangular part of
74 * A is not referenced.
75 * Note that when DIAG = 'U' or 'u', the diagonal elements of
76 * A are not referenced either, but are assumed to be unity.
77 * Unchanged on exit.
78 *
79 * LDA - INTEGER.
80 * On entry, LDA specifies the first dimension of A as declared
81 * in the calling (sub) program. LDA must be at least
82 * max( 1, n ).
83 * Unchanged on exit.
84 *
85 * X - DOUBLE PRECISION array of dimension at least
86 * ( 1 + ( n - 1 )*abs( INCX ) ).
87 * Before entry, the incremented array X must contain the n
88 * element right-hand side vector b. On exit, X is overwritten
89 * with the solution vector x.
90 *
91 * INCX - INTEGER.
92 * On entry, INCX specifies the increment for the elements of
93 * X. INCX must not be zero.
94 * Unchanged on exit.
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 * .. Parameters ..
107 DOUBLE PRECISION ZERO
108 PARAMETER ( ZERO = 0.0D+0 )
109 * .. Local Scalars ..
110 DOUBLE PRECISION TEMP
111 INTEGER I, INFO, IX, J, JX, KX
112 LOGICAL NOUNIT
113 * .. External Functions ..
114 LOGICAL LSAME
115 EXTERNAL LSAME
116 * .. External Subroutines ..
117 EXTERNAL XERBLA
118 * .. Intrinsic Functions ..
119 INTRINSIC MAX
120 * ..
121 * .. Executable Statements ..
122 *
123 * Test the input parameters.
124 *
125 INFO = 0
126 IF ( .NOT.LSAME( UPLO , 'U' ).AND.
127 $ .NOT.LSAME( UPLO , 'L' ) )THEN
128 INFO = 1
129 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
130 $ .NOT.LSAME( TRANS, 'T' ).AND.
131 $ .NOT.LSAME( TRANS, 'C' ) )THEN
132 INFO = 2
133 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
134 $ .NOT.LSAME( DIAG , 'N' ) )THEN
135 INFO = 3
136 ELSE IF( N.LT.0 )THEN
137 INFO = 4
138 ELSE IF( LDA.LT.MAX( 1, N ) )THEN
139 INFO = 6
140 ELSE IF( INCX.EQ.0 )THEN
141 INFO = 8
142 END IF
143 IF( INFO.NE.0 )THEN
144 CALL XERBLA( 'DTRSV ', INFO )
145 RETURN
146 END IF
147 *
148 * Quick return if possible.
149 *
150 IF( N.EQ.0 )
151 $ 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 )
176 $ 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 )
188 $ X( JX ) = X( JX )/A( J, J )
189 TEMP = X( JX )
190 IX = JX
191 DO 30, I = J - 1, 1, -1
192 IX = IX - INCX
193 X( IX ) = X( IX ) - TEMP*A( I, J )
194 30 CONTINUE
195 END IF
196 JX = JX - INCX
197 40 CONTINUE
198 END IF
199 ELSE
200 IF( INCX.EQ.1 )THEN
201 DO 60, J = 1, N
202 IF( X( J ).NE.ZERO )THEN
203 IF( NOUNIT )
204 $ X( J ) = X( J )/A( J, J )
205 TEMP = X( J )
206 DO 50, I = J + 1, N
207 X( I ) = X( I ) - TEMP*A( I, J )
208 50 CONTINUE
209 END IF
210 60 CONTINUE
211 ELSE
212 JX = KX
213 DO 80, J = 1, N
214 IF( X( JX ).NE.ZERO )THEN
215 IF( NOUNIT )
216 $ X( JX ) = X( JX )/A( J, J )
217 TEMP = X( JX )
218 IX = JX
219 DO 70, I = J + 1, N
220 IX = IX + INCX
221 X( IX ) = X( IX ) - TEMP*A( I, J )
222 70 CONTINUE
223 END IF
224 JX = JX + INCX
225 80 CONTINUE
226 END IF
227 END IF
228 ELSE
229 *
230 * Form x := inv( A' )*x.
231 *
232 IF( LSAME( UPLO, 'U' ) )THEN
233 IF( INCX.EQ.1 )THEN
234 DO 100, J = 1, N
235 TEMP = X( J )
236 DO 90, I = 1, J - 1
237 TEMP = TEMP - A( I, J )*X( I )
238 90 CONTINUE
239 IF( NOUNIT )
240 $ TEMP = TEMP/A( J, J )
241 X( J ) = TEMP
242 100 CONTINUE
243 ELSE
244 JX = KX
245 DO 120, J = 1, N
246 TEMP = X( JX )
247 IX = KX
248 DO 110, I = 1, J - 1
249 TEMP = TEMP - A( I, J )*X( IX )
250 IX = IX + INCX
251 110 CONTINUE
252 IF( NOUNIT )
253 $ TEMP = TEMP/A( J, J )
254 X( JX ) = TEMP
255 JX = JX + INCX
256 120 CONTINUE
257 END IF
258 ELSE
259 IF( INCX.EQ.1 )THEN
260 DO 140, J = N, 1, -1
261 TEMP = X( J )
262 DO 130, I = N, J + 1, -1
263 TEMP = TEMP - A( I, J )*X( I )
264 130 CONTINUE
265 IF( NOUNIT )
266 $ TEMP = TEMP/A( J, J )
267 X( J ) = TEMP
268 140 CONTINUE
269 ELSE
270 KX = KX + ( N - 1 )*INCX
271 JX = KX
272 DO 160, J = N, 1, -1
273 TEMP = X( JX )
274 IX = KX
275 DO 150, I = N, J + 1, -1
276 TEMP = TEMP - A( I, J )*X( IX )
277 IX = IX - INCX
278 150 CONTINUE
279 IF( NOUNIT )
280 $ TEMP = TEMP/A( J, J )
281 X( JX ) = TEMP
282 JX = JX - INCX
283 160 CONTINUE
284 END IF
285 END IF
286 END IF
287 *
288 RETURN
289 *
290 * End of DTRSV .
291 *
292 END

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22