1: *> \brief \b DTPSV
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: * Definition:
9: * ===========
10: *
11: * SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
12: *
13: * .. Scalar Arguments ..
14: * INTEGER INCX,N
15: * CHARACTER DIAG,TRANS,UPLO
16: * ..
17: * .. Array Arguments ..
18: * DOUBLE PRECISION AP(*),X(*)
19: * ..
20: *
21: *
22: *> \par Purpose:
23: * =============
24: *>
25: *> \verbatim
26: *>
27: *> DTPSV solves one of the systems of equations
28: *>
29: *> A*x = b, or A**T*x = b,
30: *>
31: *> where b and x are n element vectors and A is an n by n unit, or
32: *> non-unit, upper or lower triangular matrix, supplied in packed form.
33: *>
34: *> No test for singularity or near-singularity is included in this
35: *> routine. Such tests must be performed before calling this routine.
36: *> \endverbatim
37: *
38: * Arguments:
39: * ==========
40: *
41: *> \param[in] UPLO
42: *> \verbatim
43: *> UPLO is CHARACTER*1
44: *> On entry, UPLO specifies whether the matrix is an upper or
45: *> lower triangular matrix as follows:
46: *>
47: *> UPLO = 'U' or 'u' A is an upper triangular matrix.
48: *>
49: *> UPLO = 'L' or 'l' A is a lower triangular matrix.
50: *> \endverbatim
51: *>
52: *> \param[in] TRANS
53: *> \verbatim
54: *> TRANS is CHARACTER*1
55: *> On entry, TRANS specifies the equations to be solved as
56: *> follows:
57: *>
58: *> TRANS = 'N' or 'n' A*x = b.
59: *>
60: *> TRANS = 'T' or 't' A**T*x = b.
61: *>
62: *> TRANS = 'C' or 'c' A**T*x = b.
63: *> \endverbatim
64: *>
65: *> \param[in] DIAG
66: *> \verbatim
67: *> DIAG is CHARACTER*1
68: *> On entry, DIAG specifies whether or not A is unit
69: *> triangular as follows:
70: *>
71: *> DIAG = 'U' or 'u' A is assumed to be unit triangular.
72: *>
73: *> DIAG = 'N' or 'n' A is not assumed to be unit
74: *> triangular.
75: *> \endverbatim
76: *>
77: *> \param[in] N
78: *> \verbatim
79: *> N is INTEGER
80: *> On entry, N specifies the order of the matrix A.
81: *> N must be at least zero.
82: *> \endverbatim
83: *>
84: *> \param[in] AP
85: *> \verbatim
86: *> AP is DOUBLE PRECISION array, dimension at least
87: *> ( ( n*( n + 1 ) )/2 ).
88: *> Before entry with UPLO = 'U' or 'u', the array AP must
89: *> contain the upper triangular matrix packed sequentially,
90: *> column by column, so that AP( 1 ) contains a( 1, 1 ),
91: *> AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
92: *> respectively, and so on.
93: *> Before entry with UPLO = 'L' or 'l', the array AP must
94: *> contain the lower triangular matrix packed sequentially,
95: *> column by column, so that AP( 1 ) contains a( 1, 1 ),
96: *> AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
97: *> respectively, and so on.
98: *> Note that when DIAG = 'U' or 'u', the diagonal elements of
99: *> A are not referenced, but are assumed to be unity.
100: *> \endverbatim
101: *>
102: *> \param[in,out] X
103: *> \verbatim
104: *> X is DOUBLE PRECISION array, dimension at least
105: *> ( 1 + ( n - 1 )*abs( INCX ) ).
106: *> Before entry, the incremented array X must contain the n
107: *> element right-hand side vector b. On exit, X is overwritten
108: *> with the solution vector x.
109: *> \endverbatim
110: *>
111: *> \param[in] INCX
112: *> \verbatim
113: *> INCX is INTEGER
114: *> On entry, INCX specifies the increment for the elements of
115: *> X. INCX must not be zero.
116: *> \endverbatim
117: *
118: * Authors:
119: * ========
120: *
121: *> \author Univ. of Tennessee
122: *> \author Univ. of California Berkeley
123: *> \author Univ. of Colorado Denver
124: *> \author NAG Ltd.
125: *
126: *> \ingroup double_blas_level2
127: *
128: *> \par Further Details:
129: * =====================
130: *>
131: *> \verbatim
132: *>
133: *> Level 2 Blas routine.
134: *>
135: *> -- Written on 22-October-1986.
136: *> Jack Dongarra, Argonne National Lab.
137: *> Jeremy Du Croz, Nag Central Office.
138: *> Sven Hammarling, Nag Central Office.
139: *> Richard Hanson, Sandia National Labs.
140: *> \endverbatim
141: *>
142: * =====================================================================
143: SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
144: *
145: * -- Reference BLAS level2 routine --
146: * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
147: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148: *
149: * .. Scalar Arguments ..
150: INTEGER INCX,N
151: CHARACTER DIAG,TRANS,UPLO
152: * ..
153: * .. Array Arguments ..
154: DOUBLE PRECISION AP(*),X(*)
155: * ..
156: *
157: * =====================================================================
158: *
159: * .. Parameters ..
160: DOUBLE PRECISION ZERO
161: PARAMETER (ZERO=0.0D+0)
162: * ..
163: * .. Local Scalars ..
164: DOUBLE PRECISION TEMP
165: INTEGER I,INFO,IX,J,JX,K,KK,KX
166: LOGICAL NOUNIT
167: * ..
168: * .. External Functions ..
169: LOGICAL LSAME
170: EXTERNAL LSAME
171: * ..
172: * .. External Subroutines ..
173: EXTERNAL XERBLA
174: * ..
175: *
176: * Test the input parameters.
177: *
178: INFO = 0
179: IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
180: INFO = 1
181: ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
182: + .NOT.LSAME(TRANS,'C')) THEN
183: INFO = 2
184: ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
185: INFO = 3
186: ELSE IF (N.LT.0) THEN
187: INFO = 4
188: ELSE IF (INCX.EQ.0) THEN
189: INFO = 7
190: END IF
191: IF (INFO.NE.0) THEN
192: CALL XERBLA('DTPSV ',INFO)
193: RETURN
194: END IF
195: *
196: * Quick return if possible.
197: *
198: IF (N.EQ.0) RETURN
199: *
200: NOUNIT = LSAME(DIAG,'N')
201: *
202: * Set up the start point in X if the increment is not unity. This
203: * will be ( N - 1 )*INCX too small for descending loops.
204: *
205: IF (INCX.LE.0) THEN
206: KX = 1 - (N-1)*INCX
207: ELSE IF (INCX.NE.1) THEN
208: KX = 1
209: END IF
210: *
211: * Start the operations. In this version the elements of AP are
212: * accessed sequentially with one pass through AP.
213: *
214: IF (LSAME(TRANS,'N')) THEN
215: *
216: * Form x := inv( A )*x.
217: *
218: IF (LSAME(UPLO,'U')) THEN
219: KK = (N* (N+1))/2
220: IF (INCX.EQ.1) THEN
221: DO 20 J = N,1,-1
222: IF (X(J).NE.ZERO) THEN
223: IF (NOUNIT) X(J) = X(J)/AP(KK)
224: TEMP = X(J)
225: K = KK - 1
226: DO 10 I = J - 1,1,-1
227: X(I) = X(I) - TEMP*AP(K)
228: K = K - 1
229: 10 CONTINUE
230: END IF
231: KK = KK - J
232: 20 CONTINUE
233: ELSE
234: JX = KX + (N-1)*INCX
235: DO 40 J = N,1,-1
236: IF (X(JX).NE.ZERO) THEN
237: IF (NOUNIT) X(JX) = X(JX)/AP(KK)
238: TEMP = X(JX)
239: IX = JX
240: DO 30 K = KK - 1,KK - J + 1,-1
241: IX = IX - INCX
242: X(IX) = X(IX) - TEMP*AP(K)
243: 30 CONTINUE
244: END IF
245: JX = JX - INCX
246: KK = KK - J
247: 40 CONTINUE
248: END IF
249: ELSE
250: KK = 1
251: IF (INCX.EQ.1) THEN
252: DO 60 J = 1,N
253: IF (X(J).NE.ZERO) THEN
254: IF (NOUNIT) X(J) = X(J)/AP(KK)
255: TEMP = X(J)
256: K = KK + 1
257: DO 50 I = J + 1,N
258: X(I) = X(I) - TEMP*AP(K)
259: K = K + 1
260: 50 CONTINUE
261: END IF
262: KK = KK + (N-J+1)
263: 60 CONTINUE
264: ELSE
265: JX = KX
266: DO 80 J = 1,N
267: IF (X(JX).NE.ZERO) THEN
268: IF (NOUNIT) X(JX) = X(JX)/AP(KK)
269: TEMP = X(JX)
270: IX = JX
271: DO 70 K = KK + 1,KK + N - J
272: IX = IX + INCX
273: X(IX) = X(IX) - TEMP*AP(K)
274: 70 CONTINUE
275: END IF
276: JX = JX + INCX
277: KK = KK + (N-J+1)
278: 80 CONTINUE
279: END IF
280: END IF
281: ELSE
282: *
283: * Form x := inv( A**T )*x.
284: *
285: IF (LSAME(UPLO,'U')) THEN
286: KK = 1
287: IF (INCX.EQ.1) THEN
288: DO 100 J = 1,N
289: TEMP = X(J)
290: K = KK
291: DO 90 I = 1,J - 1
292: TEMP = TEMP - AP(K)*X(I)
293: K = K + 1
294: 90 CONTINUE
295: IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
296: X(J) = TEMP
297: KK = KK + J
298: 100 CONTINUE
299: ELSE
300: JX = KX
301: DO 120 J = 1,N
302: TEMP = X(JX)
303: IX = KX
304: DO 110 K = KK,KK + J - 2
305: TEMP = TEMP - AP(K)*X(IX)
306: IX = IX + INCX
307: 110 CONTINUE
308: IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
309: X(JX) = TEMP
310: JX = JX + INCX
311: KK = KK + J
312: 120 CONTINUE
313: END IF
314: ELSE
315: KK = (N* (N+1))/2
316: IF (INCX.EQ.1) THEN
317: DO 140 J = N,1,-1
318: TEMP = X(J)
319: K = KK
320: DO 130 I = N,J + 1,-1
321: TEMP = TEMP - AP(K)*X(I)
322: K = K - 1
323: 130 CONTINUE
324: IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
325: X(J) = TEMP
326: KK = KK - (N-J+1)
327: 140 CONTINUE
328: ELSE
329: KX = KX + (N-1)*INCX
330: JX = KX
331: DO 160 J = N,1,-1
332: TEMP = X(JX)
333: IX = KX
334: DO 150 K = KK,KK - (N- (J+1)),-1
335: TEMP = TEMP - AP(K)*X(IX)
336: IX = IX - INCX
337: 150 CONTINUE
338: IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
339: X(JX) = TEMP
340: JX = JX - INCX
341: KK = KK - (N-J+1)
342: 160 CONTINUE
343: END IF
344: END IF
345: END IF
346: *
347: RETURN
348: *
349: * End of DTPSV
350: *
351: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtpsv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas