share/lapack/blas/fortran/zher2.f

   1       SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
   2 *     .. Scalar Arguments ..
   3       COMPLEX*16         ALPHA
   4       INTEGER            INCX, INCY, LDA, N
   5       CHARACTER*1        UPLO
   6 *     .. Array Arguments ..
   7       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
   8 *     ..
   9 *
  10 *  Purpose
  11 *  =======
  12 *
  13 *  ZHER2  performs the hermitian rank 2 operation
  14 *
  15 *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  16 *
  17 *  where alpha is a scalar, x and y are n element vectors and A is an n
  18 *  by n hermitian matrix.
  19 *
  20 *  Parameters
  21 *  ==========
  22 *
  23 *  UPLO   - CHARACTER*1.
  24 *           On entry, UPLO specifies whether the upper or lower
  25 *           triangular part of the array A is to be referenced as
  26 *           follows:
  27 *
  28 *              UPLO = 'U' or 'u'   Only the upper triangular part of A
  29 *                                  is to be referenced.
  30 *
  31 *              UPLO = 'L' or 'l'   Only the lower triangular part of A
  32 *                                  is to be referenced.
  33 *
  34 *           Unchanged on exit.
  35 *
  36 *  N      - INTEGER.
  37 *           On entry, N specifies the order of the matrix A.
  38 *           N must be at least zero.
  39 *           Unchanged on exit.
  40 *
  41 *  ALPHA  - COMPLEX*16      .
  42 *           On entry, ALPHA specifies the scalar alpha.
  43 *           Unchanged on exit.
  44 *
  45 *  X      - COMPLEX*16       array of dimension at least
  46 *           ( 1 + ( n - 1 )*abs( INCX ) ).
  47 *           Before entry, the incremented array X must contain the n
  48 *           element vector x.
  49 *           Unchanged on exit.
  50 *
  51 *  INCX   - INTEGER.
  52 *           On entry, INCX specifies the increment for the elements of
  53 *           X. INCX must not be zero.
  54 *           Unchanged on exit.
  55 *
  56 *  Y      - COMPLEX*16       array of dimension at least
  57 *           ( 1 + ( n - 1 )*abs( INCY ) ).
  58 *           Before entry, the incremented array Y must contain the n
  59 *           element vector y.
  60 *           Unchanged on exit.
  61 *
  62 *  INCY   - INTEGER.
  63 *           On entry, INCY specifies the increment for the elements of
  64 *           Y. INCY must not be zero.
  65 *           Unchanged on exit.
  66 *
  67 *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
  68 *           Before entry with  UPLO = 'U' or 'u', the leading n by n
  69 *           upper triangular part of the array A must contain the upper
  70 *           triangular part of the hermitian matrix and the strictly
  71 *           lower triangular part of A is not referenced. On exit, the
  72 *           upper triangular part of the array A is overwritten by the
  73 *           upper triangular part of the updated matrix.
  74 *           Before entry with UPLO = 'L' or 'l', the leading n by n
  75 *           lower triangular part of the array A must contain the lower
  76 *           triangular part of the hermitian matrix and the strictly
  77 *           upper triangular part of A is not referenced. On exit, the
  78 *           lower triangular part of the array A is overwritten by the
  79 *           lower triangular part of the updated matrix.
  80 *           Note that the imaginary parts of the diagonal elements need
  81 *           not be set, they are assumed to be zero, and on exit they
  82 *           are set to zero.
  83 *
  84 *  LDA    - INTEGER.
  85 *           On entry, LDA specifies the first dimension of A as declared
  86 *           in the calling (sub) program. LDA must be at least
  87 *           max( 1, n ).
  88 *           Unchanged on exit.
  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 *     .. Parameters ..
 101       COMPLEX*16         ZERO
 102       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 103 *     .. Local Scalars ..
 104       COMPLEX*16         TEMP1, TEMP2
 105       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
 106 *     .. External Functions ..
 107       LOGICAL            LSAME
 108       EXTERNAL           LSAME
 109 *     .. External Subroutines ..
 110       EXTERNAL           XERBLA
 111 *     .. Intrinsic Functions ..
 112       INTRINSIC          DCONJG, MAX, DBLE
 113 *     ..
 114 *     .. Executable Statements ..
 115 *
 116 *     Test the input parameters.
 117 *
 118       INFO = 0
 119       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
 120      $         .NOT.LSAME( UPLO, 'L' )      )THEN
 121          INFO = 1
 122       ELSE IF( N.LT.0 )THEN
 123          INFO = 2
 124       ELSE IF( INCX.EQ.0 )THEN
 125          INFO = 5
 126       ELSE IF( INCY.EQ.0 )THEN
 127          INFO = 7
 128       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
 129          INFO = 9
 130       END IF
 131       IF( INFO.NE.0 )THEN
 132          CALL XERBLA( 'ZHER2 ', INFO )
 133          RETURN
 134       END IF
 135 *
 136 *     Quick return if possible.
 137 *
 138       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
 139      $   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 A are
 160 *     accessed sequentially with one pass through the triangular part
 161 *     of A.
 162 *
 163       IF( LSAME( UPLO, 'U' ) )THEN
 164 *
 165 *        Form  A  when A is stored in the upper triangle.
 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                   DO 10, I = 1, J - 1
 173                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
 174    10             CONTINUE
 175                   A( J, J ) = DBLE( A( J, J ) ) +
 176      $                        DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
 177                ELSE
 178                   A( J, J ) = DBLE( A( J, J ) )
 179                END IF
 180    20       CONTINUE
 181          ELSE
 182             DO 40, J = 1, N
 183                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
 184                   TEMP1 = ALPHA*DCONJG( Y( JY ) )
 185                   TEMP2 = DCONJG( ALPHA*X( JX ) )
 186                   IX    = KX
 187                   IY    = KY
 188                   DO 30, I = 1, J - 1
 189                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
 190      $                                     + Y( IY )*TEMP2
 191                      IX        = IX        + INCX
 192                      IY        = IY        + INCY
 193    30             CONTINUE
 194                   A( J, J ) = DBLE( A( J, J ) ) +
 195      $                        DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
 196                ELSE
 197                   A( J, J ) = DBLE( A( J, J ) )
 198                END IF
 199                JX = JX + INCX
 200                JY = JY + INCY
 201    40       CONTINUE
 202          END IF
 203       ELSE
 204 *
 205 *        Form  A  when A is stored in the lower triangle.
 206 *
 207          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
 208             DO 60, J = 1, N
 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                   A( J, J ) = DBLE( A( J, J ) ) +
 213      $                        DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
 214                   DO 50, I = J + 1, N
 215                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
 216    50             CONTINUE
 217                ELSE
 218                   A( J, J ) = DBLE( A( J, J ) )
 219                END IF
 220    60       CONTINUE
 221          ELSE
 222             DO 80, J = 1, N
 223                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
 224                   TEMP1     = ALPHA*DCONJG( Y( JY ) )
 225                   TEMP2     = DCONJG( ALPHA*X( JX ) )
 226                   A( J, J ) = DBLE( A( J, J ) ) +
 227      $                        DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
 228                   IX        = JX
 229                   IY        = JY
 230                   DO 70, I = J + 1, N
 231                      IX        = IX        + INCX
 232                      IY        = IY        + INCY
 233                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
 234      $                                     + Y( IY )*TEMP2
 235    70             CONTINUE
 236                ELSE
 237                   A( J, J ) = DBLE( A( J, J ) )
 238                END IF
 239                JX = JX + INCX
 240                JY = JY + INCY
 241    80       CONTINUE
 242          END IF
 243       END IF
 244 *
 245       RETURN
 246 *
 247 *     End of ZHER2 .
 248 *
 249       END