LAPACK: BLAS/SRC/zhpr2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
00002 *     .. Scalar Arguments ..
00003       DOUBLE COMPLEX ALPHA
00004       INTEGER INCX,INCY,N
00005       CHARACTER UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE COMPLEX AP(*),X(*),Y(*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  ZHPR2  performs the hermitian rank 2 operation
00015 *
00016 *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
00017 *
00018 *  where alpha is a scalar, x and y are n element vectors and A is an
00019 *  n by n hermitian matrix, supplied in packed form.
00020 *
00021 *  Arguments
00022 *  ==========
00023 *
00024 *  UPLO   - CHARACTER*1.
00025 *           On entry, UPLO specifies whether the upper or lower
00026 *           triangular part of the matrix A is supplied in the packed
00027 *           array AP as follows:
00028 *
00029 *              UPLO = 'U' or 'u'   The upper triangular part of A is
00030 *                                  supplied in AP.
00031 *
00032 *              UPLO = 'L' or 'l'   The lower triangular part of A is
00033 *                                  supplied in AP.
00034 *
00035 *           Unchanged on exit.
00036 *
00037 *  N      - INTEGER.
00038 *           On entry, N specifies the order of the matrix A.
00039 *           N must be at least zero.
00040 *           Unchanged on exit.
00041 *
00042 *  ALPHA  - COMPLEX*16      .
00043 *           On entry, ALPHA specifies the scalar alpha.
00044 *           Unchanged on exit.
00045 *
00046 *  X      - COMPLEX*16       array of dimension at least
00047 *           ( 1 + ( n - 1 )*abs( INCX ) ).
00048 *           Before entry, the incremented array X must contain the n
00049 *           element vector x.
00050 *           Unchanged on exit.
00051 *
00052 *  INCX   - INTEGER.
00053 *           On entry, INCX specifies the increment for the elements of
00054 *           X. INCX must not be zero.
00055 *           Unchanged on exit.
00056 *
00057 *  Y      - COMPLEX*16       array of dimension at least
00058 *           ( 1 + ( n - 1 )*abs( INCY ) ).
00059 *           Before entry, the incremented array Y must contain the n
00060 *           element vector y.
00061 *           Unchanged on exit.
00062 *
00063 *  INCY   - INTEGER.
00064 *           On entry, INCY specifies the increment for the elements of
00065 *           Y. INCY must not be zero.
00066 *           Unchanged on exit.
00067 *
00068 *  AP     - COMPLEX*16       array of DIMENSION at least
00069 *           ( ( n*( n + 1 ) )/2 ).
00070 *           Before entry with  UPLO = 'U' or 'u', the array AP must
00071 *           contain the upper triangular part of the hermitian matrix
00072 *           packed sequentially, column by column, so that AP( 1 )
00073 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
00074 *           and a( 2, 2 ) respectively, and so on. On exit, the array
00075 *           AP is overwritten by the upper triangular part of the
00076 *           updated matrix.
00077 *           Before entry with UPLO = 'L' or 'l', the array AP must
00078 *           contain the lower triangular part of the hermitian matrix
00079 *           packed sequentially, column by column, so that AP( 1 )
00080 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
00081 *           and a( 3, 1 ) respectively, and so on. On exit, the array
00082 *           AP is overwritten by the lower triangular part of the
00083 *           updated matrix.
00084 *           Note that the imaginary parts of the diagonal elements need
00085 *           not be set, they are assumed to be zero, and on exit they
00086 *           are set to zero.
00087 *
00088 *  Further Details
00089 *  ===============
00090 *
00091 *  Level 2 Blas routine.
00092 *
00093 *  -- Written on 22-October-1986.
00094 *     Jack Dongarra, Argonne National Lab.
00095 *     Jeremy Du Croz, Nag Central Office.
00096 *     Sven Hammarling, Nag Central Office.
00097 *     Richard Hanson, Sandia National Labs.
00098 *
00099 *  =====================================================================
00100 *
00101 *     .. Parameters ..
00102       DOUBLE COMPLEX ZERO
00103       PARAMETER (ZERO= (0.0D+0,0.0D+0))
00104 *     ..
00105 *     .. Local Scalars ..
00106       DOUBLE COMPLEX TEMP1,TEMP2
00107       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
00108 *     ..
00109 *     .. External Functions ..
00110       LOGICAL LSAME
00111       EXTERNAL LSAME
00112 *     ..
00113 *     .. External Subroutines ..
00114       EXTERNAL XERBLA
00115 *     ..
00116 *     .. Intrinsic Functions ..
00117       INTRINSIC DBLE,DCONJG
00118 *     ..
00119 *
00120 *     Test the input parameters.
00121 *
00122       INFO = 0
00123       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
00124           INFO = 1
00125       ELSE IF (N.LT.0) THEN
00126           INFO = 2
00127       ELSE IF (INCX.EQ.0) THEN
00128           INFO = 5
00129       ELSE IF (INCY.EQ.0) THEN
00130           INFO = 7
00131       END IF
00132       IF (INFO.NE.0) THEN
00133           CALL XERBLA('ZHPR2 ',INFO)
00134           RETURN
00135       END IF
00136 *
00137 *     Quick return if possible.
00138 *
00139       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
00140 *
00141 *     Set up the start points in X and Y if the increments are not both
00142 *     unity.
00143 *
00144       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
00145           IF (INCX.GT.0) THEN
00146               KX = 1
00147           ELSE
00148               KX = 1 - (N-1)*INCX
00149           END IF
00150           IF (INCY.GT.0) THEN
00151               KY = 1
00152           ELSE
00153               KY = 1 - (N-1)*INCY
00154           END IF
00155           JX = KX
00156           JY = KY
00157       END IF
00158 *
00159 *     Start the operations. In this version the elements of the array AP
00160 *     are accessed sequentially with one pass through AP.
00161 *
00162       KK = 1
00163       IF (LSAME(UPLO,'U')) THEN
00164 *
00165 *        Form  A  when upper triangle is stored in AP.
00166 *
00167           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00168               DO 20 J = 1,N
00169                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
00170                       TEMP1 = ALPHA*DCONJG(Y(J))
00171                       TEMP2 = DCONJG(ALPHA*X(J))
00172                       K = KK
00173                       DO 10 I = 1,J - 1
00174                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
00175                           K = K + 1
00176    10                 CONTINUE
00177                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
00178      +                             DBLE(X(J)*TEMP1+Y(J)*TEMP2)
00179                   ELSE
00180                       AP(KK+J-1) = DBLE(AP(KK+J-1))
00181                   END IF
00182                   KK = KK + J
00183    20         CONTINUE
00184           ELSE
00185               DO 40 J = 1,N
00186                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
00187                       TEMP1 = ALPHA*DCONJG(Y(JY))
00188                       TEMP2 = DCONJG(ALPHA*X(JX))
00189                       IX = KX
00190                       IY = KY
00191                       DO 30 K = KK,KK + J - 2
00192                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
00193                           IX = IX + INCX
00194                           IY = IY + INCY
00195    30                 CONTINUE
00196                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
00197      +                             DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
00198                   ELSE
00199                       AP(KK+J-1) = DBLE(AP(KK+J-1))
00200                   END IF
00201                   JX = JX + INCX
00202                   JY = JY + INCY
00203                   KK = KK + J
00204    40         CONTINUE
00205           END IF
00206       ELSE
00207 *
00208 *        Form  A  when lower triangle is stored in AP.
00209 *
00210           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00211               DO 60 J = 1,N
00212                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
00213                       TEMP1 = ALPHA*DCONJG(Y(J))
00214                       TEMP2 = DCONJG(ALPHA*X(J))
00215                       AP(KK) = DBLE(AP(KK)) +
00216      +                         DBLE(X(J)*TEMP1+Y(J)*TEMP2)
00217                       K = KK + 1
00218                       DO 50 I = J + 1,N
00219                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
00220                           K = K + 1
00221    50                 CONTINUE
00222                   ELSE
00223                       AP(KK) = DBLE(AP(KK))
00224                   END IF
00225                   KK = KK + N - J + 1
00226    60         CONTINUE
00227           ELSE
00228               DO 80 J = 1,N
00229                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
00230                       TEMP1 = ALPHA*DCONJG(Y(JY))
00231                       TEMP2 = DCONJG(ALPHA*X(JX))
00232                       AP(KK) = DBLE(AP(KK)) +
00233      +                         DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
00234                       IX = JX
00235                       IY = JY
00236                       DO 70 K = KK + 1,KK + N - J
00237                           IX = IX + INCX
00238                           IY = IY + INCY
00239                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
00240    70                 CONTINUE
00241                   ELSE
00242                       AP(KK) = DBLE(AP(KK))
00243                   END IF
00244                   JX = JX + INCX
00245                   JY = JY + INCY
00246                   KK = KK + N - J + 1
00247    80         CONTINUE
00248           END IF
00249       END IF
00250 *
00251       RETURN
00252 *
00253 *     End of ZHPR2 .
00254 *
00255       END