1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
|
SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
* .. Scalar Arguments ..
DOUBLE COMPLEX ALPHA
INTEGER INCX,INCY,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE COMPLEX AP(*),X(*),Y(*)
* ..
*
* Purpose
* =======
*
* ZHPR2 performs the hermitian rank 2 operation
*
* A := alpha*x*y**H + conjg( alpha )*y*x**H + A,
*
* where alpha is a scalar, x and y are n element vectors and A is an
* n by n hermitian matrix, supplied in packed form.
*
* Arguments
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the matrix A is supplied in the packed
* array AP as follows:
*
* UPLO = 'U' or 'u' The upper triangular part of A is
* supplied in AP.
*
* UPLO = 'L' or 'l' The lower triangular part of A is
* supplied in AP.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y.
* Unchanged on exit.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
* and a( 2, 2 ) respectively, and so on. On exit, the array
* AP is overwritten by the upper triangular part of the
* updated matrix.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
* and a( 3, 1 ) respectively, and so on. On exit, the array
* AP is overwritten by the lower triangular part of the
* updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
* Further Details
* ===============
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE COMPLEX ZERO
PARAMETER (ZERO= (0.0D+0,0.0D+0))
* ..
* .. Local Scalars ..
DOUBLE COMPLEX TEMP1,TEMP2
INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLI
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE,DCONJG
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
ELSE IF (INCY.EQ.0) THEN
INFO = 7
END IF
IF (INFO.NE.0) THEN
CALL XERBLI('ZHPR2 ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Set up the start points in X and Y if the increments are not both
* unity.
*
IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
JX = KX
JY = KY
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
KK = 1
IF (LSAME(UPLO,'U')) THEN
*
* Form A when upper triangle is stored in AP.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 20 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
TEMP1 = ALPHA*DCONJG(Y(J))
TEMP2 = DCONJG(ALPHA*X(J))
K = KK
DO 10 I = 1,J - 1
AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
K = K + 1
10 CONTINUE
AP(KK+J-1) = DBLE(AP(KK+J-1)) +
+ DBLE(X(J)*TEMP1+Y(J)*TEMP2)
ELSE
AP(KK+J-1) = DBLE(AP(KK+J-1))
END IF
KK = KK + J
20 CONTINUE
ELSE
DO 40 J = 1,N
IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
TEMP1 = ALPHA*DCONJG(Y(JY))
TEMP2 = DCONJG(ALPHA*X(JX))
IX = KX
IY = KY
DO 30 K = KK,KK + J - 2
AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
IX = IX + INCX
IY = IY + INCY
30 CONTINUE
AP(KK+J-1) = DBLE(AP(KK+J-1)) +
+ DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
ELSE
AP(KK+J-1) = DBLE(AP(KK+J-1))
END IF
JX = JX + INCX
JY = JY + INCY
KK = KK + J
40 CONTINUE
END IF
ELSE
*
* Form A when lower triangle is stored in AP.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
TEMP1 = ALPHA*DCONJG(Y(J))
TEMP2 = DCONJG(ALPHA*X(J))
AP(KK) = DBLE(AP(KK)) +
+ DBLE(X(J)*TEMP1+Y(J)*TEMP2)
K = KK + 1
DO 50 I = J + 1,N
AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
K = K + 1
50 CONTINUE
ELSE
AP(KK) = DBLE(AP(KK))
END IF
KK = KK + N - J + 1
60 CONTINUE
ELSE
DO 80 J = 1,N
IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
TEMP1 = ALPHA*DCONJG(Y(JY))
TEMP2 = DCONJG(ALPHA*X(JX))
AP(KK) = DBLE(AP(KK)) +
+ DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
IX = JX
IY = JY
DO 70 K = KK + 1,KK + N - J
IX = IX + INCX
IY = IY + INCY
AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
70 CONTINUE
ELSE
AP(KK) = DBLE(AP(KK))
END IF
JX = JX + INCX
JY = JY + INCY
KK = KK + N - J + 1
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHPR2 .
*
END
|