1       SUBROUTINE ZHFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
  2      $                  C )
  3 *
  4 *  -- LAPACK routine (version 3.3.1)                                    --
  5 *
  6 *  -- Contributed by Julien Langou of the Univ. of Colorado Denver    --
  7 *  -- April 2011                                                      --
  8 *
  9 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 10 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 11 *
 12 *     ..
 13 *     .. Scalar Arguments ..
 14       DOUBLE PRECISION   ALPHA, BETA
 15       INTEGER            K, LDA, N
 16       CHARACTER          TRANS, TRANSR, UPLO
 17 *     ..
 18 *     .. Array Arguments ..
 19       COMPLEX*16         A( LDA, * ), C( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  Level 3 BLAS like routine for C in RFP Format.
 26 *
 27 *  ZHFRK performs one of the Hermitian rank--k operations
 28 *
 29 *     C := alpha*A*A**H + beta*C,
 30 *
 31 *  or
 32 *
 33 *     C := alpha*A**H*A + beta*C,
 34 *
 35 *  where alpha and beta are real scalars, C is an n--by--n Hermitian
 36 *  matrix and A is an n--by--k matrix in the first case and a k--by--n
 37 *  matrix in the second case.
 38 *
 39 *  Arguments
 40 *  ==========
 41 *
 42 *  TRANSR  (input) CHARACTER*1
 43 *          = 'N':  The Normal Form of RFP A is stored;
 44 *          = 'C':  The Conjugate-transpose Form of RFP A is stored.
 45 *
 46 *  UPLO    (input) CHARACTER*1
 47 *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
 48 *           triangular  part  of the  array  C  is to be  referenced  as
 49 *           follows:
 50 *
 51 *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
 52 *                                  is to be referenced.
 53 *
 54 *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
 55 *                                  is to be referenced.
 56 *
 57 *           Unchanged on exit.
 58 *
 59 *  TRANS   (input) CHARACTER*1
 60 *           On entry,  TRANS  specifies the operation to be performed as
 61 *           follows:
 62 *
 63 *              TRANS = 'N' or 'n'   C := alpha*A*A**H + beta*C.
 64 *
 65 *              TRANS = 'C' or 'c'   C := alpha*A**H*A + beta*C.
 66 *
 67 *           Unchanged on exit.
 68 *
 69 *  N       (input) INTEGER
 70 *           On entry,  N specifies the order of the matrix C.  N must be
 71 *           at least zero.
 72 *           Unchanged on exit.
 73 *
 74 *  K       (input) INTEGER
 75 *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
 76 *           of  columns   of  the   matrix   A,   and  on   entry   with
 77 *           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
 78 *           matrix A.  K must be at least zero.
 79 *           Unchanged on exit.
 80 *
 81 *  ALPHA   (input) DOUBLE PRECISION
 82 *           On entry, ALPHA specifies the scalar alpha.
 83 *           Unchanged on exit.
 84 *
 85 *  A       (input) COMPLEX*16 array of DIMENSION (LDA,ka)
 86 *           where KA
 87 *           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
 88 *           entry with TRANS = 'N' or 'n', the leading N--by--K part of
 89 *           the array A must contain the matrix A, otherwise the leading
 90 *           K--by--N part of the array A must contain the matrix A.
 91 *           Unchanged on exit.
 92 *
 93 *  LDA     (input) INTEGER
 94 *           On entry, LDA specifies the first dimension of A as declared
 95 *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 96 *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
 97 *           be at least  max( 1, k ).
 98 *           Unchanged on exit.
 99 *
100 *  BETA    (input) DOUBLE PRECISION
101 *           On entry, BETA specifies the scalar beta.
102 *           Unchanged on exit.
103 *
104 *  C       (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
105 *           On entry, the matrix A in RFP Format. RFP Format is
106 *           described by TRANSR, UPLO and N. Note that the imaginary
107 *           parts of the diagonal elements need not be set, they are
108 *           assumed to be zero, and on exit they are set to zero.
109 *
110 *  =====================================================================
111 *
112 *     .. Parameters ..
113       DOUBLE PRECISION   ONE, ZERO
114       COMPLEX*16         CZERO
115       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
116       PARAMETER          ( CZERO = ( 0.0D+00.0D+0 ) )
117 *     ..
118 *     .. Local Scalars ..
119       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
120       INTEGER            INFO, NROWA, J, NK, N1, N2
121       COMPLEX*16         CALPHA, CBETA
122 *     ..
123 *     .. External Functions ..
124       LOGICAL            LSAME
125       EXTERNAL           LSAME
126 *     ..
127 *     .. External Subroutines ..
128       EXTERNAL           XERBLA, ZGEMM, ZHERK
129 *     ..
130 *     .. Intrinsic Functions ..
131       INTRINSIC          MAXDCMPLX
132 *     ..
133 *     .. Executable Statements ..
134 *
135 *
136 *     Test the input parameters.
137 *
138       INFO = 0
139       NORMALTRANSR = LSAME( TRANSR, 'N' )
140       LOWER = LSAME( UPLO, 'L' )
141       NOTRANS = LSAME( TRANS, 'N' )
142 *
143       IF( NOTRANS ) THEN
144          NROWA = N
145       ELSE
146          NROWA = K
147       END IF
148 *
149       IF.NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
150          INFO = -1
151       ELSE IF.NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
152          INFO = -2
153       ELSE IF.NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
154          INFO = -3
155       ELSE IF( N.LT.0 ) THEN
156          INFO = -4
157       ELSE IF( K.LT.0 ) THEN
158          INFO = -5
159       ELSE IF( LDA.LT.MAX1, NROWA ) ) THEN
160          INFO = -8
161       END IF
162       IF( INFO.NE.0 ) THEN
163          CALL XERBLA( 'ZHFRK '-INFO )
164          RETURN
165       END IF
166 *
167 *     Quick return if possible.
168 *
169 *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
170 *     done (it is in ZHERK for example) and left in the general case.
171 *
172       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
173      $    ( BETA.EQ.ONE ) ) )RETURN
174 *
175       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
176          DO J = 1, ( ( N*( N+1 ) ) / 2 )
177             C( J ) = CZERO
178          END DO
179          RETURN
180       END IF
181 *
182       CALPHA = DCMPLX( ALPHA, ZERO )
183       CBETA = DCMPLX( BETA, ZERO )
184 *
185 *     C is N-by-N.
186 *     If N is odd, set NISODD = .TRUE., and N1 and N2.
187 *     If N is even, NISODD = .FALSE., and NK.
188 *
189       IFMOD( N, 2 ).EQ.0 ) THEN
190          NISODD = .FALSE.
191          NK = N / 2
192       ELSE
193          NISODD = .TRUE.
194          IF( LOWER ) THEN
195             N2 = N / 2
196             N1 = N - N2
197          ELSE
198             N1 = N / 2
199             N2 = N - N1
200          END IF
201       END IF
202 *
203       IF( NISODD ) THEN
204 *
205 *        N is odd
206 *
207          IF( NORMALTRANSR ) THEN
208 *
209 *           N is odd and TRANSR = 'N'
210 *
211             IF( LOWER ) THEN
212 *
213 *              N is odd, TRANSR = 'N', and UPLO = 'L'
214 *
215                IF( NOTRANS ) THEN
216 *
217 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
218 *
219                   CALL ZHERK( 'L''N', N1, K, ALPHA, A( 11 ), LDA,
220      $                        BETA, C( 1 ), N )
221                   CALL ZHERK( 'U''N', N2, K, ALPHA, A( N1+11 ), LDA,
222      $                        BETA, C( N+1 ), N )
223                   CALL ZGEMM( 'N''C', N2, N1, K, CALPHA, A( N1+11 ),
224      $                        LDA, A( 11 ), LDA, CBETA, C( N1+1 ), N )
225 *
226                ELSE
227 *
228 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
229 *
230                   CALL ZHERK( 'L''C', N1, K, ALPHA, A( 11 ), LDA,
231      $                        BETA, C( 1 ), N )
232                   CALL ZHERK( 'U''C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
233      $                        BETA, C( N+1 ), N )
234                   CALL ZGEMM( 'C''N', N2, N1, K, CALPHA, A( 1, N1+1 ),
235      $                        LDA, A( 11 ), LDA, CBETA, C( N1+1 ), N )
236 *
237                END IF
238 *
239             ELSE
240 *
241 *              N is odd, TRANSR = 'N', and UPLO = 'U'
242 *
243                IF( NOTRANS ) THEN
244 *
245 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
246 *
247                   CALL ZHERK( 'L''N', N1, K, ALPHA, A( 11 ), LDA,
248      $                        BETA, C( N2+1 ), N )
249                   CALL ZHERK( 'U''N', N2, K, ALPHA, A( N2, 1 ), LDA,
250      $                        BETA, C( N1+1 ), N )
251                   CALL ZGEMM( 'N''C', N1, N2, K, CALPHA, A( 11 ),
252      $                        LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N )
253 *
254                ELSE
255 *
256 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
257 *
258                   CALL ZHERK( 'L''C', N1, K, ALPHA, A( 11 ), LDA,
259      $                        BETA, C( N2+1 ), N )
260                   CALL ZHERK( 'U''C', N2, K, ALPHA, A( 1, N2 ), LDA,
261      $                        BETA, C( N1+1 ), N )
262                   CALL ZGEMM( 'C''N', N1, N2, K, CALPHA, A( 11 ),
263      $                        LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N )
264 *
265                END IF
266 *
267             END IF
268 *
269          ELSE
270 *
271 *           N is odd, and TRANSR = 'C'
272 *
273             IF( LOWER ) THEN
274 *
275 *              N is odd, TRANSR = 'C', and UPLO = 'L'
276 *
277                IF( NOTRANS ) THEN
278 *
279 *                 N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
280 *
281                   CALL ZHERK( 'U''N', N1, K, ALPHA, A( 11 ), LDA,
282      $                        BETA, C( 1 ), N1 )
283                   CALL ZHERK( 'L''N', N2, K, ALPHA, A( N1+11 ), LDA,
284      $                        BETA, C( 2 ), N1 )
285                   CALL ZGEMM( 'N''C', N1, N2, K, CALPHA, A( 11 ),
286      $                        LDA, A( N1+11 ), LDA, CBETA,
287      $                        C( N1*N1+1 ), N1 )
288 *
289                ELSE
290 *
291 *                 N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
292 *
293                   CALL ZHERK( 'U''C', N1, K, ALPHA, A( 11 ), LDA,
294      $                        BETA, C( 1 ), N1 )
295                   CALL ZHERK( 'L''C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
296      $                        BETA, C( 2 ), N1 )
297                   CALL ZGEMM( 'C''N', N1, N2, K, CALPHA, A( 11 ),
298      $                        LDA, A( 1, N1+1 ), LDA, CBETA,
299      $                        C( N1*N1+1 ), N1 )
300 *
301                END IF
302 *
303             ELSE
304 *
305 *              N is odd, TRANSR = 'C', and UPLO = 'U'
306 *
307                IF( NOTRANS ) THEN
308 *
309 *                 N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
310 *
311                   CALL ZHERK( 'U''N', N1, K, ALPHA, A( 11 ), LDA,
312      $                        BETA, C( N2*N2+1 ), N2 )
313                   CALL ZHERK( 'L''N', N2, K, ALPHA, A( N1+11 ), LDA,
314      $                        BETA, C( N1*N2+1 ), N2 )
315                   CALL ZGEMM( 'N''C', N2, N1, K, CALPHA, A( N1+11 ),
316      $                        LDA, A( 11 ), LDA, CBETA, C( 1 ), N2 )
317 *
318                ELSE
319 *
320 *                 N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
321 *
322                   CALL ZHERK( 'U''C', N1, K, ALPHA, A( 11 ), LDA,
323      $                        BETA, C( N2*N2+1 ), N2 )
324                   CALL ZHERK( 'L''C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
325      $                        BETA, C( N1*N2+1 ), N2 )
326                   CALL ZGEMM( 'C''N', N2, N1, K, CALPHA, A( 1, N1+1 ),
327      $                        LDA, A( 11 ), LDA, CBETA, C( 1 ), N2 )
328 *
329                END IF
330 *
331             END IF
332 *
333          END IF
334 *
335       ELSE
336 *
337 *        N is even
338 *
339          IF( NORMALTRANSR ) THEN
340 *
341 *           N is even and TRANSR = 'N'
342 *
343             IF( LOWER ) THEN
344 *
345 *              N is even, TRANSR = 'N', and UPLO = 'L'
346 *
347                IF( NOTRANS ) THEN
348 *
349 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
350 *
351                   CALL ZHERK( 'L''N', NK, K, ALPHA, A( 11 ), LDA,
352      $                        BETA, C( 2 ), N+1 )
353                   CALL ZHERK( 'U''N', NK, K, ALPHA, A( NK+11 ), LDA,
354      $                        BETA, C( 1 ), N+1 )
355                   CALL ZGEMM( 'N''C', NK, NK, K, CALPHA, A( NK+11 ),
356      $                        LDA, A( 11 ), LDA, CBETA, C( NK+2 ),
357      $                        N+1 )
358 *
359                ELSE
360 *
361 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
362 *
363                   CALL ZHERK( 'L''C', NK, K, ALPHA, A( 11 ), LDA,
364      $                        BETA, C( 2 ), N+1 )
365                   CALL ZHERK( 'U''C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
366      $                        BETA, C( 1 ), N+1 )
367                   CALL ZGEMM( 'C''N', NK, NK, K, CALPHA, A( 1, NK+1 ),
368      $                        LDA, A( 11 ), LDA, CBETA, C( NK+2 ),
369      $                        N+1 )
370 *
371                END IF
372 *
373             ELSE
374 *
375 *              N is even, TRANSR = 'N', and UPLO = 'U'
376 *
377                IF( NOTRANS ) THEN
378 *
379 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
380 *
381                   CALL ZHERK( 'L''N', NK, K, ALPHA, A( 11 ), LDA,
382      $                        BETA, C( NK+2 ), N+1 )
383                   CALL ZHERK( 'U''N', NK, K, ALPHA, A( NK+11 ), LDA,
384      $                        BETA, C( NK+1 ), N+1 )
385                   CALL ZGEMM( 'N''C', NK, NK, K, CALPHA, A( 11 ),
386      $                        LDA, A( NK+11 ), LDA, CBETA, C( 1 ),
387      $                        N+1 )
388 *
389                ELSE
390 *
391 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
392 *
393                   CALL ZHERK( 'L''C', NK, K, ALPHA, A( 11 ), LDA,
394      $                        BETA, C( NK+2 ), N+1 )
395                   CALL ZHERK( 'U''C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
396      $                        BETA, C( NK+1 ), N+1 )
397                   CALL ZGEMM( 'C''N', NK, NK, K, CALPHA, A( 11 ),
398      $                        LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ),
399      $                        N+1 )
400 *
401                END IF
402 *
403             END IF
404 *
405          ELSE
406 *
407 *           N is even, and TRANSR = 'C'
408 *
409             IF( LOWER ) THEN
410 *
411 *              N is even, TRANSR = 'C', and UPLO = 'L'
412 *
413                IF( NOTRANS ) THEN
414 *
415 *                 N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
416 *
417                   CALL ZHERK( 'U''N', NK, K, ALPHA, A( 11 ), LDA,
418      $                        BETA, C( NK+1 ), NK )
419                   CALL ZHERK( 'L''N', NK, K, ALPHA, A( NK+11 ), LDA,
420      $                        BETA, C( 1 ), NK )
421                   CALL ZGEMM( 'N''C', NK, NK, K, CALPHA, A( 11 ),
422      $                        LDA, A( NK+11 ), LDA, CBETA,
423      $                        C( ( ( NK+1 )*NK )+1 ), NK )
424 *
425                ELSE
426 *
427 *                 N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
428 *
429                   CALL ZHERK( 'U''C', NK, K, ALPHA, A( 11 ), LDA,
430      $                        BETA, C( NK+1 ), NK )
431                   CALL ZHERK( 'L''C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
432      $                        BETA, C( 1 ), NK )
433                   CALL ZGEMM( 'C''N', NK, NK, K, CALPHA, A( 11 ),
434      $                        LDA, A( 1, NK+1 ), LDA, CBETA,
435      $                        C( ( ( NK+1 )*NK )+1 ), NK )
436 *
437                END IF
438 *
439             ELSE
440 *
441 *              N is even, TRANSR = 'C', and UPLO = 'U'
442 *
443                IF( NOTRANS ) THEN
444 *
445 *                 N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
446 *
447                   CALL ZHERK( 'U''N', NK, K, ALPHA, A( 11 ), LDA,
448      $                        BETA, C( NK*( NK+1 )+1 ), NK )
449                   CALL ZHERK( 'L''N', NK, K, ALPHA, A( NK+11 ), LDA,
450      $                        BETA, C( NK*NK+1 ), NK )
451                   CALL ZGEMM( 'N''C', NK, NK, K, CALPHA, A( NK+11 ),
452      $                        LDA, A( 11 ), LDA, CBETA, C( 1 ), NK )
453 *
454                ELSE
455 *
456 *                 N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
457 *
458                   CALL ZHERK( 'U''C', NK, K, ALPHA, A( 11 ), LDA,
459      $                        BETA, C( NK*( NK+1 )+1 ), NK )
460                   CALL ZHERK( 'L''C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
461      $                        BETA, C( NK*NK+1 ), NK )
462                   CALL ZGEMM( 'C''N', NK, NK, K, CALPHA, A( 1, NK+1 ),
463      $                        LDA, A( 11 ), LDA, CBETA, C( 1 ), NK )
464 *
465                END IF
466 *
467             END IF
468 *
469          END IF
470 *
471       END IF
472 *
473       RETURN
474 *
475 *     End of ZHFRK
476 *
477       END