1       SUBROUTINE CHERK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
  2 *     .. Scalar Arguments ..
  3       REAL ALPHA,BETA
  4       INTEGER K,LDA,LDC,N
  5       CHARACTER TRANS,UPLO
  6 *     ..
  7 *     .. Array Arguments ..
  8       COMPLEX A(LDA,*),C(LDC,*)
  9 *     ..
 10 *
 11 *  Purpose
 12 *  =======
 13 *
 14 *  CHERK  performs one of the hermitian rank k operations
 15 *
 16 *     C := alpha*A*A**H + beta*C,
 17 *
 18 *  or
 19 *
 20 *     C := alpha*A**H*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*A**H + beta*C.
 47 *
 48 *              TRANS = 'C' or 'c'   C := alpha*A**H*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  - REAL            .
 65 *           On entry, ALPHA specifies the scalar alpha.
 66 *           Unchanged on exit.
 67 *
 68 *  A      - COMPLEX          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   - REAL            .
 84 *           On entry, BETA specifies the scalar beta.
 85 *           Unchanged on exit.
 86 *
 87 *  C      - COMPLEX          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 *  Further Details
111 *  ===============
112 *
113 *  Level 3 Blas routine.
114 *
115 *  -- Written on 8-February-1989.
116 *     Jack Dongarra, Argonne National Laboratory.
117 *     Iain Duff, AERE Harwell.
118 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
119 *     Sven Hammarling, Numerical Algorithms Group Ltd.
120 *
121 *  -- Modified 8-Nov-93 to set C(J,J) to REAL( C(J,J) ) when BETA = 1.
122 *     Ed Anderson, Cray Research Inc.
123 *
124 *  =====================================================================
125 *
126 *     .. External Functions ..
127       LOGICAL LSAME
128       EXTERNAL LSAME
129 *     ..
130 *     .. External Subroutines ..
131       EXTERNAL XERBLA
132 *     ..
133 *     .. Intrinsic Functions ..
134       INTRINSIC CMPLX,CONJG,MAX,REAL
135 *     ..
136 *     .. Local Scalars ..
137       COMPLEX TEMP
138       REAL RTEMP
139       INTEGER I,INFO,J,L,NROWA
140       LOGICAL UPPER
141 *     ..
142 *     .. Parameters ..
143       REAL ONE,ZERO
144       PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
145 *     ..
146 *
147 *     Test the input parameters.
148 *
149       IF (LSAME(TRANS,'N')) THEN
150           NROWA = N
151       ELSE
152           NROWA = K
153       END IF
154       UPPER = LSAME(UPLO,'U')
155 *
156       INFO = 0
157       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
158           INFO = 1
159       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
160      +         (.NOT.LSAME(TRANS,'C'))) THEN
161           INFO = 2
162       ELSE IF (N.LT.0THEN
163           INFO = 3
164       ELSE IF (K.LT.0THEN
165           INFO = 4
166       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
167           INFO = 7
168       ELSE IF (LDC.LT.MAX(1,N)) THEN
169           INFO = 10
170       END IF
171       IF (INFO.NE.0THEN
172           CALL XERBLA('CHERK ',INFO)
173           RETURN
174       END IF
175 *
176 *     Quick return if possible.
177 *
178       IF ((N.EQ.0.OR. (((ALPHA.EQ.ZERO).OR.
179      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
180 *
181 *     And when  alpha.eq.zero.
182 *
183       IF (ALPHA.EQ.ZERO) THEN
184           IF (UPPER) THEN
185               IF (BETA.EQ.ZERO) THEN
186                   DO 20 J = 1,N
187                       DO 10 I = 1,J
188                           C(I,J) = ZERO
189    10                 CONTINUE
190    20             CONTINUE
191               ELSE
192                   DO 40 J = 1,N
193                       DO 30 I = 1,J - 1
194                           C(I,J) = BETA*C(I,J)
195    30                 CONTINUE
196                       C(J,J) = BETA*REAL(C(J,J))
197    40             CONTINUE
198               END IF
199           ELSE
200               IF (BETA.EQ.ZERO) THEN
201                   DO 60 J = 1,N
202                       DO 50 I = J,N
203                           C(I,J) = ZERO
204    50                 CONTINUE
205    60             CONTINUE
206               ELSE
207                   DO 80 J = 1,N
208                       C(J,J) = BETA*REAL(C(J,J))
209                       DO 70 I = J + 1,N
210                           C(I,J) = BETA*C(I,J)
211    70                 CONTINUE
212    80             CONTINUE
213               END IF
214           END IF
215           RETURN
216       END IF
217 *
218 *     Start the operations.
219 *
220       IF (LSAME(TRANS,'N')) THEN
221 *
222 *        Form  C := alpha*A*A**H + beta*C.
223 *
224           IF (UPPER) THEN
225               DO 130 J = 1,N
226                   IF (BETA.EQ.ZERO) THEN
227                       DO 90 I = 1,J
228                           C(I,J) = ZERO
229    90                 CONTINUE
230                   ELSE IF (BETA.NE.ONE) THEN
231                       DO 100 I = 1,J - 1
232                           C(I,J) = BETA*C(I,J)
233   100                 CONTINUE
234                       C(J,J) = BETA*REAL(C(J,J))
235                   ELSE
236                       C(J,J) = REAL(C(J,J))
237                   END IF
238                   DO 120 L = 1,K
239                       IF (A(J,L).NE.CMPLX(ZERO)) THEN
240                           TEMP = ALPHA*CONJG(A(J,L))
241                           DO 110 I = 1,J - 1
242                               C(I,J) = C(I,J) + TEMP*A(I,L)
243   110                     CONTINUE
244                           C(J,J) = REAL(C(J,J)) + REAL(TEMP*A(I,L))
245                       END IF
246   120             CONTINUE
247   130         CONTINUE
248           ELSE
249               DO 180 J = 1,N
250                   IF (BETA.EQ.ZERO) THEN
251                       DO 140 I = J,N
252                           C(I,J) = ZERO
253   140                 CONTINUE
254                   ELSE IF (BETA.NE.ONE) THEN
255                       C(J,J) = BETA*REAL(C(J,J))
256                       DO 150 I = J + 1,N
257                           C(I,J) = BETA*C(I,J)
258   150                 CONTINUE
259                   ELSE
260                       C(J,J) = REAL(C(J,J))
261                   END IF
262                   DO 170 L = 1,K
263                       IF (A(J,L).NE.CMPLX(ZERO)) THEN
264                           TEMP = ALPHA*CONJG(A(J,L))
265                           C(J,J) = REAL(C(J,J)) + REAL(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*A**H*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 + CONJG(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 + CONJG(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*REAL(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 + CONJG(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*REAL(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 + CONJG(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 CHERK .
329 *
330       END