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