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