BLAS/SRC/zher2k.f

   1 *> \brief \b ZHER2K
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE ZHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  12 *
  13 *       .. Scalar Arguments ..
  14 *       COMPLEX*16 ALPHA
  15 *       DOUBLE PRECISION BETA
  16 *       INTEGER K,LDA,LDB,LDC,N
  17 *       CHARACTER TRANS,UPLO
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *> ZHER2K  performs one of the hermitian rank 2k operations
  30 *>
  31 *>    C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C,
  32 *>
  33 *> or
  34 *>
  35 *>    C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C,
  36 *>
  37 *> where  alpha and beta  are scalars with  beta  real,  C is an  n by n
  38 *> hermitian matrix and  A and B  are  n by k matrices in the first case
  39 *> and  k by n  matrices in the second case.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] UPLO
  46 *> \verbatim
  47 *>          UPLO is CHARACTER*1
  48 *>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  49 *>           triangular  part  of the  array  C  is to be  referenced  as
  50 *>           follows:
  51 *>
  52 *>              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  53 *>                                  is to be referenced.
  54 *>
  55 *>              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  56 *>                                  is to be referenced.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] TRANS
  60 *> \verbatim
  61 *>          TRANS is CHARACTER*1
  62 *>           On entry,  TRANS  specifies the operation to be performed as
  63 *>           follows:
  64 *>
  65 *>              TRANS = 'N' or 'n'    C := alpha*A*B**H          +
  66 *>                                         conjg( alpha )*B*A**H +
  67 *>                                         beta*C.
  68 *>
  69 *>              TRANS = 'C' or 'c'    C := alpha*A**H*B          +
  70 *>                                         conjg( alpha )*B**H*A +
  71 *>                                         beta*C.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] N
  75 *> \verbatim
  76 *>          N is INTEGER
  77 *>           On entry,  N specifies the order of the matrix C.  N must be
  78 *>           at least zero.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] K
  82 *> \verbatim
  83 *>          K is INTEGER
  84 *>           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  85 *>           of  columns  of the  matrices  A and B,  and on  entry  with
  86 *>           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
  87 *>           matrices  A and B.  K must be at least zero.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] ALPHA
  91 *> \verbatim
  92 *>          ALPHA is COMPLEX*16 .
  93 *>           On entry, ALPHA specifies the scalar alpha.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] A
  97 *> \verbatim
  98 *>          A is COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
  99 *>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
 100 *>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
 101 *>           part of the array  A  must contain the matrix  A,  otherwise
 102 *>           the leading  k by n  part of the array  A  must contain  the
 103 *>           matrix A.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] LDA
 107 *> \verbatim
 108 *>          LDA is INTEGER
 109 *>           On entry, LDA specifies the first dimension of A as declared
 110 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 111 *>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
 112 *>           be at least  max( 1, k ).
 113 *> \endverbatim
 114 *>
 115 *> \param[in] B
 116 *> \verbatim
 117 *>          B is COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
 118 *>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
 119 *>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
 120 *>           part of the array  B  must contain the matrix  B,  otherwise
 121 *>           the leading  k by n  part of the array  B  must contain  the
 122 *>           matrix B.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] LDB
 126 *> \verbatim
 127 *>          LDB is INTEGER
 128 *>           On entry, LDB specifies the first dimension of B as declared
 129 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 130 *>           then  LDB must be at least  max( 1, n ), otherwise  LDB must
 131 *>           be at least  max( 1, k ).
 132 *>           Unchanged on exit.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] BETA
 136 *> \verbatim
 137 *>          BETA is DOUBLE PRECISION .
 138 *>           On entry, BETA specifies the scalar beta.
 139 *> \endverbatim
 140 *>
 141 *> \param[in,out] C
 142 *> \verbatim
 143 *>          C is COMPLEX*16 array of DIMENSION ( LDC, n ).
 144 *>           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
 145 *>           upper triangular part of the array C must contain the upper
 146 *>           triangular part  of the  hermitian matrix  and the strictly
 147 *>           lower triangular part of C is not referenced.  On exit, the
 148 *>           upper triangular part of the array  C is overwritten by the
 149 *>           upper triangular part of the updated matrix.
 150 *>           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
 151 *>           lower triangular part of the array C must contain the lower
 152 *>           triangular part  of the  hermitian matrix  and the strictly
 153 *>           upper triangular part of C is not referenced.  On exit, the
 154 *>           lower triangular part of the array  C is overwritten by the
 155 *>           lower triangular part of the updated matrix.
 156 *>           Note that the imaginary parts of the diagonal elements need
 157 *>           not be set,  they are assumed to be zero,  and on exit they
 158 *>           are set to zero.
 159 *> \endverbatim
 160 *>
 161 *> \param[in] LDC
 162 *> \verbatim
 163 *>          LDC is INTEGER
 164 *>           On entry, LDC specifies the first dimension of C as declared
 165 *>           in  the  calling  (sub)  program.   LDC  must  be  at  least
 166 *>           max( 1, n ).
 167 *> \endverbatim
 168 *
 169 *  Authors:
 170 *  ========
 171 *
 172 *> \author Univ. of Tennessee
 173 *> \author Univ. of California Berkeley
 174 *> \author Univ. of Colorado Denver
 175 *> \author NAG Ltd.
 176 *
 177 *> \date November 2011
 178 *
 179 *> \ingroup complex16_blas_level3
 180 *
 181 *> \par Further Details:
 182 *  =====================
 183 *>
 184 *> \verbatim
 185 *>
 186 *>  Level 3 Blas routine.
 187 *>
 188 *>  -- Written on 8-February-1989.
 189 *>     Jack Dongarra, Argonne National Laboratory.
 190 *>     Iain Duff, AERE Harwell.
 191 *>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 192 *>     Sven Hammarling, Numerical Algorithms Group Ltd.
 193 *>
 194 *>  -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
 195 *>     Ed Anderson, Cray Research Inc.
 196 *> \endverbatim
 197 *>
 198 *  =====================================================================
 199       SUBROUTINE ZHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
 200 *
 201 *  -- Reference BLAS level3 routine (version 3.4.0) --
 202 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 203 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 204 *     November 2011
 205 *
 206 *     .. Scalar Arguments ..
 207       COMPLEX*16 ALPHA
 208       DOUBLE PRECISION BETA
 209       INTEGER K,LDA,LDB,LDC,N
 210       CHARACTER TRANS,UPLO
 211 *     ..
 212 *     .. Array Arguments ..
 213       COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
 214 *     ..
 215 *
 216 *  =====================================================================
 217 *
 218 *     .. External Functions ..
 219       LOGICAL LSAME
 220       EXTERNAL LSAME
 221 *     ..
 222 *     .. External Subroutines ..
 223       EXTERNAL XERBLA
 224 *     ..
 225 *     .. Intrinsic Functions ..
 226       INTRINSIC DBLE,DCONJG,MAX
 227 *     ..
 228 *     .. Local Scalars ..
 229       COMPLEX*16 TEMP1,TEMP2
 230       INTEGER I,INFO,J,L,NROWA
 231       LOGICAL UPPER
 232 *     ..
 233 *     .. Parameters ..
 234       DOUBLE PRECISION ONE
 235       PARAMETER (ONE=1.0D+0)
 236       COMPLEX*16 ZERO
 237       PARAMETER (ZERO= (0.0D+0,0.0D+0))
 238 *     ..
 239 *
 240 *     Test the input parameters.
 241 *
 242       IF (LSAME(TRANS,'N')) THEN
 243           NROWA = N
 244       ELSE
 245           NROWA = K
 246       END IF
 247       UPPER = LSAME(UPLO,'U')
 248 *
 249       INFO = 0
 250       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
 251           INFO = 1
 252       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
 253      +         (.NOT.LSAME(TRANS,'C'))) THEN
 254           INFO = 2
 255       ELSE IF (N.LT.0) THEN
 256           INFO = 3
 257       ELSE IF (K.LT.0) THEN
 258           INFO = 4
 259       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
 260           INFO = 7
 261       ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
 262           INFO = 9
 263       ELSE IF (LDC.LT.MAX(1,N)) THEN
 264           INFO = 12
 265       END IF
 266       IF (INFO.NE.0) THEN
 267           CALL XERBLA('ZHER2K',INFO)
 268           RETURN
 269       END IF
 270 *
 271 *     Quick return if possible.
 272 *
 273       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
 274      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
 275 *
 276 *     And when  alpha.eq.zero.
 277 *
 278       IF (ALPHA.EQ.ZERO) THEN
 279           IF (UPPER) THEN
 280               IF (BETA.EQ.DBLE(ZERO)) THEN
 281                   DO 20 J = 1,N
 282                       DO 10 I = 1,J
 283                           C(I,J) = ZERO
 284    10                 CONTINUE
 285    20             CONTINUE
 286               ELSE
 287                   DO 40 J = 1,N
 288                       DO 30 I = 1,J - 1
 289                           C(I,J) = BETA*C(I,J)
 290    30                 CONTINUE
 291                       C(J,J) = BETA*DBLE(C(J,J))
 292    40             CONTINUE
 293               END IF
 294           ELSE
 295               IF (BETA.EQ.DBLE(ZERO)) THEN
 296                   DO 60 J = 1,N
 297                       DO 50 I = J,N
 298                           C(I,J) = ZERO
 299    50                 CONTINUE
 300    60             CONTINUE
 301               ELSE
 302                   DO 80 J = 1,N
 303                       C(J,J) = BETA*DBLE(C(J,J))
 304                       DO 70 I = J + 1,N
 305                           C(I,J) = BETA*C(I,J)
 306    70                 CONTINUE
 307    80             CONTINUE
 308               END IF
 309           END IF
 310           RETURN
 311       END IF
 312 *
 313 *     Start the operations.
 314 *
 315       IF (LSAME(TRANS,'N')) THEN
 316 *
 317 *        Form  C := alpha*A*B**H + conjg( alpha )*B*A**H +
 318 *                   C.
 319 *
 320           IF (UPPER) THEN
 321               DO 130 J = 1,N
 322                   IF (BETA.EQ.DBLE(ZERO)) THEN
 323                       DO 90 I = 1,J
 324                           C(I,J) = ZERO
 325    90                 CONTINUE
 326                   ELSE IF (BETA.NE.ONE) THEN
 327                       DO 100 I = 1,J - 1
 328                           C(I,J) = BETA*C(I,J)
 329   100                 CONTINUE
 330                       C(J,J) = BETA*DBLE(C(J,J))
 331                   ELSE
 332                       C(J,J) = DBLE(C(J,J))
 333                   END IF
 334                   DO 120 L = 1,K
 335                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
 336                           TEMP1 = ALPHA*DCONJG(B(J,L))
 337                           TEMP2 = DCONJG(ALPHA*A(J,L))
 338                           DO 110 I = 1,J - 1
 339                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
 340      +                                 B(I,L)*TEMP2
 341   110                     CONTINUE
 342                           C(J,J) = DBLE(C(J,J)) +
 343      +                             DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2)
 344                       END IF
 345   120             CONTINUE
 346   130         CONTINUE
 347           ELSE
 348               DO 180 J = 1,N
 349                   IF (BETA.EQ.DBLE(ZERO)) THEN
 350                       DO 140 I = J,N
 351                           C(I,J) = ZERO
 352   140                 CONTINUE
 353                   ELSE IF (BETA.NE.ONE) THEN
 354                       DO 150 I = J + 1,N
 355                           C(I,J) = BETA*C(I,J)
 356   150                 CONTINUE
 357                       C(J,J) = BETA*DBLE(C(J,J))
 358                   ELSE
 359                       C(J,J) = DBLE(C(J,J))
 360                   END IF
 361                   DO 170 L = 1,K
 362                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
 363                           TEMP1 = ALPHA*DCONJG(B(J,L))
 364                           TEMP2 = DCONJG(ALPHA*A(J,L))
 365                           DO 160 I = J + 1,N
 366                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
 367      +                                 B(I,L)*TEMP2
 368   160                     CONTINUE
 369                           C(J,J) = DBLE(C(J,J)) +
 370      +                             DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2)
 371                       END IF
 372   170             CONTINUE
 373   180         CONTINUE
 374           END IF
 375       ELSE
 376 *
 377 *        Form  C := alpha*A**H*B + conjg( alpha )*B**H*A +
 378 *                   C.
 379 *
 380           IF (UPPER) THEN
 381               DO 210 J = 1,N
 382                   DO 200 I = 1,J
 383                       TEMP1 = ZERO
 384                       TEMP2 = ZERO
 385                       DO 190 L = 1,K
 386                           TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J)
 387                           TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J)
 388   190                 CONTINUE
 389                       IF (I.EQ.J) THEN
 390                           IF (BETA.EQ.DBLE(ZERO)) THEN
 391                               C(J,J) = DBLE(ALPHA*TEMP1+
 392      +                                 DCONJG(ALPHA)*TEMP2)
 393                           ELSE
 394                               C(J,J) = BETA*DBLE(C(J,J)) +
 395      +                                 DBLE(ALPHA*TEMP1+
 396      +                                 DCONJG(ALPHA)*TEMP2)
 397                           END IF
 398                       ELSE
 399                           IF (BETA.EQ.DBLE(ZERO)) THEN
 400                               C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2
 401                           ELSE
 402                               C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
 403      +                                 DCONJG(ALPHA)*TEMP2
 404                           END IF
 405                       END IF
 406   200             CONTINUE
 407   210         CONTINUE
 408           ELSE
 409               DO 240 J = 1,N
 410                   DO 230 I = J,N
 411                       TEMP1 = ZERO
 412                       TEMP2 = ZERO
 413                       DO 220 L = 1,K
 414                           TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J)
 415                           TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J)
 416   220                 CONTINUE
 417                       IF (I.EQ.J) THEN
 418                           IF (BETA.EQ.DBLE(ZERO)) THEN
 419                               C(J,J) = DBLE(ALPHA*TEMP1+
 420      +                                 DCONJG(ALPHA)*TEMP2)
 421                           ELSE
 422                               C(J,J) = BETA*DBLE(C(J,J)) +
 423      +                                 DBLE(ALPHA*TEMP1+
 424      +                                 DCONJG(ALPHA)*TEMP2)
 425                           END IF
 426                       ELSE
 427                           IF (BETA.EQ.DBLE(ZERO)) THEN
 428                               C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2
 429                           ELSE
 430                               C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
 431      +                                 DCONJG(ALPHA)*TEMP2
 432                           END IF
 433                       END IF
 434   230             CONTINUE
 435   240         CONTINUE
 436           END IF
 437       END IF
 438 *
 439       RETURN
 440 *
 441 *     End of ZHER2K.
 442 *
 443       END