src/BLAS/zherk.f

   1       SUBROUTINE ZHERK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
   2 *     .. Scalar Arguments ..
   3       DOUBLE PRECISION ALPHA,BETA
   4       INTEGER K,LDA,LDC,N
   5       CHARACTER TRANS,UPLO
   6 *     ..
   7 *     .. Array Arguments ..
   8       DOUBLE COMPLEX 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 *  Arguments
  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       DOUBLE COMPLEX TEMP
 135       DOUBLE PRECISION RTEMP
 136       INTEGER I,INFO,J,L,NROWA
 137       LOGICAL UPPER
 138 *     ..
 139 *     .. Parameters ..
 140       DOUBLE PRECISION ONE,ZERO
 141       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 142 *     ..
 143 *
 144 *     Test the input parameters.
 145 *
 146       IF (LSAME(TRANS,'N')) THEN
 147           NROWA = N
 148       ELSE
 149           NROWA = K
 150       END IF
 151       UPPER = LSAME(UPLO,'U')
 152 *
 153       INFO = 0
 154       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
 155           INFO = 1
 156       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
 157      +         (.NOT.LSAME(TRANS,'C'))) THEN
 158           INFO = 2
 159       ELSE IF (N.LT.0) THEN
 160           INFO = 3
 161       ELSE IF (K.LT.0) THEN
 162           INFO = 4
 163       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
 164           INFO = 7
 165       ELSE IF (LDC.LT.MAX(1,N)) THEN
 166           INFO = 10
 167       END IF
 168       IF (INFO.NE.0) THEN
 169           CALL XERBLA('ZHERK ',INFO)
 170           RETURN
 171       END IF
 172 *
 173 *     Quick return if possible.
 174 *
 175       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
 176      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
 177 *
 178 *     And when  alpha.eq.zero.
 179 *
 180       IF (ALPHA.EQ.ZERO) THEN
 181           IF (UPPER) THEN
 182               IF (BETA.EQ.ZERO) THEN
 183                   DO 20 J = 1,N
 184                       DO 10 I = 1,J
 185                           C(I,J) = ZERO
 186    10                 CONTINUE
 187    20             CONTINUE
 188               ELSE
 189                   DO 40 J = 1,N
 190                       DO 30 I = 1,J - 1
 191                           C(I,J) = BETA*C(I,J)
 192    30                 CONTINUE
 193                       C(J,J) = BETA*DBLE(C(J,J))
 194    40             CONTINUE
 195               END IF
 196           ELSE
 197               IF (BETA.EQ.ZERO) THEN
 198                   DO 60 J = 1,N
 199                       DO 50 I = J,N
 200                           C(I,J) = ZERO
 201    50                 CONTINUE
 202    60             CONTINUE
 203               ELSE
 204                   DO 80 J = 1,N
 205                       C(J,J) = BETA*DBLE(C(J,J))
 206                       DO 70 I = J + 1,N
 207                           C(I,J) = BETA*C(I,J)
 208    70                 CONTINUE
 209    80             CONTINUE
 210               END IF
 211           END IF
 212           RETURN
 213       END IF
 214 *
 215 *     Start the operations.
 216 *
 217       IF (LSAME(TRANS,'N')) THEN
 218 *
 219 *        Form  C := alpha*A*conjg( A' ) + beta*C.
 220 *
 221           IF (UPPER) THEN
 222               DO 130 J = 1,N
 223                   IF (BETA.EQ.ZERO) THEN
 224                       DO 90 I = 1,J
 225                           C(I,J) = ZERO
 226    90                 CONTINUE
 227                   ELSE IF (BETA.NE.ONE) THEN
 228                       DO 100 I = 1,J - 1
 229                           C(I,J) = BETA*C(I,J)
 230   100                 CONTINUE
 231                       C(J,J) = BETA*DBLE(C(J,J))
 232                   ELSE
 233                       C(J,J) = DBLE(C(J,J))
 234                   END IF
 235                   DO 120 L = 1,K
 236                       IF (A(J,L).NE.DCMPLX(ZERO)) THEN
 237                           TEMP = ALPHA*DCONJG(A(J,L))
 238                           DO 110 I = 1,J - 1
 239                               C(I,J) = C(I,J) + TEMP*A(I,L)
 240   110                     CONTINUE
 241                           C(J,J) = DBLE(C(J,J)) + DBLE(TEMP*A(I,L))
 242                       END IF
 243   120             CONTINUE
 244   130         CONTINUE
 245           ELSE
 246               DO 180 J = 1,N
 247                   IF (BETA.EQ.ZERO) THEN
 248                       DO 140 I = J,N
 249                           C(I,J) = ZERO
 250   140                 CONTINUE
 251                   ELSE IF (BETA.NE.ONE) THEN
 252                       C(J,J) = BETA*DBLE(C(J,J))
 253                       DO 150 I = J + 1,N
 254                           C(I,J) = BETA*C(I,J)
 255   150                 CONTINUE
 256                   ELSE
 257                       C(J,J) = DBLE(C(J,J))
 258                   END IF
 259                   DO 170 L = 1,K
 260                       IF (A(J,L).NE.DCMPLX(ZERO)) THEN
 261                           TEMP = ALPHA*DCONJG(A(J,L))
 262                           C(J,J) = DBLE(C(J,J)) + DBLE(TEMP*A(J,L))
 263                           DO 160 I = J + 1,N
 264                               C(I,J) = C(I,J) + TEMP*A(I,L)
 265   160                     CONTINUE
 266                       END IF
 267   170             CONTINUE
 268   180         CONTINUE
 269           END IF
 270       ELSE
 271 *
 272 *        Form  C := alpha*conjg( A' )*A + beta*C.
 273 *
 274           IF (UPPER) THEN
 275               DO 220 J = 1,N
 276                   DO 200 I = 1,J - 1
 277                       TEMP = ZERO
 278                       DO 190 L = 1,K
 279                           TEMP = TEMP + DCONJG(A(L,I))*A(L,J)
 280   190                 CONTINUE
 281                       IF (BETA.EQ.ZERO) THEN
 282                           C(I,J) = ALPHA*TEMP
 283                       ELSE
 284                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
 285                       END IF
 286   200             CONTINUE
 287                   RTEMP = ZERO
 288                   DO 210 L = 1,K
 289                       RTEMP = RTEMP + DCONJG(A(L,J))*A(L,J)
 290   210             CONTINUE
 291                   IF (BETA.EQ.ZERO) THEN
 292                       C(J,J) = ALPHA*RTEMP
 293                   ELSE
 294                       C(J,J) = ALPHA*RTEMP + BETA*DBLE(C(J,J))
 295                   END IF
 296   220         CONTINUE
 297           ELSE
 298               DO 260 J = 1,N
 299                   RTEMP = ZERO
 300                   DO 230 L = 1,K
 301                       RTEMP = RTEMP + DCONJG(A(L,J))*A(L,J)
 302   230             CONTINUE
 303                   IF (BETA.EQ.ZERO) THEN
 304                       C(J,J) = ALPHA*RTEMP
 305                   ELSE
 306                       C(J,J) = ALPHA*RTEMP + BETA*DBLE(C(J,J))
 307                   END IF
 308                   DO 250 I = J + 1,N
 309                       TEMP = ZERO
 310                       DO 240 L = 1,K
 311                           TEMP = TEMP + DCONJG(A(L,I))*A(L,J)
 312   240                 CONTINUE
 313                       IF (BETA.EQ.ZERO) THEN
 314                           C(I,J) = ALPHA*TEMP
 315                       ELSE
 316                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
 317                       END IF
 318   250             CONTINUE
 319   260         CONTINUE
 320           END IF
 321       END IF
 322 *
 323       RETURN
 324 *
 325 *     End of ZHERK .
 326 *
 327       END