001:       SUBROUTINE CHFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
002:      +                  C )
003: *
004: *  -- LAPACK routine (version 3.2)                                    --
005: *
006: *  -- Contributed by Julien Langou of the Univ. of Colorado Denver    --
007: *  -- November 2008                                                   --
008: *
009: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
010: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
011: *
012: *     ..
013: *     .. Scalar Arguments ..
014:       REAL               ALPHA, BETA
015:       INTEGER            K, LDA, N
016:       CHARACTER          TRANS, TRANSR, UPLO
017: *     ..
018: *     .. Array Arguments ..
019:       COMPLEX            A( LDA, * ), C( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  Level 3 BLAS like routine for C in RFP Format.
026: *
027: *  CHFRK performs one of the Hermitian rank--k operations
028: *
029: *     C := alpha*A*conjg( A' ) + beta*C,
030: *
031: *  or
032: *
033: *     C := alpha*conjg( A' )*A + beta*C,
034: *
035: *  where alpha and beta are real scalars, C is an n--by--n Hermitian
036: *  matrix and A is an n--by--k matrix in the first case and a k--by--n
037: *  matrix in the second case.
038: *
039: *  Arguments
040: *  ==========
041: *
042: *  TRANSR - (input) CHARACTER.
043: *          = 'N':  The Normal Form of RFP A is stored;
044: *          = 'C':  The Conjugate-transpose Form of RFP A is stored.
045: *
046: *  UPLO   - (input) CHARACTER.
047: *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
048: *           triangular  part  of the  array  C  is to be  referenced  as
049: *           follows:
050: *
051: *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
052: *                                  is to be referenced.
053: *
054: *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
055: *                                  is to be referenced.
056: *
057: *           Unchanged on exit.
058: *
059: *  TRANS  - (input) CHARACTER.
060: *           On entry,  TRANS  specifies the operation to be performed as
061: *           follows:
062: *
063: *              TRANS = 'N' or 'n'   C := alpha*A*conjg( A' ) + beta*C.
064: *
065: *              TRANS = 'C' or 'c'   C := alpha*conjg( A' )*A + beta*C.
066: *
067: *           Unchanged on exit.
068: *
069: *  N      - (input) INTEGER.
070: *           On entry,  N specifies the order of the matrix C.  N must be
071: *           at least zero.
072: *           Unchanged on exit.
073: *
074: *  K      - (input) INTEGER.
075: *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
076: *           of  columns   of  the   matrix   A,   and  on   entry   with
077: *           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
078: *           matrix A.  K must be at least zero.
079: *           Unchanged on exit.
080: *
081: *  ALPHA  - (input) REAL.
082: *           On entry, ALPHA specifies the scalar alpha.
083: *           Unchanged on exit.
084: *
085: *  A      - (input) COMPLEX array of DIMENSION ( LDA, ka ), where KA
086: *           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
087: *           entry with TRANS = 'N' or 'n', the leading N--by--K part of
088: *           the array A must contain the matrix A, otherwise the leading
089: *           K--by--N part of the array A must contain the matrix A.
090: *           Unchanged on exit.
091: *
092: *  LDA    - (input) INTEGER.
093: *           On entry, LDA specifies the first dimension of A as declared
094: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
095: *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
096: *           be at least  max( 1, k ).
097: *           Unchanged on exit.
098: *
099: *  BETA   - (input) REAL.
100: *           On entry, BETA specifies the scalar beta.
101: *           Unchanged on exit.
102: *
103: *  C      - (input/output) COMPLEX array, dimension ( N*(N+1)/2 ).
104: *           On entry, the matrix A in RFP Format. RFP Format is
105: *           described by TRANSR, UPLO and N. Note that the imaginary
106: *           parts of the diagonal elements need not be set, they are
107: *           assumed to be zero, and on exit they are set to zero.
108: *
109: *  Arguments
110: *  ==========
111: *
112: *     ..
113: *     .. Parameters ..
114:       REAL               ONE, ZERO
115:       COMPLEX            CZERO
116:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
117:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
118: *     ..
119: *     .. Local Scalars ..
120:       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
121:       INTEGER            INFO, NROWA, J, NK, N1, N2
122:       COMPLEX            CALPHA, CBETA
123: *     ..
124: *     .. External Functions ..
125:       LOGICAL            LSAME
126:       EXTERNAL           LSAME
127: *     ..
128: *     .. External Subroutines ..
129:       EXTERNAL           CGEMM, CHERK, XERBLA
130: *     ..
131: *     .. Intrinsic Functions ..
132:       INTRINSIC          MAX, CMPLX
133: *     ..
134: *     .. Executable Statements ..
135: *
136: *
137: *     Test the input parameters.
138: *
139:       INFO = 0
140:       NORMALTRANSR = LSAME( TRANSR, 'N' )
141:       LOWER = LSAME( UPLO, 'L' )
142:       NOTRANS = LSAME( TRANS, 'N' )
143: *
144:       IF( NOTRANS ) THEN
145:          NROWA = N
146:       ELSE
147:          NROWA = K
148:       END IF
149: *
150:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
151:          INFO = -1
152:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
153:          INFO = -2
154:       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
155:          INFO = -3
156:       ELSE IF( N.LT.0 ) THEN
157:          INFO = -4
158:       ELSE IF( K.LT.0 ) THEN
159:          INFO = -5
160:       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
161:          INFO = -8
162:       END IF
163:       IF( INFO.NE.0 ) THEN
164:          CALL XERBLA( 'CHFRK ', -INFO )
165:          RETURN
166:       END IF
167: *
168: *     Quick return if possible.
169: *
170: *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
171: *     done (it is in CHERK for example) and left in the general case.
172: *
173:       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
174:      +    ( BETA.EQ.ONE ) ) )RETURN
175: *
176:       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
177:          DO J = 1, ( ( N*( N+1 ) ) / 2 )
178:             C( J ) = CZERO
179:          END DO
180:          RETURN
181:       END IF
182: *
183:       CALPHA = CMPLX( ALPHA, ZERO )
184:       CBETA = CMPLX( BETA, ZERO )
185: *
186: *     C is N-by-N.
187: *     If N is odd, set NISODD = .TRUE., and N1 and N2.
188: *     If N is even, NISODD = .FALSE., and NK.
189: *
190:       IF( MOD( N, 2 ).EQ.0 ) THEN
191:          NISODD = .FALSE.
192:          NK = N / 2
193:       ELSE
194:          NISODD = .TRUE.
195:          IF( LOWER ) THEN
196:             N2 = N / 2
197:             N1 = N - N2
198:          ELSE
199:             N1 = N / 2
200:             N2 = N - N1
201:          END IF
202:       END IF
203: *
204:       IF( NISODD ) THEN
205: *
206: *        N is odd
207: *
208:          IF( NORMALTRANSR ) THEN
209: *
210: *           N is odd and TRANSR = 'N'
211: *
212:             IF( LOWER ) THEN
213: *
214: *              N is odd, TRANSR = 'N', and UPLO = 'L'
215: *
216:                IF( NOTRANS ) THEN
217: *
218: *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
219: *
220:                   CALL CHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
221:      +                        BETA, C( 1 ), N )
222:                   CALL CHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
223:      +                        BETA, C( N+1 ), N )
224:                   CALL CGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
225:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N )
226: *
227:                ELSE
228: *
229: *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
230: *
231:                   CALL CHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
232:      +                        BETA, C( 1 ), N )
233:                   CALL CHERK( 'U', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
234:      +                        BETA, C( N+1 ), N )
235:                   CALL CGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ),
236:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N )
237: *
238:                END IF
239: *
240:             ELSE
241: *
242: *              N is odd, TRANSR = 'N', and UPLO = 'U'
243: *
244:                IF( NOTRANS ) THEN
245: *
246: *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
247: *
248:                   CALL CHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
249:      +                        BETA, C( N2+1 ), N )
250:                   CALL CHERK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
251:      +                        BETA, C( N1+1 ), N )
252:                   CALL CGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ),
253:      +                        LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N )
254: *
255:                ELSE
256: *
257: *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
258: *
259:                   CALL CHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
260:      +                        BETA, C( N2+1 ), N )
261:                   CALL CHERK( 'U', 'C', N2, K, ALPHA, A( 1, N2 ), LDA,
262:      +                        BETA, C( N1+1 ), N )
263:                   CALL CGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ),
264:      +                        LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N )
265: *
266:                END IF
267: *
268:             END IF
269: *
270:          ELSE
271: *
272: *           N is odd, and TRANSR = 'C'
273: *
274:             IF( LOWER ) THEN
275: *
276: *              N is odd, TRANSR = 'C', and UPLO = 'L'
277: *
278:                IF( NOTRANS ) THEN
279: *
280: *                 N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
281: *
282:                   CALL CHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
283:      +                        BETA, C( 1 ), N1 )
284:                   CALL CHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
285:      +                        BETA, C( 2 ), N1 )
286:                   CALL CGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ),
287:      +                        LDA, A( N1+1, 1 ), LDA, CBETA,
288:      +                        C( N1*N1+1 ), N1 )
289: *
290:                ELSE
291: *
292: *                 N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
293: *
294:                   CALL CHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
295:      +                        BETA, C( 1 ), N1 )
296:                   CALL CHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
297:      +                        BETA, C( 2 ), N1 )
298:                   CALL CGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ),
299:      +                        LDA, A( 1, N1+1 ), LDA, CBETA,
300:      +                        C( N1*N1+1 ), N1 )
301: *
302:                END IF
303: *
304:             ELSE
305: *
306: *              N is odd, TRANSR = 'C', and UPLO = 'U'
307: *
308:                IF( NOTRANS ) THEN
309: *
310: *                 N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
311: *
312:                   CALL CHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
313:      +                        BETA, C( N2*N2+1 ), N2 )
314:                   CALL CHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
315:      +                        BETA, C( N1*N2+1 ), N2 )
316:                   CALL CGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
317:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 )
318: *
319:                ELSE
320: *
321: *                 N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
322: *
323:                   CALL CHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
324:      +                        BETA, C( N2*N2+1 ), N2 )
325:                   CALL CHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
326:      +                        BETA, C( N1*N2+1 ), N2 )
327:                   CALL CGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ),
328:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 )
329: *
330:                END IF
331: *
332:             END IF
333: *
334:          END IF
335: *
336:       ELSE
337: *
338: *        N is even
339: *
340:          IF( NORMALTRANSR ) THEN
341: *
342: *           N is even and TRANSR = 'N'
343: *
344:             IF( LOWER ) THEN
345: *
346: *              N is even, TRANSR = 'N', and UPLO = 'L'
347: *
348:                IF( NOTRANS ) THEN
349: *
350: *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
351: *
352:                   CALL CHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
353:      +                        BETA, C( 2 ), N+1 )
354:                   CALL CHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
355:      +                        BETA, C( 1 ), N+1 )
356:                   CALL CGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
357:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ),
358:      +                        N+1 )
359: *
360:                ELSE
361: *
362: *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
363: *
364:                   CALL CHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
365:      +                        BETA, C( 2 ), N+1 )
366:                   CALL CHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
367:      +                        BETA, C( 1 ), N+1 )
368:                   CALL CGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ),
369:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ),
370:      +                        N+1 )
371: *
372:                END IF
373: *
374:             ELSE
375: *
376: *              N is even, TRANSR = 'N', and UPLO = 'U'
377: *
378:                IF( NOTRANS ) THEN
379: *
380: *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
381: *
382:                   CALL CHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
383:      +                        BETA, C( NK+2 ), N+1 )
384:                   CALL CHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
385:      +                        BETA, C( NK+1 ), N+1 )
386:                   CALL CGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
387:      +                        LDA, A( NK+1, 1 ), LDA, CBETA, C( 1 ),
388:      +                        N+1 )
389: *
390:                ELSE
391: *
392: *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
393: *
394:                   CALL CHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
395:      +                        BETA, C( NK+2 ), N+1 )
396:                   CALL CHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
397:      +                        BETA, C( NK+1 ), N+1 )
398:                   CALL CGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ),
399:      +                        LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ),
400:      +                        N+1 )
401: *
402:                END IF
403: *
404:             END IF
405: *
406:          ELSE
407: *
408: *           N is even, and TRANSR = 'C'
409: *
410:             IF( LOWER ) THEN
411: *
412: *              N is even, TRANSR = 'C', and UPLO = 'L'
413: *
414:                IF( NOTRANS ) THEN
415: *
416: *                 N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
417: *
418:                   CALL CHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
419:      +                        BETA, C( NK+1 ), NK )
420:                   CALL CHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
421:      +                        BETA, C( 1 ), NK )
422:                   CALL CGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
423:      +                        LDA, A( NK+1, 1 ), LDA, CBETA,
424:      +                        C( ( ( NK+1 )*NK )+1 ), NK )
425: *
426:                ELSE
427: *
428: *                 N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
429: *
430:                   CALL CHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
431:      +                        BETA, C( NK+1 ), NK )
432:                   CALL CHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
433:      +                        BETA, C( 1 ), NK )
434:                   CALL CGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ),
435:      +                        LDA, A( 1, NK+1 ), LDA, CBETA,
436:      +                        C( ( ( NK+1 )*NK )+1 ), NK )
437: *
438:                END IF
439: *
440:             ELSE
441: *
442: *              N is even, TRANSR = 'C', and UPLO = 'U'
443: *
444:                IF( NOTRANS ) THEN
445: *
446: *                 N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
447: *
448:                   CALL CHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
449:      +                        BETA, C( NK*( NK+1 )+1 ), NK )
450:                   CALL CHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
451:      +                        BETA, C( NK*NK+1 ), NK )
452:                   CALL CGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
453:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK )
454: *
455:                ELSE
456: *
457: *                 N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
458: *
459:                   CALL CHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
460:      +                        BETA, C( NK*( NK+1 )+1 ), NK )
461:                   CALL CHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
462:      +                        BETA, C( NK*NK+1 ), NK )
463:                   CALL CGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ),
464:      +                        LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK )
465: *
466:                END IF
467: *
468:             END IF
469: *
470:          END IF
471: *
472:       END IF
473: *
474:       RETURN
475: *
476: *     End of CHFRK
477: *
478:       END
479: