2
SUBROUTINE DTPMV (UPLO, TRANS, DIAG, N, AP, X, INCX)
3
C***BEGIN PROLOGUE DTPMV
4
C***PURPOSE Perform one of the matrix-vector operations.
5
C***LIBRARY SLATEC (BLAS)
7
C***TYPE DOUBLE PRECISION (STPMV-S, DTPMV-D, CTPMV-C)
8
C***KEYWORDS LEVEL 2 BLAS, LINEAR ALGEBRA
9
C***AUTHOR Dongarra, J. J., (ANL)
11
C Hammarling, S., (NAG)
12
C Hanson, R. J., (SNLA)
15
C DTPMV performs one of the matrix-vector operations
17
C x := A*x, or x := A'*x,
19
C where x is an n element vector and A is an n by n unit, or non-unit,
20
C upper or lower triangular matrix, supplied in packed form.
26
C On entry, UPLO specifies whether the matrix is an upper or
27
C lower triangular matrix as follows:
29
C UPLO = 'U' or 'u' A is an upper triangular matrix.
31
C UPLO = 'L' or 'l' A is a lower triangular matrix.
35
C TRANS - CHARACTER*1.
36
C On entry, TRANS specifies the operation to be performed as
39
C TRANS = 'N' or 'n' x := A*x.
41
C TRANS = 'T' or 't' x := A'*x.
43
C TRANS = 'C' or 'c' x := A'*x.
48
C On entry, DIAG specifies whether or not A is unit
49
C triangular as follows:
51
C DIAG = 'U' or 'u' A is assumed to be unit triangular.
53
C DIAG = 'N' or 'n' A is not assumed to be unit
59
C On entry, N specifies the order of the matrix A.
60
C N must be at least zero.
63
C AP - DOUBLE PRECISION array of DIMENSION at least
65
C Before entry with UPLO = 'U' or 'u', the array AP must
66
C contain the upper triangular matrix packed sequentially,
67
C column by column, so that AP( 1 ) contains a( 1, 1 ),
68
C AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
69
C respectively, and so on.
70
C Before entry with UPLO = 'L' or 'l', the array AP must
71
C contain the lower triangular matrix packed sequentially,
72
C column by column, so that AP( 1 ) contains a( 1, 1 ),
73
C AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
74
C respectively, and so on.
75
C Note that when DIAG = 'U' or 'u', the diagonal elements of
76
C A are not referenced, but are assumed to be unity.
79
C X - DOUBLE PRECISION array of dimension at least
80
C ( 1 + ( n - 1 )*abs( INCX ) ).
81
C Before entry, the incremented array X must contain the n
82
C element vector x. On exit, X is overwritten with the
83
C transformed vector x.
86
C On entry, INCX specifies the increment for the elements of
87
C X. INCX must not be zero.
90
C***REFERENCES Dongarra, J. J., Du Croz, J., Hammarling, S., and
91
C Hanson, R. J. An extended set of Fortran basic linear
92
C algebra subprograms. ACM TOMS, Vol. 14, No. 1,
93
C pp. 1-17, March 1988.
94
C***ROUTINES CALLED LSAME, XERBLA
95
C***REVISION HISTORY (YYMMDD)
97
C 910605 Modified to meet SLATEC prologue standards. Only comment
98
C lines were modified. (BKS)
99
C***END PROLOGUE DTPMV
100
C .. Scalar Arguments ..
102
CHARACTER*1 DIAG, TRANS, UPLO
103
C .. Array Arguments ..
104
DOUBLE PRECISION AP( * ), X( * )
106
DOUBLE PRECISION ZERO
107
PARAMETER ( ZERO = 0.0D+0 )
108
C .. Local Scalars ..
109
DOUBLE PRECISION TEMP
110
INTEGER I, INFO, IX, J, JX, K, KK, KX
112
C .. External Functions ..
115
C .. External Subroutines ..
117
C***FIRST EXECUTABLE STATEMENT DTPMV
119
C Test the input parameters.
122
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
123
$ .NOT.LSAME( UPLO , 'L' ) )THEN
125
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
126
$ .NOT.LSAME( TRANS, 'T' ).AND.
127
$ .NOT.LSAME( TRANS, 'C' ) )THEN
129
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
130
$ .NOT.LSAME( DIAG , 'N' ) )THEN
132
ELSE IF( N.LT.0 )THEN
134
ELSE IF( INCX.EQ.0 )THEN
138
CALL XERBLA( 'DTPMV ', INFO )
142
C Quick return if possible.
147
NOUNIT = LSAME( DIAG, 'N' )
149
C Set up the start point in X if the increment is not unity. This
150
C will be ( N - 1 )*INCX too small for descending loops.
153
KX = 1 - ( N - 1 )*INCX
154
ELSE IF( INCX.NE.1 )THEN
158
C Start the operations. In this version the elements of AP are
159
C accessed sequentially with one pass through AP.
161
IF( LSAME( TRANS, 'N' ) )THEN
165
IF( LSAME( UPLO, 'U' ) )THEN
169
IF( X( J ).NE.ZERO )THEN
173
X( I ) = X( I ) + TEMP*AP( K )
177
$ X( J ) = X( J )*AP( KK + J - 1 )
184
IF( X( JX ).NE.ZERO )THEN
187
DO 30, K = KK, KK + J - 2
188
X( IX ) = X( IX ) + TEMP*AP( K )
192
$ X( JX ) = X( JX )*AP( KK + J - 1 )
199
KK = ( N*( N + 1 ) )/2
202
IF( X( J ).NE.ZERO )THEN
205
DO 50, I = N, J + 1, -1
206
X( I ) = X( I ) + TEMP*AP( K )
210
$ X( J ) = X( J )*AP( KK - N + J )
212
KK = KK - ( N - J + 1 )
215
KX = KX + ( N - 1 )*INCX
218
IF( X( JX ).NE.ZERO )THEN
221
DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1
222
X( IX ) = X( IX ) + TEMP*AP( K )
226
$ X( JX ) = X( JX )*AP( KK - N + J )
229
KK = KK - ( N - J + 1 )
237
IF( LSAME( UPLO, 'U' ) )THEN
238
KK = ( N*( N + 1 ) )/2
243
$ TEMP = TEMP*AP( KK )
245
DO 90, I = J - 1, 1, -1
246
TEMP = TEMP + AP( K )*X( I )
253
JX = KX + ( N - 1 )*INCX
258
$ TEMP = TEMP*AP( KK )
259
DO 110, K = KK - 1, KK - J + 1, -1
261
TEMP = TEMP + AP( K )*X( IX )
274
$ TEMP = TEMP*AP( KK )
277
TEMP = TEMP + AP( K )*X( I )
281
KK = KK + ( N - J + 1 )
289
$ TEMP = TEMP*AP( KK )
290
DO 150, K = KK + 1, KK + N - J
292
TEMP = TEMP + AP( K )*X( IX )
296
KK = KK + ( N - J + 1 )