![[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|
Return to zhfrk.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZHFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
2: + C )
3: *
4: * -- LAPACK routine (version 3.2.2) --
5: *
6: * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
7: * -- June 2010 --
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*conjg( A' ) + beta*C,
30: *
31: * or
32: *
33: * C := alpha*conjg( A' )*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
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
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
60: * On entry, TRANS specifies the operation to be performed as
61: * follows:
62: *
63: * TRANS = 'N' or 'n' C := alpha*A*conjg( A' ) + beta*C.
64: *
65: * TRANS = 'C' or 'c' C := alpha*conjg( A' )*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: * Arguments
111: * ==========
112: *
113: * ..
114: * .. Parameters ..
115: DOUBLE PRECISION ONE, ZERO
116: COMPLEX*16 CZERO
117: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
118: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
119: * ..
120: * .. Local Scalars ..
121: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
122: INTEGER INFO, NROWA, J, NK, N1, N2
123: COMPLEX*16 CALPHA, CBETA
124: * ..
125: * .. External Functions ..
126: LOGICAL LSAME
127: EXTERNAL LSAME
128: * ..
129: * .. External Subroutines ..
130: EXTERNAL XERBLA, ZGEMM, ZHERK
131: * ..
132: * .. Intrinsic Functions ..
133: INTRINSIC MAX, DCMPLX
134: * ..
135: * .. Executable Statements ..
136: *
137: *
138: * Test the input parameters.
139: *
140: INFO = 0
141: NORMALTRANSR = LSAME( TRANSR, 'N' )
142: LOWER = LSAME( UPLO, 'L' )
143: NOTRANS = LSAME( TRANS, 'N' )
144: *
145: IF( NOTRANS ) THEN
146: NROWA = N
147: ELSE
148: NROWA = K
149: END IF
150: *
151: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
152: INFO = -1
153: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
154: INFO = -2
155: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
156: INFO = -3
157: ELSE IF( N.LT.0 ) THEN
158: INFO = -4
159: ELSE IF( K.LT.0 ) THEN
160: INFO = -5
161: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
162: INFO = -8
163: END IF
164: IF( INFO.NE.0 ) THEN
165: CALL XERBLA( 'ZHFRK ', -INFO )
166: RETURN
167: END IF
168: *
169: * Quick return if possible.
170: *
171: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
172: * done (it is in ZHERK for example) and left in the general case.
173: *
174: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
175: + ( BETA.EQ.ONE ) ) )RETURN
176: *
177: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
178: DO J = 1, ( ( N*( N+1 ) ) / 2 )
179: C( J ) = CZERO
180: END DO
181: RETURN
182: END IF
183: *
184: CALPHA = DCMPLX( ALPHA, ZERO )
185: CBETA = DCMPLX( BETA, ZERO )
186: *
187: * C is N-by-N.
188: * If N is odd, set NISODD = .TRUE., and N1 and N2.
189: * If N is even, NISODD = .FALSE., and NK.
190: *
191: IF( MOD( N, 2 ).EQ.0 ) THEN
192: NISODD = .FALSE.
193: NK = N / 2
194: ELSE
195: NISODD = .TRUE.
196: IF( LOWER ) THEN
197: N2 = N / 2
198: N1 = N - N2
199: ELSE
200: N1 = N / 2
201: N2 = N - N1
202: END IF
203: END IF
204: *
205: IF( NISODD ) THEN
206: *
207: * N is odd
208: *
209: IF( NORMALTRANSR ) THEN
210: *
211: * N is odd and TRANSR = 'N'
212: *
213: IF( LOWER ) THEN
214: *
215: * N is odd, TRANSR = 'N', and UPLO = 'L'
216: *
217: IF( NOTRANS ) THEN
218: *
219: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
220: *
221: CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
222: + BETA, C( 1 ), N )
223: CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
224: + BETA, C( N+1 ), N )
225: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
226: + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N )
227: *
228: ELSE
229: *
230: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
231: *
232: CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
233: + BETA, C( 1 ), N )
234: CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
235: + BETA, C( N+1 ), N )
236: CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ),
237: + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N )
238: *
239: END IF
240: *
241: ELSE
242: *
243: * N is odd, TRANSR = 'N', and UPLO = 'U'
244: *
245: IF( NOTRANS ) THEN
246: *
247: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
248: *
249: CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
250: + BETA, C( N2+1 ), N )
251: CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
252: + BETA, C( N1+1 ), N )
253: CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ),
254: + LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N )
255: *
256: ELSE
257: *
258: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
259: *
260: CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
261: + BETA, C( N2+1 ), N )
262: CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N2 ), LDA,
263: + BETA, C( N1+1 ), N )
264: CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ),
265: + LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N )
266: *
267: END IF
268: *
269: END IF
270: *
271: ELSE
272: *
273: * N is odd, and TRANSR = 'C'
274: *
275: IF( LOWER ) THEN
276: *
277: * N is odd, TRANSR = 'C', and UPLO = 'L'
278: *
279: IF( NOTRANS ) THEN
280: *
281: * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
282: *
283: CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
284: + BETA, C( 1 ), N1 )
285: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
286: + BETA, C( 2 ), N1 )
287: CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ),
288: + LDA, A( N1+1, 1 ), LDA, CBETA,
289: + C( N1*N1+1 ), N1 )
290: *
291: ELSE
292: *
293: * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
294: *
295: CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
296: + BETA, C( 1 ), N1 )
297: CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
298: + BETA, C( 2 ), N1 )
299: CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ),
300: + LDA, A( 1, N1+1 ), LDA, CBETA,
301: + C( N1*N1+1 ), N1 )
302: *
303: END IF
304: *
305: ELSE
306: *
307: * N is odd, TRANSR = 'C', and UPLO = 'U'
308: *
309: IF( NOTRANS ) THEN
310: *
311: * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
312: *
313: CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
314: + BETA, C( N2*N2+1 ), N2 )
315: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
316: + BETA, C( N1*N2+1 ), N2 )
317: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ),
318: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 )
319: *
320: ELSE
321: *
322: * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
323: *
324: CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA,
325: + BETA, C( N2*N2+1 ), N2 )
326: CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA,
327: + BETA, C( N1*N2+1 ), N2 )
328: CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ),
329: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 )
330: *
331: END IF
332: *
333: END IF
334: *
335: END IF
336: *
337: ELSE
338: *
339: * N is even
340: *
341: IF( NORMALTRANSR ) THEN
342: *
343: * N is even and TRANSR = 'N'
344: *
345: IF( LOWER ) THEN
346: *
347: * N is even, TRANSR = 'N', and UPLO = 'L'
348: *
349: IF( NOTRANS ) THEN
350: *
351: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
352: *
353: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
354: + BETA, C( 2 ), N+1 )
355: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
356: + BETA, C( 1 ), N+1 )
357: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
358: + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ),
359: + N+1 )
360: *
361: ELSE
362: *
363: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
364: *
365: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
366: + BETA, C( 2 ), N+1 )
367: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
368: + BETA, C( 1 ), N+1 )
369: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ),
370: + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ),
371: + N+1 )
372: *
373: END IF
374: *
375: ELSE
376: *
377: * N is even, TRANSR = 'N', and UPLO = 'U'
378: *
379: IF( NOTRANS ) THEN
380: *
381: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
382: *
383: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
384: + BETA, C( NK+2 ), N+1 )
385: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
386: + BETA, C( NK+1 ), N+1 )
387: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
388: + LDA, A( NK+1, 1 ), LDA, CBETA, C( 1 ),
389: + N+1 )
390: *
391: ELSE
392: *
393: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
394: *
395: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
396: + BETA, C( NK+2 ), N+1 )
397: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
398: + BETA, C( NK+1 ), N+1 )
399: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ),
400: + LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ),
401: + N+1 )
402: *
403: END IF
404: *
405: END IF
406: *
407: ELSE
408: *
409: * N is even, and TRANSR = 'C'
410: *
411: IF( LOWER ) THEN
412: *
413: * N is even, TRANSR = 'C', and UPLO = 'L'
414: *
415: IF( NOTRANS ) THEN
416: *
417: * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
418: *
419: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
420: + BETA, C( NK+1 ), NK )
421: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
422: + BETA, C( 1 ), NK )
423: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ),
424: + LDA, A( NK+1, 1 ), LDA, CBETA,
425: + C( ( ( NK+1 )*NK )+1 ), NK )
426: *
427: ELSE
428: *
429: * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
430: *
431: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
432: + BETA, C( NK+1 ), NK )
433: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
434: + BETA, C( 1 ), NK )
435: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ),
436: + LDA, A( 1, NK+1 ), LDA, CBETA,
437: + C( ( ( NK+1 )*NK )+1 ), NK )
438: *
439: END IF
440: *
441: ELSE
442: *
443: * N is even, TRANSR = 'C', and UPLO = 'U'
444: *
445: IF( NOTRANS ) THEN
446: *
447: * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
448: *
449: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
450: + BETA, C( NK*( NK+1 )+1 ), NK )
451: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
452: + BETA, C( NK*NK+1 ), NK )
453: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ),
454: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK )
455: *
456: ELSE
457: *
458: * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
459: *
460: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA,
461: + BETA, C( NK*( NK+1 )+1 ), NK )
462: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA,
463: + BETA, C( NK*NK+1 ), NK )
464: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ),
465: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK )
466: *
467: END IF
468: *
469: END IF
470: *
471: END IF
472: *
473: END IF
474: *
475: RETURN
476: *
477: * End of ZHFRK
478: *
479: END