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'*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' 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 - 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' )*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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas