LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 of 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 *> \date September 2012
164 *
165 *> \ingroup complex16OTHERcomputational
166 *
167 * =====================================================================
168  SUBROUTINE zhfrk( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
169  $ c )
170 *
171 * -- LAPACK computational routine (version 3.4.2) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * September 2012
175 *
176 * .. Scalar Arguments ..
177  DOUBLE PRECISION alpha, beta
178  INTEGER k, lda, n
179  CHARACTER trans, transr, uplo
180 * ..
181 * .. Array Arguments ..
182  COMPLEX*16 a( lda, * ), c( * )
183 * ..
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188  DOUBLE PRECISION one, zero
189  COMPLEX*16 czero
190  parameter( one = 1.0d+0, zero = 0.0d+0 )
191  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
192 * ..
193 * .. Local Scalars ..
194  LOGICAL lower, normaltransr, nisodd, notrans
195  INTEGER info, nrowa, j, nk, n1, n2
196  COMPLEX*16 calpha, cbeta
197 * ..
198 * .. External Functions ..
199  LOGICAL lsame
200  EXTERNAL lsame
201 * ..
202 * .. External Subroutines ..
203  EXTERNAL xerbla, zgemm, zherk
204 * ..
205 * .. Intrinsic Functions ..
206  INTRINSIC max, dcmplx
207 * ..
208 * .. Executable Statements ..
209 *
210 *
211 * Test the input parameters.
212 *
213  info = 0
214  normaltransr = lsame( transr, 'N' )
215  lower = lsame( uplo, 'L' )
216  notrans = lsame( trans, 'N' )
217 *
218  IF( notrans ) THEN
219  nrowa = n
220  ELSE
221  nrowa = k
222  END IF
223 *
224  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
225  info = -1
226  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
227  info = -2
228  ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'C' ) ) THEN
229  info = -3
230  ELSE IF( n.LT.0 ) THEN
231  info = -4
232  ELSE IF( k.LT.0 ) THEN
233  info = -5
234  ELSE IF( lda.LT.max( 1, nrowa ) ) THEN
235  info = -8
236  END IF
237  IF( info.NE.0 ) THEN
238  CALL xerbla( 'ZHFRK ', -info )
239  return
240  END IF
241 *
242 * Quick return if possible.
243 *
244 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
245 * done (it is in ZHERK for example) and left in the general case.
246 *
247  IF( ( n.EQ.0 ) .OR. ( ( ( alpha.EQ.zero ) .OR. ( k.EQ.0 ) ) .AND.
248  $ ( beta.EQ.one ) ) )return
249 *
250  IF( ( alpha.EQ.zero ) .AND. ( beta.EQ.zero ) ) THEN
251  DO j = 1, ( ( n*( n+1 ) ) / 2 )
252  c( j ) = czero
253  END DO
254  return
255  END IF
256 *
257  calpha = dcmplx( alpha, zero )
258  cbeta = dcmplx( beta, zero )
259 *
260 * C is N-by-N.
261 * If N is odd, set NISODD = .TRUE., and N1 and N2.
262 * If N is even, NISODD = .FALSE., and NK.
263 *
264  IF( mod( n, 2 ).EQ.0 ) THEN
265  nisodd = .false.
266  nk = n / 2
267  ELSE
268  nisodd = .true.
269  IF( lower ) THEN
270  n2 = n / 2
271  n1 = n - n2
272  ELSE
273  n1 = n / 2
274  n2 = n - n1
275  END IF
276  END IF
277 *
278  IF( nisodd ) THEN
279 *
280 * N is odd
281 *
282  IF( normaltransr ) THEN
283 *
284 * N is odd and TRANSR = 'N'
285 *
286  IF( lower ) THEN
287 *
288 * N is odd, TRANSR = 'N', and UPLO = 'L'
289 *
290  IF( notrans ) THEN
291 *
292 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
293 *
294  CALL zherk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
295  $ beta, c( 1 ), n )
296  CALL zherk( 'U', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
297  $ beta, c( n+1 ), n )
298  CALL zgemm( 'N', 'C', n2, n1, k, calpha, a( n1+1, 1 ),
299  $ lda, a( 1, 1 ), lda, cbeta, c( n1+1 ), n )
300 *
301  ELSE
302 *
303 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
304 *
305  CALL zherk( 'L', 'C', n1, k, alpha, a( 1, 1 ), lda,
306  $ beta, c( 1 ), n )
307  CALL zherk( 'U', 'C', n2, k, alpha, a( 1, n1+1 ), lda,
308  $ beta, c( n+1 ), n )
309  CALL zgemm( 'C', 'N', n2, n1, k, calpha, a( 1, n1+1 ),
310  $ lda, a( 1, 1 ), lda, cbeta, c( n1+1 ), n )
311 *
312  END IF
313 *
314  ELSE
315 *
316 * N is odd, TRANSR = 'N', and UPLO = 'U'
317 *
318  IF( notrans ) THEN
319 *
320 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
321 *
322  CALL zherk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
323  $ beta, c( n2+1 ), n )
324  CALL zherk( 'U', 'N', n2, k, alpha, a( n2, 1 ), lda,
325  $ beta, c( n1+1 ), n )
326  CALL zgemm( 'N', 'C', n1, n2, k, calpha, a( 1, 1 ),
327  $ lda, a( n2, 1 ), lda, cbeta, c( 1 ), n )
328 *
329  ELSE
330 *
331 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
332 *
333  CALL zherk( 'L', 'C', n1, k, alpha, a( 1, 1 ), lda,
334  $ beta, c( n2+1 ), n )
335  CALL zherk( 'U', 'C', n2, k, alpha, a( 1, n2 ), lda,
336  $ beta, c( n1+1 ), n )
337  CALL zgemm( 'C', 'N', n1, n2, k, calpha, a( 1, 1 ),
338  $ lda, a( 1, n2 ), lda, cbeta, c( 1 ), n )
339 *
340  END IF
341 *
342  END IF
343 *
344  ELSE
345 *
346 * N is odd, and TRANSR = 'C'
347 *
348  IF( lower ) THEN
349 *
350 * N is odd, TRANSR = 'C', and UPLO = 'L'
351 *
352  IF( notrans ) THEN
353 *
354 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
355 *
356  CALL zherk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
357  $ beta, c( 1 ), n1 )
358  CALL zherk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
359  $ beta, c( 2 ), n1 )
360  CALL zgemm( 'N', 'C', n1, n2, k, calpha, a( 1, 1 ),
361  $ lda, a( n1+1, 1 ), lda, cbeta,
362  $ c( n1*n1+1 ), n1 )
363 *
364  ELSE
365 *
366 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
367 *
368  CALL zherk( 'U', 'C', n1, k, alpha, a( 1, 1 ), lda,
369  $ beta, c( 1 ), n1 )
370  CALL zherk( 'L', 'C', n2, k, alpha, a( 1, n1+1 ), lda,
371  $ beta, c( 2 ), n1 )
372  CALL zgemm( 'C', 'N', n1, n2, k, calpha, a( 1, 1 ),
373  $ lda, a( 1, n1+1 ), lda, cbeta,
374  $ c( n1*n1+1 ), n1 )
375 *
376  END IF
377 *
378  ELSE
379 *
380 * N is odd, TRANSR = 'C', and UPLO = 'U'
381 *
382  IF( notrans ) THEN
383 *
384 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
385 *
386  CALL zherk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
387  $ beta, c( n2*n2+1 ), n2 )
388  CALL zherk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
389  $ beta, c( n1*n2+1 ), n2 )
390  CALL zgemm( 'N', 'C', n2, n1, k, calpha, a( n1+1, 1 ),
391  $ lda, a( 1, 1 ), lda, cbeta, c( 1 ), n2 )
392 *
393  ELSE
394 *
395 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
396 *
397  CALL zherk( 'U', 'C', n1, k, alpha, a( 1, 1 ), lda,
398  $ beta, c( n2*n2+1 ), n2 )
399  CALL zherk( 'L', 'C', n2, k, alpha, a( 1, n1+1 ), lda,
400  $ beta, c( n1*n2+1 ), n2 )
401  CALL zgemm( 'C', 'N', n2, n1, k, calpha, a( 1, n1+1 ),
402  $ lda, a( 1, 1 ), lda, cbeta, c( 1 ), n2 )
403 *
404  END IF
405 *
406  END IF
407 *
408  END IF
409 *
410  ELSE
411 *
412 * N is even
413 *
414  IF( normaltransr ) THEN
415 *
416 * N is even and TRANSR = 'N'
417 *
418  IF( lower ) THEN
419 *
420 * N is even, TRANSR = 'N', and UPLO = 'L'
421 *
422  IF( notrans ) THEN
423 *
424 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
425 *
426  CALL zherk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
427  $ beta, c( 2 ), n+1 )
428  CALL zherk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
429  $ beta, c( 1 ), n+1 )
430  CALL zgemm( 'N', 'C', nk, nk, k, calpha, a( nk+1, 1 ),
431  $ lda, a( 1, 1 ), lda, cbeta, c( nk+2 ),
432  $ n+1 )
433 *
434  ELSE
435 *
436 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C'
437 *
438  CALL zherk( 'L', 'C', nk, k, alpha, a( 1, 1 ), lda,
439  $ beta, c( 2 ), n+1 )
440  CALL zherk( 'U', 'C', nk, k, alpha, a( 1, nk+1 ), lda,
441  $ beta, c( 1 ), n+1 )
442  CALL zgemm( 'C', 'N', nk, nk, k, calpha, a( 1, nk+1 ),
443  $ lda, a( 1, 1 ), lda, cbeta, c( nk+2 ),
444  $ n+1 )
445 *
446  END IF
447 *
448  ELSE
449 *
450 * N is even, TRANSR = 'N', and UPLO = 'U'
451 *
452  IF( notrans ) THEN
453 *
454 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
455 *
456  CALL zherk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
457  $ beta, c( nk+2 ), n+1 )
458  CALL zherk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
459  $ beta, c( nk+1 ), n+1 )
460  CALL zgemm( 'N', 'C', nk, nk, k, calpha, a( 1, 1 ),
461  $ lda, a( nk+1, 1 ), lda, cbeta, c( 1 ),
462  $ n+1 )
463 *
464  ELSE
465 *
466 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C'
467 *
468  CALL zherk( 'L', 'C', nk, k, alpha, a( 1, 1 ), lda,
469  $ beta, c( nk+2 ), n+1 )
470  CALL zherk( 'U', 'C', nk, k, alpha, a( 1, nk+1 ), lda,
471  $ beta, c( nk+1 ), n+1 )
472  CALL zgemm( 'C', 'N', nk, nk, k, calpha, a( 1, 1 ),
473  $ lda, a( 1, nk+1 ), lda, cbeta, c( 1 ),
474  $ n+1 )
475 *
476  END IF
477 *
478  END IF
479 *
480  ELSE
481 *
482 * N is even, and TRANSR = 'C'
483 *
484  IF( lower ) THEN
485 *
486 * N is even, TRANSR = 'C', and UPLO = 'L'
487 *
488  IF( notrans ) THEN
489 *
490 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N'
491 *
492  CALL zherk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
493  $ beta, c( nk+1 ), nk )
494  CALL zherk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
495  $ beta, c( 1 ), nk )
496  CALL zgemm( 'N', 'C', nk, nk, k, calpha, a( 1, 1 ),
497  $ lda, a( nk+1, 1 ), lda, cbeta,
498  $ c( ( ( nk+1 )*nk )+1 ), nk )
499 *
500  ELSE
501 *
502 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C'
503 *
504  CALL zherk( 'U', 'C', nk, k, alpha, a( 1, 1 ), lda,
505  $ beta, c( nk+1 ), nk )
506  CALL zherk( 'L', 'C', nk, k, alpha, a( 1, nk+1 ), lda,
507  $ beta, c( 1 ), nk )
508  CALL zgemm( 'C', 'N', nk, nk, k, calpha, a( 1, 1 ),
509  $ lda, a( 1, nk+1 ), lda, cbeta,
510  $ c( ( ( nk+1 )*nk )+1 ), nk )
511 *
512  END IF
513 *
514  ELSE
515 *
516 * N is even, TRANSR = 'C', and UPLO = 'U'
517 *
518  IF( notrans ) THEN
519 *
520 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N'
521 *
522  CALL zherk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
523  $ beta, c( nk*( nk+1 )+1 ), nk )
524  CALL zherk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
525  $ beta, c( nk*nk+1 ), nk )
526  CALL zgemm( 'N', 'C', nk, nk, k, calpha, a( nk+1, 1 ),
527  $ lda, a( 1, 1 ), lda, cbeta, c( 1 ), nk )
528 *
529  ELSE
530 *
531 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C'
532 *
533  CALL zherk( 'U', 'C', nk, k, alpha, a( 1, 1 ), lda,
534  $ beta, c( nk*( nk+1 )+1 ), nk )
535  CALL zherk( 'L', 'C', nk, k, alpha, a( 1, nk+1 ), lda,
536  $ beta, c( nk*nk+1 ), nk )
537  CALL zgemm( 'C', 'N', nk, nk, k, calpha, a( 1, nk+1 ),
538  $ lda, a( 1, 1 ), lda, cbeta, c( 1 ), nk )
539 *
540  END IF
541 *
542  END IF
543 *
544  END IF
545 *
546  END IF
547 *
548  return
549 *
550 * End of ZHFRK
551 *
552  END