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**T*x, or x := A**H*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**T*x.
40: *
41: * TRANS = 'C' or 'c' x := A**H*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: * The vector and matrix arguments are not referenced when N = 0, or M = 0
93: *
94: * -- Written on 22-October-1986.
95: * Jack Dongarra, Argonne National Lab.
96: * Jeremy Du Croz, Nag Central Office.
97: * Sven Hammarling, Nag Central Office.
98: * Richard Hanson, Sandia National Labs.
99: *
100: * =====================================================================
101: *
102: * .. Parameters ..
103: DOUBLE COMPLEX ZERO
104: PARAMETER (ZERO= (0.0D+0,0.0D+0))
105: * ..
106: * .. Local Scalars ..
107: DOUBLE COMPLEX TEMP
108: INTEGER I,INFO,IX,J,JX,K,KK,KX
109: LOGICAL NOCONJ,NOUNIT
110: * ..
111: * .. External Functions ..
112: LOGICAL LSAME
113: EXTERNAL LSAME
114: * ..
115: * .. External Subroutines ..
116: EXTERNAL XERBLA
117: * ..
118: * .. Intrinsic Functions ..
119: INTRINSIC DCONJG
120: * ..
121: *
122: * Test the input parameters.
123: *
124: INFO = 0
125: IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
126: INFO = 1
127: ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
128: + .NOT.LSAME(TRANS,'C')) THEN
129: INFO = 2
130: ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
131: INFO = 3
132: ELSE IF (N.LT.0) THEN
133: INFO = 4
134: ELSE IF (INCX.EQ.0) THEN
135: INFO = 7
136: END IF
137: IF (INFO.NE.0) THEN
138: CALL XERBLA('ZTPMV ',INFO)
139: RETURN
140: END IF
141: *
142: * Quick return if possible.
143: *
144: IF (N.EQ.0) RETURN
145: *
146: NOCONJ = LSAME(TRANS,'T')
147: NOUNIT = LSAME(DIAG,'N')
148: *
149: * Set up the start point in X if the increment is not unity. This
150: * will be ( N - 1 )*INCX too small for descending loops.
151: *
152: IF (INCX.LE.0) THEN
153: KX = 1 - (N-1)*INCX
154: ELSE IF (INCX.NE.1) THEN
155: KX = 1
156: END IF
157: *
158: * Start the operations. In this version the elements of AP are
159: * accessed sequentially with one pass through AP.
160: *
161: IF (LSAME(TRANS,'N')) THEN
162: *
163: * Form x:= A*x.
164: *
165: IF (LSAME(UPLO,'U')) THEN
166: KK = 1
167: IF (INCX.EQ.1) THEN
168: DO 20 J = 1,N
169: IF (X(J).NE.ZERO) THEN
170: TEMP = X(J)
171: K = KK
172: DO 10 I = 1,J - 1
173: X(I) = X(I) + TEMP*AP(K)
174: K = K + 1
175: 10 CONTINUE
176: IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
177: END IF
178: KK = KK + J
179: 20 CONTINUE
180: ELSE
181: JX = KX
182: DO 40 J = 1,N
183: IF (X(JX).NE.ZERO) THEN
184: TEMP = X(JX)
185: IX = KX
186: DO 30 K = KK,KK + J - 2
187: X(IX) = X(IX) + TEMP*AP(K)
188: IX = IX + INCX
189: 30 CONTINUE
190: IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
191: END IF
192: JX = JX + INCX
193: KK = KK + J
194: 40 CONTINUE
195: END IF
196: ELSE
197: KK = (N* (N+1))/2
198: IF (INCX.EQ.1) THEN
199: DO 60 J = N,1,-1
200: IF (X(J).NE.ZERO) THEN
201: TEMP = X(J)
202: K = KK
203: DO 50 I = N,J + 1,-1
204: X(I) = X(I) + TEMP*AP(K)
205: K = K - 1
206: 50 CONTINUE
207: IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
208: END IF
209: KK = KK - (N-J+1)
210: 60 CONTINUE
211: ELSE
212: KX = KX + (N-1)*INCX
213: JX = KX
214: DO 80 J = N,1,-1
215: IF (X(JX).NE.ZERO) THEN
216: TEMP = X(JX)
217: IX = KX
218: DO 70 K = KK,KK - (N- (J+1)),-1
219: X(IX) = X(IX) + TEMP*AP(K)
220: IX = IX - INCX
221: 70 CONTINUE
222: IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
223: END IF
224: JX = JX - INCX
225: KK = KK - (N-J+1)
226: 80 CONTINUE
227: END IF
228: END IF
229: ELSE
230: *
231: * Form x := A**T*x or x := A**H*x.
232: *
233: IF (LSAME(UPLO,'U')) THEN
234: KK = (N* (N+1))/2
235: IF (INCX.EQ.1) THEN
236: DO 110 J = N,1,-1
237: TEMP = X(J)
238: K = KK - 1
239: IF (NOCONJ) THEN
240: IF (NOUNIT) TEMP = TEMP*AP(KK)
241: DO 90 I = J - 1,1,-1
242: TEMP = TEMP + AP(K)*X(I)
243: K = K - 1
244: 90 CONTINUE
245: ELSE
246: IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
247: DO 100 I = J - 1,1,-1
248: TEMP = TEMP + DCONJG(AP(K))*X(I)
249: K = K - 1
250: 100 CONTINUE
251: END IF
252: X(J) = TEMP
253: KK = KK - J
254: 110 CONTINUE
255: ELSE
256: JX = KX + (N-1)*INCX
257: DO 140 J = N,1,-1
258: TEMP = X(JX)
259: IX = JX
260: IF (NOCONJ) THEN
261: IF (NOUNIT) TEMP = TEMP*AP(KK)
262: DO 120 K = KK - 1,KK - J + 1,-1
263: IX = IX - INCX
264: TEMP = TEMP + AP(K)*X(IX)
265: 120 CONTINUE
266: ELSE
267: IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
268: DO 130 K = KK - 1,KK - J + 1,-1
269: IX = IX - INCX
270: TEMP = TEMP + DCONJG(AP(K))*X(IX)
271: 130 CONTINUE
272: END IF
273: X(JX) = TEMP
274: JX = JX - INCX
275: KK = KK - J
276: 140 CONTINUE
277: END IF
278: ELSE
279: KK = 1
280: IF (INCX.EQ.1) THEN
281: DO 170 J = 1,N
282: TEMP = X(J)
283: K = KK + 1
284: IF (NOCONJ) THEN
285: IF (NOUNIT) TEMP = TEMP*AP(KK)
286: DO 150 I = J + 1,N
287: TEMP = TEMP + AP(K)*X(I)
288: K = K + 1
289: 150 CONTINUE
290: ELSE
291: IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
292: DO 160 I = J + 1,N
293: TEMP = TEMP + DCONJG(AP(K))*X(I)
294: K = K + 1
295: 160 CONTINUE
296: END IF
297: X(J) = TEMP
298: KK = KK + (N-J+1)
299: 170 CONTINUE
300: ELSE
301: JX = KX
302: DO 200 J = 1,N
303: TEMP = X(JX)
304: IX = JX
305: IF (NOCONJ) THEN
306: IF (NOUNIT) TEMP = TEMP*AP(KK)
307: DO 180 K = KK + 1,KK + N - J
308: IX = IX + INCX
309: TEMP = TEMP + AP(K)*X(IX)
310: 180 CONTINUE
311: ELSE
312: IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
313: DO 190 K = KK + 1,KK + N - J
314: IX = IX + INCX
315: TEMP = TEMP + DCONJG(AP(K))*X(IX)
316: 190 CONTINUE
317: END IF
318: X(JX) = TEMP
319: JX = JX + INCX
320: KK = KK + (N-J+1)
321: 200 CONTINUE
322: END IF
323: END IF
324: END IF
325: *
326: RETURN
327: *
328: * End of ZTPMV .
329: *
330: END
CVSweb interface <joel.bertrand@systella.fr>