LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zhfrk.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:173
subroutine zhfrk(TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C)
ZHFRK performs a Hermitian rank-k operation for matrix in RFP format.
Definition: zhfrk.f:168