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+0, 0.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 MAX, DCMPLX
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.MAX( 1, 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: IF( MOD( 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( 1, 1 ), LDA,
220: $ BETA, C( 1 ), N )
221: CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
222: $ BETA, C( N+1 ), N )
223: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
224: $ LDA, A( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ),
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( 1, 1 ), 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( 1, 1 ),
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( 1, 1 ), LDA,
282: $ BETA, C( 1 ), N1 )
283: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
284: $ BETA, C( 2 ), N1 )
285: CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ),
286: $ LDA, A( N1+1, 1 ), 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( 1, 1 ), 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( 1, 1 ),
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( 1, 1 ), LDA,
312: $ BETA, C( N2*N2+1 ), N2 )
313: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
314: $ BETA, C( N1*N2+1 ), N2 )
315: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
316: $ LDA, A( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), LDA,
352: $ BETA, C( 2 ), N+1 )
353: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
354: $ BETA, C( 1 ), N+1 )
355: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
356: $ LDA, A( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), LDA,
382: $ BETA, C( NK+2 ), N+1 )
383: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
384: $ BETA, C( NK+1 ), N+1 )
385: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
386: $ LDA, A( NK+1, 1 ), 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( 1, 1 ), 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( 1, 1 ),
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( 1, 1 ), LDA,
418: $ BETA, C( NK+1 ), NK )
419: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
420: $ BETA, C( 1 ), NK )
421: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
422: $ LDA, A( NK+1, 1 ), 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( 1, 1 ), 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( 1, 1 ),
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( 1, 1 ), LDA,
448: $ BETA, C( NK*( NK+1 )+1 ), NK )
449: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
450: $ BETA, C( NK*NK+1 ), NK )
451: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
452: $ LDA, A( 1, 1 ), 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( 1, 1 ), 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( 1, 1 ), 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
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