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