LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cher2k.f
Go to the documentation of this file.
1 *> \brief \b CHER2K
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
12 *
13 * .. Scalar Arguments ..
14 * COMPLEX ALPHA
15 * REAL BETA
16 * INTEGER K,LDA,LDB,LDC,N
17 * CHARACTER TRANS,UPLO
18 * ..
19 * .. Array Arguments ..
20 * COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> CHER2K performs one of the hermitian rank 2k operations
30 *>
31 *> C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C,
32 *>
33 *> or
34 *>
35 *> C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C,
36 *>
37 *> where alpha and beta are scalars with beta real, C is an n by n
38 *> hermitian matrix and A and B are n by k matrices in the first case
39 *> and k by n matrices in the second case.
40 *> \endverbatim
41 *
42 * Arguments:
43 * ==========
44 *
45 *> \param[in] UPLO
46 *> \verbatim
47 *> UPLO is CHARACTER*1
48 *> On entry, UPLO specifies whether the upper or lower
49 *> triangular part of the array C is to be referenced as
50 *> follows:
51 *>
52 *> UPLO = 'U' or 'u' Only the upper triangular part of C
53 *> is to be referenced.
54 *>
55 *> UPLO = 'L' or 'l' Only the lower triangular part of C
56 *> is to be referenced.
57 *> \endverbatim
58 *>
59 *> \param[in] TRANS
60 *> \verbatim
61 *> TRANS is CHARACTER*1
62 *> On entry, TRANS specifies the operation to be performed as
63 *> follows:
64 *>
65 *> TRANS = 'N' or 'n' C := alpha*A*B**H +
66 *> conjg( alpha )*B*A**H +
67 *> beta*C.
68 *>
69 *> TRANS = 'C' or 'c' C := alpha*A**H*B +
70 *> conjg( alpha )*B**H*A +
71 *> beta*C.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> On entry, N specifies the order of the matrix C. N must be
78 *> at least zero.
79 *> \endverbatim
80 *>
81 *> \param[in] K
82 *> \verbatim
83 *> K is INTEGER
84 *> On entry with TRANS = 'N' or 'n', K specifies the number
85 *> of columns of the matrices A and B, and on entry with
86 *> TRANS = 'C' or 'c', K specifies the number of rows of the
87 *> matrices A and B. K must be at least zero.
88 *> \endverbatim
89 *>
90 *> \param[in] ALPHA
91 *> \verbatim
92 *> ALPHA is COMPLEX
93 *> On entry, ALPHA specifies the scalar alpha.
94 *> \endverbatim
95 *>
96 *> \param[in] A
97 *> \verbatim
98 *> A is COMPLEX array of DIMENSION ( LDA, ka ), where ka is
99 *> k when TRANS = 'N' or 'n', and is n otherwise.
100 *> Before entry with TRANS = 'N' or 'n', the leading n by k
101 *> part of the array A must contain the matrix A, otherwise
102 *> the leading k by n part of the array A must contain the
103 *> matrix A.
104 *> \endverbatim
105 *>
106 *> \param[in] LDA
107 *> \verbatim
108 *> LDA is INTEGER
109 *> On entry, LDA specifies the first dimension of A as declared
110 *> in the calling (sub) program. When TRANS = 'N' or 'n'
111 *> then LDA must be at least max( 1, n ), otherwise LDA must
112 *> be at least max( 1, k ).
113 *> \endverbatim
114 *>
115 *> \param[in] B
116 *> \verbatim
117 *> B is COMPLEX array of DIMENSION ( LDB, kb ), where kb is
118 *> k when TRANS = 'N' or 'n', and is n otherwise.
119 *> Before entry with TRANS = 'N' or 'n', the leading n by k
120 *> part of the array B must contain the matrix B, otherwise
121 *> the leading k by n part of the array B must contain the
122 *> matrix B.
123 *> \endverbatim
124 *>
125 *> \param[in] LDB
126 *> \verbatim
127 *> LDB is INTEGER
128 *> On entry, LDB specifies the first dimension of B as declared
129 *> in the calling (sub) program. When TRANS = 'N' or 'n'
130 *> then LDB must be at least max( 1, n ), otherwise LDB must
131 *> be at least max( 1, k ).
132 *> \endverbatim
133 *>
134 *> \param[in] BETA
135 *> \verbatim
136 *> BETA is REAL
137 *> On entry, BETA specifies the scalar beta.
138 *> \endverbatim
139 *>
140 *> \param[in,out] C
141 *> \verbatim
142 *> C is COMPLEX array of DIMENSION ( LDC, n ).
143 *> Before entry with UPLO = 'U' or 'u', the leading n by n
144 *> upper triangular part of the array C must contain the upper
145 *> triangular part of the hermitian matrix and the strictly
146 *> lower triangular part of C is not referenced. On exit, the
147 *> upper triangular part of the array C is overwritten by the
148 *> upper triangular part of the updated matrix.
149 *> Before entry with UPLO = 'L' or 'l', the leading n by n
150 *> lower triangular part of the array C must contain the lower
151 *> triangular part of the hermitian matrix and the strictly
152 *> upper triangular part of C is not referenced. On exit, the
153 *> lower triangular part of the array C is overwritten by the
154 *> lower triangular part of the updated matrix.
155 *> Note that the imaginary parts of the diagonal elements need
156 *> not be set, they are assumed to be zero, and on exit they
157 *> are set to zero.
158 *> \endverbatim
159 *>
160 *> \param[in] LDC
161 *> \verbatim
162 *> LDC is INTEGER
163 *> On entry, LDC specifies the first dimension of C as declared
164 *> in the calling (sub) program. LDC must be at least
165 *> max( 1, n ).
166 *> \endverbatim
167 *
168 * Authors:
169 * ========
170 *
171 *> \author Univ. of Tennessee
172 *> \author Univ. of California Berkeley
173 *> \author Univ. of Colorado Denver
174 *> \author NAG Ltd.
175 *
176 *> \date November 2011
177 *
178 *> \ingroup complex_blas_level3
179 *
180 *> \par Further Details:
181 * =====================
182 *>
183 *> \verbatim
184 *>
185 *> Level 3 Blas routine.
186 *>
187 *> -- Written on 8-February-1989.
188 *> Jack Dongarra, Argonne National Laboratory.
189 *> Iain Duff, AERE Harwell.
190 *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
191 *> Sven Hammarling, Numerical Algorithms Group Ltd.
192 *>
193 *> -- Modified 8-Nov-93 to set C(J,J) to REAL( C(J,J) ) when BETA = 1.
194 *> Ed Anderson, Cray Research Inc.
195 *> \endverbatim
196 *>
197 * =====================================================================
198  SUBROUTINE cher2k(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
199 *
200 * -- Reference BLAS level3 routine (version 3.4.0) --
201 * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
202 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203 * November 2011
204 *
205 * .. Scalar Arguments ..
206  COMPLEX alpha
207  REAL beta
208  INTEGER k,lda,ldb,ldc,n
209  CHARACTER trans,uplo
210 * ..
211 * .. Array Arguments ..
212  COMPLEX a(lda,*),b(ldb,*),c(ldc,*)
213 * ..
214 *
215 * =====================================================================
216 *
217 * .. External Functions ..
218  LOGICAL lsame
219  EXTERNAL lsame
220 * ..
221 * .. External Subroutines ..
222  EXTERNAL xerbla
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC conjg,max,real
226 * ..
227 * .. Local Scalars ..
228  COMPLEX temp1,temp2
229  INTEGER i,info,j,l,nrowa
230  LOGICAL upper
231 * ..
232 * .. Parameters ..
233  REAL one
234  parameter(one=1.0e+0)
235  COMPLEX zero
236  parameter(zero= (0.0e+0,0.0e+0))
237 * ..
238 *
239 * Test the input parameters.
240 *
241  IF (lsame(trans,'N')) THEN
242  nrowa = n
243  ELSE
244  nrowa = k
245  END IF
246  upper = lsame(uplo,'U')
247 *
248  info = 0
249  IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
250  info = 1
251  ELSE IF ((.NOT.lsame(trans,'N')) .AND.
252  + (.NOT.lsame(trans,'C'))) THEN
253  info = 2
254  ELSE IF (n.LT.0) THEN
255  info = 3
256  ELSE IF (k.LT.0) THEN
257  info = 4
258  ELSE IF (lda.LT.max(1,nrowa)) THEN
259  info = 7
260  ELSE IF (ldb.LT.max(1,nrowa)) THEN
261  info = 9
262  ELSE IF (ldc.LT.max(1,n)) THEN
263  info = 12
264  END IF
265  IF (info.NE.0) THEN
266  CALL xerbla('CHER2K',info)
267  return
268  END IF
269 *
270 * Quick return if possible.
271 *
272  IF ((n.EQ.0) .OR. (((alpha.EQ.zero).OR.
273  + (k.EQ.0)).AND. (beta.EQ.one))) return
274 *
275 * And when alpha.eq.zero.
276 *
277  IF (alpha.EQ.zero) THEN
278  IF (upper) THEN
279  IF (beta.EQ.REAL(zero)) then
280  DO 20 j = 1,n
281  DO 10 i = 1,j
282  c(i,j) = zero
283  10 continue
284  20 continue
285  ELSE
286  DO 40 j = 1,n
287  DO 30 i = 1,j - 1
288  c(i,j) = beta*c(i,j)
289  30 continue
290  c(j,j) = beta*REAL(c(j,j))
291  40 continue
292  END IF
293  ELSE
294  IF (beta.EQ.REAL(zero)) then
295  DO 60 j = 1,n
296  DO 50 i = j,n
297  c(i,j) = zero
298  50 continue
299  60 continue
300  ELSE
301  DO 80 j = 1,n
302  c(j,j) = beta*REAL(c(j,j))
303  DO 70 i = j + 1,n
304  c(i,j) = beta*c(i,j)
305  70 continue
306  80 continue
307  END IF
308  END IF
309  return
310  END IF
311 *
312 * Start the operations.
313 *
314  IF (lsame(trans,'N')) THEN
315 *
316 * Form C := alpha*A*B**H + conjg( alpha )*B*A**H +
317 * C.
318 *
319  IF (upper) THEN
320  DO 130 j = 1,n
321  IF (beta.EQ.REAL(zero)) then
322  DO 90 i = 1,j
323  c(i,j) = zero
324  90 continue
325  ELSE IF (beta.NE.one) THEN
326  DO 100 i = 1,j - 1
327  c(i,j) = beta*c(i,j)
328  100 continue
329  c(j,j) = beta*REAL(c(j,j))
330  ELSE
331  c(j,j) = REAL(c(j,j))
332  END IF
333  DO 120 l = 1,k
334  IF ((a(j,l).NE.zero) .OR. (b(j,l).NE.zero)) THEN
335  temp1 = alpha*conjg(b(j,l))
336  temp2 = conjg(alpha*a(j,l))
337  DO 110 i = 1,j - 1
338  c(i,j) = c(i,j) + a(i,l)*temp1 +
339  + b(i,l)*temp2
340  110 continue
341  c(j,j) = REAL(C(J,J)) +
342  + REAL(a(j,l)*temp1+b(j,l)*temp2)
343  END IF
344  120 continue
345  130 continue
346  ELSE
347  DO 180 j = 1,n
348  IF (beta.EQ.REAL(zero)) then
349  DO 140 i = j,n
350  c(i,j) = zero
351  140 continue
352  ELSE IF (beta.NE.one) THEN
353  DO 150 i = j + 1,n
354  c(i,j) = beta*c(i,j)
355  150 continue
356  c(j,j) = beta*REAL(c(j,j))
357  ELSE
358  c(j,j) = REAL(c(j,j))
359  END IF
360  DO 170 l = 1,k
361  IF ((a(j,l).NE.zero) .OR. (b(j,l).NE.zero)) THEN
362  temp1 = alpha*conjg(b(j,l))
363  temp2 = conjg(alpha*a(j,l))
364  DO 160 i = j + 1,n
365  c(i,j) = c(i,j) + a(i,l)*temp1 +
366  + b(i,l)*temp2
367  160 continue
368  c(j,j) = REAL(C(J,J)) +
369  + REAL(a(j,l)*temp1+b(j,l)*temp2)
370  END IF
371  170 continue
372  180 continue
373  END IF
374  ELSE
375 *
376 * Form C := alpha*A**H*B + conjg( alpha )*B**H*A +
377 * C.
378 *
379  IF (upper) THEN
380  DO 210 j = 1,n
381  DO 200 i = 1,j
382  temp1 = zero
383  temp2 = zero
384  DO 190 l = 1,k
385  temp1 = temp1 + conjg(a(l,i))*b(l,j)
386  temp2 = temp2 + conjg(b(l,i))*a(l,j)
387  190 continue
388  IF (i.EQ.j) THEN
389  IF (beta.EQ.REAL(zero)) then
390  c(j,j) = REAL(alpha*temp1+
391  + conjg(alpha)*temp2)
392  ELSE
393  c(j,j) = beta*REAL(C(J,J)) +
394  + REAL(alpha*temp1+
395  + conjg(alpha)*temp2)
396  END IF
397  ELSE
398  IF (beta.EQ.REAL(zero)) then
399  c(i,j) = alpha*temp1 + conjg(alpha)*temp2
400  ELSE
401  c(i,j) = beta*c(i,j) + alpha*temp1 +
402  + conjg(alpha)*temp2
403  END IF
404  END IF
405  200 continue
406  210 continue
407  ELSE
408  DO 240 j = 1,n
409  DO 230 i = j,n
410  temp1 = zero
411  temp2 = zero
412  DO 220 l = 1,k
413  temp1 = temp1 + conjg(a(l,i))*b(l,j)
414  temp2 = temp2 + conjg(b(l,i))*a(l,j)
415  220 continue
416  IF (i.EQ.j) THEN
417  IF (beta.EQ.REAL(zero)) then
418  c(j,j) = REAL(alpha*temp1+
419  + conjg(alpha)*temp2)
420  ELSE
421  c(j,j) = beta*REAL(C(J,J)) +
422  + REAL(alpha*temp1+
423  + conjg(alpha)*temp2)
424  END IF
425  ELSE
426  IF (beta.EQ.REAL(zero)) then
427  c(i,j) = alpha*temp1 + conjg(alpha)*temp2
428  ELSE
429  c(i,j) = beta*c(i,j) + alpha*temp1 +
430  + conjg(alpha)*temp2
431  END IF
432  END IF
433  230 continue
434  240 continue
435  END IF
436  END IF
437 *
438  return
439 *
440 * End of CHER2K.
441 *
442  END