1
SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
2
* .. Scalar Arguments ..
7
* .. Array Arguments ..
8
DOUBLE COMPLEX AP(*),X(*),Y(*)
14
* ZHPR2 performs the hermitian rank 2 operation
16
* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
18
* where alpha is a scalar, x and y are n element vectors and A is an
19
* n by n hermitian matrix, supplied in packed form.
25
* On entry, UPLO specifies whether the upper or lower
26
* triangular part of the matrix A is supplied in the packed
27
* array AP as follows:
29
* UPLO = 'U' or 'u' The upper triangular part of A is
32
* UPLO = 'L' or 'l' The lower triangular part of A is
38
* On entry, N specifies the order of the matrix A.
39
* N must be at least zero.
42
* ALPHA - COMPLEX*16 .
43
* On entry, ALPHA specifies the scalar alpha.
46
* X - COMPLEX*16 array of dimension at least
47
* ( 1 + ( n - 1 )*abs( INCX ) ).
48
* Before entry, the incremented array X must contain the n
53
* On entry, INCX specifies the increment for the elements of
54
* X. INCX must not be zero.
57
* Y - COMPLEX*16 array of dimension at least
58
* ( 1 + ( n - 1 )*abs( INCY ) ).
59
* Before entry, the incremented array Y must contain the n
64
* On entry, INCY specifies the increment for the elements of
65
* Y. INCY must not be zero.
68
* AP - COMPLEX*16 array of DIMENSION at least
69
* ( ( n*( n + 1 ) )/2 ).
70
* Before entry with UPLO = 'U' or 'u', the array AP must
71
* contain the upper triangular part of the hermitian matrix
72
* packed sequentially, column by column, so that AP( 1 )
73
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
74
* and a( 2, 2 ) respectively, and so on. On exit, the array
75
* AP is overwritten by the upper triangular part of the
77
* Before entry with UPLO = 'L' or 'l', the array AP must
78
* contain the lower triangular part of the hermitian matrix
79
* packed sequentially, column by column, so that AP( 1 )
80
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
81
* and a( 3, 1 ) respectively, and so on. On exit, the array
82
* AP is overwritten by the lower triangular part of the
84
* Note that the imaginary parts of the diagonal elements need
85
* not be set, they are assumed to be zero, and on exit they
89
* Level 2 Blas routine.
91
* -- Written on 22-October-1986.
92
* Jack Dongarra, Argonne National Lab.
93
* Jeremy Du Croz, Nag Central Office.
94
* Sven Hammarling, Nag Central Office.
95
* Richard Hanson, Sandia National Labs.
100
PARAMETER (ZERO= (0.0D+0,0.0D+0))
102
* .. Local Scalars ..
103
DOUBLE COMPLEX TEMP1,TEMP2
104
INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
106
* .. External Functions ..
110
* .. External Subroutines ..
113
* .. Intrinsic Functions ..
114
INTRINSIC DBLE,DCONJG
117
* Test the input parameters.
120
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
122
ELSE IF (N.LT.0) THEN
124
ELSE IF (INCX.EQ.0) THEN
126
ELSE IF (INCY.EQ.0) THEN
130
CALL XERBLA('ZHPR2 ',INFO)
134
* Quick return if possible.
136
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
138
* Set up the start points in X and Y if the increments are not both
141
IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
156
* Start the operations. In this version the elements of the array AP
157
* are accessed sequentially with one pass through AP.
160
IF (LSAME(UPLO,'U')) THEN
162
* Form A when upper triangle is stored in AP.
164
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
166
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
167
TEMP1 = ALPHA*DCONJG(Y(J))
168
TEMP2 = DCONJG(ALPHA*X(J))
171
AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
174
AP(KK+J-1) = DBLE(AP(KK+J-1)) +
175
+ DBLE(X(J)*TEMP1+Y(J)*TEMP2)
177
AP(KK+J-1) = DBLE(AP(KK+J-1))
183
IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
184
TEMP1 = ALPHA*DCONJG(Y(JY))
185
TEMP2 = DCONJG(ALPHA*X(JX))
188
DO 30 K = KK,KK + J - 2
189
AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
193
AP(KK+J-1) = DBLE(AP(KK+J-1)) +
194
+ DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
196
AP(KK+J-1) = DBLE(AP(KK+J-1))
205
* Form A when lower triangle is stored in AP.
207
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
209
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
210
TEMP1 = ALPHA*DCONJG(Y(J))
211
TEMP2 = DCONJG(ALPHA*X(J))
212
AP(KK) = DBLE(AP(KK)) +
213
+ DBLE(X(J)*TEMP1+Y(J)*TEMP2)
216
AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
220
AP(KK) = DBLE(AP(KK))
226
IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
227
TEMP1 = ALPHA*DCONJG(Y(JY))
228
TEMP2 = DCONJG(ALPHA*X(JX))
229
AP(KK) = DBLE(AP(KK)) +
230
+ DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
233
DO 70 K = KK + 1,KK + N - J
236
AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
239
AP(KK) = DBLE(AP(KK))