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