share/lapack/blas/fortran/zherk.f

   1       SUBROUTINE ZHERK( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC )
   2 *     .. Scalar Arguments ..
   3       CHARACTER          TRANS, UPLO
   4       INTEGER            K, LDA, LDC, N
   5       DOUBLE PRECISION   ALPHA, BETA
   6 *     ..
   7 *     .. Array Arguments ..
   8       COMPLEX*16         A( LDA, * ), C( LDC, * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  ZHERK  performs one of the hermitian rank k operations
  15 *
  16 *     C := alpha*A*conjg( A' ) + beta*C,
  17 *
  18 *  or
  19 *
  20 *     C := alpha*conjg( A' )*A + beta*C,
  21 *
  22 *  where  alpha and beta  are  real scalars,  C is an  n by n  hermitian
  23 *  matrix and  A  is an  n by k  matrix in the  first case and a  k by n
  24 *  matrix in the second case.
  25 *
  26 *  Parameters
  27 *  ==========
  28 *
  29 *  UPLO   - CHARACTER*1.
  30 *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  31 *           triangular  part  of the  array  C  is to be  referenced  as
  32 *           follows:
  33 *
  34 *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  35 *                                  is to be referenced.
  36 *
  37 *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  38 *                                  is to be referenced.
  39 *
  40 *           Unchanged on exit.
  41 *
  42 *  TRANS  - CHARACTER*1.
  43 *           On entry,  TRANS  specifies the operation to be performed as
  44 *           follows:
  45 *
  46 *              TRANS = 'N' or 'n'   C := alpha*A*conjg( A' ) + beta*C.
  47 *
  48 *              TRANS = 'C' or 'c'   C := alpha*conjg( A' )*A + beta*C.
  49 *
  50 *           Unchanged on exit.
  51 *
  52 *  N      - INTEGER.
  53 *           On entry,  N specifies the order of the matrix C.  N must be
  54 *           at least zero.
  55 *           Unchanged on exit.
  56 *
  57 *  K      - INTEGER.
  58 *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  59 *           of  columns   of  the   matrix   A,   and  on   entry   with
  60 *           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
  61 *           matrix A.  K must be at least zero.
  62 *           Unchanged on exit.
  63 *
  64 *  ALPHA  - DOUBLE PRECISION            .
  65 *           On entry, ALPHA specifies the scalar alpha.
  66 *           Unchanged on exit.
  67 *
  68 *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  69 *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  70 *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  71 *           part of the array  A  must contain the matrix  A,  otherwise
  72 *           the leading  k by n  part of the array  A  must contain  the
  73 *           matrix A.
  74 *           Unchanged on exit.
  75 *
  76 *  LDA    - INTEGER.
  77 *           On entry, LDA specifies the first dimension of A as declared
  78 *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  79 *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  80 *           be at least  max( 1, k ).
  81 *           Unchanged on exit.
  82 *
  83 *  BETA   - DOUBLE PRECISION.
  84 *           On entry, BETA specifies the scalar beta.
  85 *           Unchanged on exit.
  86 *
  87 *  C      - COMPLEX*16          array of DIMENSION ( LDC, n ).
  88 *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  89 *           upper triangular part of the array C must contain the upper
  90 *           triangular part  of the  hermitian matrix  and the strictly
  91 *           lower triangular part of C is not referenced.  On exit, the
  92 *           upper triangular part of the array  C is overwritten by the
  93 *           upper triangular part of the updated matrix.
  94 *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  95 *           lower triangular part of the array C must contain the lower
  96 *           triangular part  of the  hermitian matrix  and the strictly
  97 *           upper triangular part of C is not referenced.  On exit, the
  98 *           lower triangular part of the array  C is overwritten by the
  99 *           lower triangular part of the updated matrix.
 100 *           Note that the imaginary parts of the diagonal elements need
 101 *           not be set,  they are assumed to be zero,  and on exit they
 102 *           are set to zero.
 103 *
 104 *  LDC    - INTEGER.
 105 *           On entry, LDC specifies the first dimension of C as declared
 106 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
 107 *           max( 1, n ).
 108 *           Unchanged on exit.
 109 *
 110 *
 111 *  Level 3 Blas routine.
 112 *
 113 *  -- Written on 8-February-1989.
 114 *     Jack Dongarra, Argonne National Laboratory.
 115 *     Iain Duff, AERE Harwell.
 116 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 117 *     Sven Hammarling, Numerical Algorithms Group Ltd.
 118 *
 119 *  -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
 120 *     Ed Anderson, Cray Research Inc.
 121 *
 122 *
 123 *     .. External Functions ..
 124       LOGICAL            LSAME
 125       EXTERNAL           LSAME
 126 *     ..
 127 *     .. External Subroutines ..
 128       EXTERNAL           XERBLA
 129 *     ..
 130 *     .. Intrinsic Functions ..
 131       INTRINSIC          DBLE, DCMPLX, DCONJG, MAX
 132 *     ..
 133 *     .. Local Scalars ..
 134       LOGICAL            UPPER
 135       INTEGER            I, INFO, J, L, NROWA
 136       DOUBLE PRECISION   RTEMP
 137       COMPLEX*16         TEMP
 138 *     ..
 139 *     .. Parameters ..
 140       DOUBLE PRECISION   ONE, ZERO
 141       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 142 *     ..
 143 *     .. Executable Statements ..
 144 *
 145 *     Test the input parameters.
 146 *
 147       IF( LSAME( TRANS, 'N' ) ) THEN
 148          NROWA = N
 149       ELSE
 150          NROWA = K
 151       END IF
 152       UPPER = LSAME( UPLO, 'U' )
 153 *
 154       INFO = 0
 155       IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
 156          INFO = 1
 157       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND.
 158      $         ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN
 159          INFO = 2
 160       ELSE IF( N.LT.0 ) THEN
 161          INFO = 3
 162       ELSE IF( K.LT.0 ) THEN
 163          INFO = 4
 164       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
 165          INFO = 7
 166       ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
 167          INFO = 10
 168       END IF
 169       IF( INFO.NE.0 ) THEN
 170          CALL XERBLA( 'ZHERK ', INFO )
 171          RETURN
 172       END IF
 173 *
 174 *     Quick return if possible.
 175 *
 176       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
 177      $    ( BETA.EQ.ONE ) ) )RETURN
 178 *
 179 *     And when  alpha.eq.zero.
 180 *
 181       IF( ALPHA.EQ.ZERO ) THEN
 182          IF( UPPER ) THEN
 183             IF( BETA.EQ.ZERO ) THEN
 184                DO 20 J = 1, N
 185                   DO 10 I = 1, J
 186                      C( I, J ) = ZERO
 187    10             CONTINUE
 188    20          CONTINUE
 189             ELSE
 190                DO 40 J = 1, N
 191                   DO 30 I = 1, J - 1
 192                      C( I, J ) = BETA*C( I, J )
 193    30             CONTINUE
 194                   C( J, J ) = BETA*DBLE( C( J, J ) )
 195    40          CONTINUE
 196             END IF
 197          ELSE
 198             IF( BETA.EQ.ZERO ) THEN
 199                DO 60 J = 1, N
 200                   DO 50 I = J, N
 201                      C( I, J ) = ZERO
 202    50             CONTINUE
 203    60          CONTINUE
 204             ELSE
 205                DO 80 J = 1, N
 206                   C( J, J ) = BETA*DBLE( C( J, J ) )
 207                   DO 70 I = J + 1, N
 208                      C( I, J ) = BETA*C( I, J )
 209    70             CONTINUE
 210    80          CONTINUE
 211             END IF
 212          END IF
 213          RETURN
 214       END IF
 215 *
 216 *     Start the operations.
 217 *
 218       IF( LSAME( TRANS, 'N' ) ) THEN
 219 *
 220 *        Form  C := alpha*A*conjg( A' ) + beta*C.
 221 *
 222          IF( UPPER ) THEN
 223             DO 130 J = 1, N
 224                IF( BETA.EQ.ZERO ) THEN
 225                   DO 90 I = 1, J
 226                      C( I, J ) = ZERO
 227    90             CONTINUE
 228                ELSE IF( BETA.NE.ONE ) THEN
 229                   DO 100 I = 1, J - 1
 230                      C( I, J ) = BETA*C( I, J )
 231   100             CONTINUE
 232                   C( J, J ) = BETA*DBLE( C( J, J ) )
 233                ELSE
 234                   C( J, J ) = DBLE( C( J, J ) )
 235                END IF
 236                DO 120 L = 1, K
 237                   IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN
 238                      TEMP = ALPHA*DCONJG( A( J, L ) )
 239                      DO 110 I = 1, J - 1
 240                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
 241   110                CONTINUE
 242                      C( J, J ) = DBLE( C( J, J ) ) +
 243      $                           DBLE( TEMP*A( I, L ) )
 244                   END IF
 245   120          CONTINUE
 246   130       CONTINUE
 247          ELSE
 248             DO 180 J = 1, N
 249                IF( BETA.EQ.ZERO ) THEN
 250                   DO 140 I = J, N
 251                      C( I, J ) = ZERO
 252   140             CONTINUE
 253                ELSE IF( BETA.NE.ONE ) THEN
 254                   C( J, J ) = BETA*DBLE( C( J, J ) )
 255                   DO 150 I = J + 1, N
 256                      C( I, J ) = BETA*C( I, J )
 257   150             CONTINUE
 258                ELSE
 259                   C( J, J ) = DBLE( C( J, J ) )
 260                END IF
 261                DO 170 L = 1, K
 262                   IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN
 263                      TEMP = ALPHA*DCONJG( A( J, L ) )
 264                      C( J, J ) = DBLE( C( J, J ) ) +
 265      $                           DBLE( TEMP*A( J, L ) )
 266                      DO 160 I = J + 1, N
 267                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
 268   160                CONTINUE
 269                   END IF
 270   170          CONTINUE
 271   180       CONTINUE
 272          END IF
 273       ELSE
 274 *
 275 *        Form  C := alpha*conjg( A' )*A + beta*C.
 276 *
 277          IF( UPPER ) THEN
 278             DO 220 J = 1, N
 279                DO 200 I = 1, J - 1
 280                   TEMP = ZERO
 281                   DO 190 L = 1, K
 282                      TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )
 283   190             CONTINUE
 284                   IF( BETA.EQ.ZERO ) THEN
 285                      C( I, J ) = ALPHA*TEMP
 286                   ELSE
 287                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
 288                   END IF
 289   200          CONTINUE
 290                RTEMP = ZERO
 291                DO 210 L = 1, K
 292                   RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )
 293   210          CONTINUE
 294                IF( BETA.EQ.ZERO ) THEN
 295                   C( J, J ) = ALPHA*RTEMP
 296                ELSE
 297                   C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )
 298                END IF
 299   220       CONTINUE
 300          ELSE
 301             DO 260 J = 1, N
 302                RTEMP = ZERO
 303                DO 230 L = 1, K
 304                   RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )
 305   230          CONTINUE
 306                IF( BETA.EQ.ZERO ) THEN
 307                   C( J, J ) = ALPHA*RTEMP
 308                ELSE
 309                   C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )
 310                END IF
 311                DO 250 I = J + 1, N
 312                   TEMP = ZERO
 313                   DO 240 L = 1, K
 314                      TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )
 315   240             CONTINUE
 316                   IF( BETA.EQ.ZERO ) THEN
 317                      C( I, J ) = ALPHA*TEMP
 318                   ELSE
 319                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
 320                   END IF
 321   250          CONTINUE
 322   260       CONTINUE
 323          END IF
 324       END IF
 325 *
 326       RETURN
 327 *
 328 *     End of ZHERK .
 329 *
 330       END