1 SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) 2* .. Scalar Arguments .. 3 DOUBLE COMPLEX ALPHA 4 INTEGER INCX,INCY,N 5 CHARACTER UPLO 6* .. 7* .. Array Arguments .. 8 DOUBLE COMPLEX AP(*),X(*),Y(*) 9* .. 10* 11* Purpose 12* ======= 13* 14* ZHPR2 performs the hermitian rank 2 operation 15* 16* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, 17* 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. 20* 21* Arguments 22* ========== 23* 24* UPLO - CHARACTER*1. 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: 28* 29* UPLO = 'U' or 'u' The upper triangular part of A is 30* supplied in AP. 31* 32* UPLO = 'L' or 'l' The lower triangular part of A is 33* supplied in AP. 34* 35* Unchanged on exit. 36* 37* N - INTEGER. 38* On entry, N specifies the order of the matrix A. 39* N must be at least zero. 40* Unchanged on exit. 41* 42* ALPHA - COMPLEX*16 . 43* On entry, ALPHA specifies the scalar alpha. 44* Unchanged on exit. 45* 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 49* element vector x. 50* Unchanged on exit. 51* 52* INCX - INTEGER. 53* On entry, INCX specifies the increment for the elements of 54* X. INCX must not be zero. 55* Unchanged on exit. 56* 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 60* element vector y. 61* Unchanged on exit. 62* 63* INCY - INTEGER. 64* On entry, INCY specifies the increment for the elements of 65* Y. INCY must not be zero. 66* Unchanged on exit. 67* 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 76* updated matrix. 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 83* updated matrix. 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 86* are set to zero. 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 TEMP1,TEMP2 107 INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY 108* .. 109* .. External Functions .. 110 LOGICAL LSAME 111 EXTERNAL LSAME 112* .. 113* .. External Subroutines .. 114 EXTERNAL XERBLA 115* .. 116* .. Intrinsic Functions .. 117 INTRINSIC DBLE,DCONJG 118* .. 119* 120* Test the input parameters. 121* 122 INFO = 0 123 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN 124 INFO = 1 125 ELSE IF (N.LT.0) THEN 126 INFO = 2 127 ELSE IF (INCX.EQ.0) THEN 128 INFO = 5 129 ELSE IF (INCY.EQ.0) THEN 130 INFO = 7 131 END IF 132 IF (INFO.NE.0) THEN 133 CALL XERBLA('ZHPR2 ',INFO) 134 RETURN 135 END IF 136* 137* Quick return if possible. 138* 139 IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN 140* 141* Set up the start points in X and Y if the increments are not both 142* unity. 143* 144 IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN 145 IF (INCX.GT.0) THEN 146 KX = 1 147 ELSE 148 KX = 1 - (N-1)*INCX 149 END IF 150 IF (INCY.GT.0) THEN 151 KY = 1 152 ELSE 153 KY = 1 - (N-1)*INCY 154 END IF 155 JX = KX 156 JY = KY 157 END IF 158* 159* Start the operations. In this version the elements of the array AP 160* are accessed sequentially with one pass through AP. 161* 162 KK = 1 163 IF (LSAME(UPLO,'U')) THEN 164* 165* Form A when upper triangle is stored in AP. 166* 167 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN 168 DO 20 J = 1,N 169 IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN 170 TEMP1 = ALPHA*DCONJG(Y(J)) 171 TEMP2 = DCONJG(ALPHA*X(J)) 172 K = KK 173 DO 10 I = 1,J - 1 174 AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 175 K = K + 1 176 10 CONTINUE 177 AP(KK+J-1) = DBLE(AP(KK+J-1)) + 178 + DBLE(X(J)*TEMP1+Y(J)*TEMP2) 179 ELSE 180 AP(KK+J-1) = DBLE(AP(KK+J-1)) 181 END IF 182 KK = KK + J 183 20 CONTINUE 184 ELSE 185 DO 40 J = 1,N 186 IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN 187 TEMP1 = ALPHA*DCONJG(Y(JY)) 188 TEMP2 = DCONJG(ALPHA*X(JX)) 189 IX = KX 190 IY = KY 191 DO 30 K = KK,KK + J - 2 192 AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 193 IX = IX + INCX 194 IY = IY + INCY 195 30 CONTINUE 196 AP(KK+J-1) = DBLE(AP(KK+J-1)) + 197 + DBLE(X(JX)*TEMP1+Y(JY)*TEMP2) 198 ELSE 199 AP(KK+J-1) = DBLE(AP(KK+J-1)) 200 END IF 201 JX = JX + INCX 202 JY = JY + INCY 203 KK = KK + J 204 40 CONTINUE 205 END IF 206 ELSE 207* 208* Form A when lower triangle is stored in AP. 209* 210 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN 211 DO 60 J = 1,N 212 IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN 213 TEMP1 = ALPHA*DCONJG(Y(J)) 214 TEMP2 = DCONJG(ALPHA*X(J)) 215 AP(KK) = DBLE(AP(KK)) + 216 + DBLE(X(J)*TEMP1+Y(J)*TEMP2) 217 K = KK + 1 218 DO 50 I = J + 1,N 219 AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 220 K = K + 1 221 50 CONTINUE 222 ELSE 223 AP(KK) = DBLE(AP(KK)) 224 END IF 225 KK = KK + N - J + 1 226 60 CONTINUE 227 ELSE 228 DO 80 J = 1,N 229 IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN 230 TEMP1 = ALPHA*DCONJG(Y(JY)) 231 TEMP2 = DCONJG(ALPHA*X(JX)) 232 AP(KK) = DBLE(AP(KK)) + 233 + DBLE(X(JX)*TEMP1+Y(JY)*TEMP2) 234 IX = JX 235 IY = JY 236 DO 70 K = KK + 1,KK + N - J 237 IX = IX + INCX 238 IY = IY + INCY 239 AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 240 70 CONTINUE 241 ELSE 242 AP(KK) = DBLE(AP(KK)) 243 END IF 244 JX = JX + INCX 245 JY = JY + INCY 246 KK = KK + N - J + 1 247 80 CONTINUE 248 END IF 249 END IF 250* 251 RETURN 252* 253* End of ZHPR2 . 254* 255 END 256