1 SUBROUTINE ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2 * .. Scalar Arguments ..
3 DOUBLE COMPLEX ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * ZHEMM performs one of the matrix-matrix operations
15 *
16 * C := alpha*A*B + beta*C,
17 *
18 * or
19 *
20 * C := alpha*B*A + beta*C,
21 *
22 * where alpha and beta are scalars, A is an hermitian matrix and B and
23 * C are m by n matrices.
24 *
25 * Arguments
26 * ==========
27 *
28 * SIDE - CHARACTER*1.
29 * On entry, SIDE specifies whether the hermitian matrix A
30 * appears on the left or right in the operation as follows:
31 *
32 * SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
33 *
34 * SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
35 *
36 * Unchanged on exit.
37 *
38 * UPLO - CHARACTER*1.
39 * On entry, UPLO specifies whether the upper or lower
40 * triangular part of the hermitian matrix A is to be
41 * referenced as follows:
42 *
43 * UPLO = 'U' or 'u' Only the upper triangular part of the
44 * hermitian matrix is to be referenced.
45 *
46 * UPLO = 'L' or 'l' Only the lower triangular part of the
47 * hermitian matrix is to be referenced.
48 *
49 * Unchanged on exit.
50 *
51 * M - INTEGER.
52 * On entry, M specifies the number of rows of the matrix C.
53 * M must be at least zero.
54 * Unchanged on exit.
55 *
56 * N - INTEGER.
57 * On entry, N specifies the number of columns of the matrix C.
58 * N must be at least zero.
59 * Unchanged on exit.
60 *
61 * ALPHA - COMPLEX*16 .
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
66 * m when SIDE = 'L' or 'l' and is n otherwise.
67 * Before entry with SIDE = 'L' or 'l', the m by m part of
68 * the array A must contain the hermitian matrix, such that
69 * when UPLO = 'U' or 'u', the leading m by m upper triangular
70 * part of the array A must contain the upper triangular part
71 * of the hermitian matrix and the strictly lower triangular
72 * part of A is not referenced, and when UPLO = 'L' or 'l',
73 * the leading m by m lower triangular part of the array A
74 * must contain the lower triangular part of the hermitian
75 * matrix and the strictly upper triangular part of A is not
76 * referenced.
77 * Before entry with SIDE = 'R' or 'r', the n by n part of
78 * the array A must contain the hermitian matrix, such that
79 * when UPLO = 'U' or 'u', the leading n by n upper triangular
80 * part of the array A must contain the upper triangular part
81 * of the hermitian matrix and the strictly lower triangular
82 * part of A is not referenced, and when UPLO = 'L' or 'l',
83 * the leading n by n lower triangular part of the array A
84 * must contain the lower triangular part of the hermitian
85 * matrix and the strictly upper triangular part of A is not
86 * referenced.
87 * Note that the imaginary parts of the diagonal elements need
88 * not be set, they are assumed to be zero.
89 * Unchanged on exit.
90 *
91 * LDA - INTEGER.
92 * On entry, LDA specifies the first dimension of A as declared
93 * in the calling (sub) program. When SIDE = 'L' or 'l' then
94 * LDA must be at least max( 1, m ), otherwise LDA must be at
95 * least max( 1, n ).
96 * Unchanged on exit.
97 *
98 * B - COMPLEX*16 array of DIMENSION ( LDB, n ).
99 * Before entry, the leading m by n part of the array B must
100 * contain the matrix B.
101 * Unchanged on exit.
102 *
103 * LDB - INTEGER.
104 * On entry, LDB specifies the first dimension of B as declared
105 * in the calling (sub) program. LDB must be at least
106 * max( 1, m ).
107 * Unchanged on exit.
108 *
109 * BETA - COMPLEX*16 .
110 * On entry, BETA specifies the scalar beta. When BETA is
111 * supplied as zero then C need not be set on input.
112 * Unchanged on exit.
113 *
114 * C - COMPLEX*16 array of DIMENSION ( LDC, n ).
115 * Before entry, the leading m by n part of the array C must
116 * contain the matrix C, except when beta is zero, in which
117 * case C need not be set on entry.
118 * On exit, the array C is overwritten by the m by n updated
119 * matrix.
120 *
121 * LDC - INTEGER.
122 * On entry, LDC specifies the first dimension of C as declared
123 * in the calling (sub) program. LDC must be at least
124 * max( 1, m ).
125 * Unchanged on exit.
126 *
127 * Further Details
128 * ===============
129 *
130 * Level 3 Blas routine.
131 *
132 * -- Written on 8-February-1989.
133 * Jack Dongarra, Argonne National Laboratory.
134 * Iain Duff, AERE Harwell.
135 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
136 * Sven Hammarling, Numerical Algorithms Group Ltd.
137 *
138 * =====================================================================
139 *
140 * .. External Functions ..
141 LOGICAL LSAME
142 EXTERNAL LSAME
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC DBLE,DCONJG,MAX
149 * ..
150 * .. Local Scalars ..
151 DOUBLE COMPLEX TEMP1,TEMP2
152 INTEGER I,INFO,J,K,NROWA
153 LOGICAL UPPER
154 * ..
155 * .. Parameters ..
156 DOUBLE COMPLEX ONE
157 PARAMETER (ONE= (1.0D+0,0.0D+0))
158 DOUBLE COMPLEX ZERO
159 PARAMETER (ZERO= (0.0D+0,0.0D+0))
160 * ..
161 *
162 * Set NROWA as the number of rows of A.
163 *
164 IF (LSAME(SIDE,'L')) THEN
165 NROWA = M
166 ELSE
167 NROWA = N
168 END IF
169 UPPER = LSAME(UPLO,'U')
170 *
171 * Test the input parameters.
172 *
173 INFO = 0
174 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
175 INFO = 1
176 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
177 INFO = 2
178 ELSE IF (M.LT.0) THEN
179 INFO = 3
180 ELSE IF (N.LT.0) THEN
181 INFO = 4
182 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
183 INFO = 7
184 ELSE IF (LDB.LT.MAX(1,M)) THEN
185 INFO = 9
186 ELSE IF (LDC.LT.MAX(1,M)) THEN
187 INFO = 12
188 END IF
189 IF (INFO.NE.0) THEN
190 CALL XERBLA('ZHEMM ',INFO)
191 RETURN
192 END IF
193 *
194 * Quick return if possible.
195 *
196 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
197 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
198 *
199 * And when alpha.eq.zero.
200 *
201 IF (ALPHA.EQ.ZERO) THEN
202 IF (BETA.EQ.ZERO) THEN
203 DO 20 J = 1,N
204 DO 10 I = 1,M
205 C(I,J) = ZERO
206 10 CONTINUE
207 20 CONTINUE
208 ELSE
209 DO 40 J = 1,N
210 DO 30 I = 1,M
211 C(I,J) = BETA*C(I,J)
212 30 CONTINUE
213 40 CONTINUE
214 END IF
215 RETURN
216 END IF
217 *
218 * Start the operations.
219 *
220 IF (LSAME(SIDE,'L')) THEN
221 *
222 * Form C := alpha*A*B + beta*C.
223 *
224 IF (UPPER) THEN
225 DO 70 J = 1,N
226 DO 60 I = 1,M
227 TEMP1 = ALPHA*B(I,J)
228 TEMP2 = ZERO
229 DO 50 K = 1,I - 1
230 C(K,J) = C(K,J) + TEMP1*A(K,I)
231 TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I))
232 50 CONTINUE
233 IF (BETA.EQ.ZERO) THEN
234 C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
235 ELSE
236 C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) +
237 + ALPHA*TEMP2
238 END IF
239 60 CONTINUE
240 70 CONTINUE
241 ELSE
242 DO 100 J = 1,N
243 DO 90 I = M,1,-1
244 TEMP1 = ALPHA*B(I,J)
245 TEMP2 = ZERO
246 DO 80 K = I + 1,M
247 C(K,J) = C(K,J) + TEMP1*A(K,I)
248 TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I))
249 80 CONTINUE
250 IF (BETA.EQ.ZERO) THEN
251 C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
252 ELSE
253 C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) +
254 + ALPHA*TEMP2
255 END IF
256 90 CONTINUE
257 100 CONTINUE
258 END IF
259 ELSE
260 *
261 * Form C := alpha*B*A + beta*C.
262 *
263 DO 170 J = 1,N
264 TEMP1 = ALPHA*DBLE(A(J,J))
265 IF (BETA.EQ.ZERO) THEN
266 DO 110 I = 1,M
267 C(I,J) = TEMP1*B(I,J)
268 110 CONTINUE
269 ELSE
270 DO 120 I = 1,M
271 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
272 120 CONTINUE
273 END IF
274 DO 140 K = 1,J - 1
275 IF (UPPER) THEN
276 TEMP1 = ALPHA*A(K,J)
277 ELSE
278 TEMP1 = ALPHA*DCONJG(A(J,K))
279 END IF
280 DO 130 I = 1,M
281 C(I,J) = C(I,J) + TEMP1*B(I,K)
282 130 CONTINUE
283 140 CONTINUE
284 DO 160 K = J + 1,N
285 IF (UPPER) THEN
286 TEMP1 = ALPHA*DCONJG(A(J,K))
287 ELSE
288 TEMP1 = ALPHA*A(K,J)
289 END IF
290 DO 150 I = 1,M
291 C(I,J) = C(I,J) + TEMP1*B(I,K)
292 150 CONTINUE
293 160 CONTINUE
294 170 CONTINUE
295 END IF
296 *
297 RETURN
298 *
299 * End of ZHEMM .
300 *
301 END
2 * .. Scalar Arguments ..
3 DOUBLE COMPLEX ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * ZHEMM performs one of the matrix-matrix operations
15 *
16 * C := alpha*A*B + beta*C,
17 *
18 * or
19 *
20 * C := alpha*B*A + beta*C,
21 *
22 * where alpha and beta are scalars, A is an hermitian matrix and B and
23 * C are m by n matrices.
24 *
25 * Arguments
26 * ==========
27 *
28 * SIDE - CHARACTER*1.
29 * On entry, SIDE specifies whether the hermitian matrix A
30 * appears on the left or right in the operation as follows:
31 *
32 * SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
33 *
34 * SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
35 *
36 * Unchanged on exit.
37 *
38 * UPLO - CHARACTER*1.
39 * On entry, UPLO specifies whether the upper or lower
40 * triangular part of the hermitian matrix A is to be
41 * referenced as follows:
42 *
43 * UPLO = 'U' or 'u' Only the upper triangular part of the
44 * hermitian matrix is to be referenced.
45 *
46 * UPLO = 'L' or 'l' Only the lower triangular part of the
47 * hermitian matrix is to be referenced.
48 *
49 * Unchanged on exit.
50 *
51 * M - INTEGER.
52 * On entry, M specifies the number of rows of the matrix C.
53 * M must be at least zero.
54 * Unchanged on exit.
55 *
56 * N - INTEGER.
57 * On entry, N specifies the number of columns of the matrix C.
58 * N must be at least zero.
59 * Unchanged on exit.
60 *
61 * ALPHA - COMPLEX*16 .
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
66 * m when SIDE = 'L' or 'l' and is n otherwise.
67 * Before entry with SIDE = 'L' or 'l', the m by m part of
68 * the array A must contain the hermitian matrix, such that
69 * when UPLO = 'U' or 'u', the leading m by m upper triangular
70 * part of the array A must contain the upper triangular part
71 * of the hermitian matrix and the strictly lower triangular
72 * part of A is not referenced, and when UPLO = 'L' or 'l',
73 * the leading m by m lower triangular part of the array A
74 * must contain the lower triangular part of the hermitian
75 * matrix and the strictly upper triangular part of A is not
76 * referenced.
77 * Before entry with SIDE = 'R' or 'r', the n by n part of
78 * the array A must contain the hermitian matrix, such that
79 * when UPLO = 'U' or 'u', the leading n by n upper triangular
80 * part of the array A must contain the upper triangular part
81 * of the hermitian matrix and the strictly lower triangular
82 * part of A is not referenced, and when UPLO = 'L' or 'l',
83 * the leading n by n lower triangular part of the array A
84 * must contain the lower triangular part of the hermitian
85 * matrix and the strictly upper triangular part of A is not
86 * referenced.
87 * Note that the imaginary parts of the diagonal elements need
88 * not be set, they are assumed to be zero.
89 * Unchanged on exit.
90 *
91 * LDA - INTEGER.
92 * On entry, LDA specifies the first dimension of A as declared
93 * in the calling (sub) program. When SIDE = 'L' or 'l' then
94 * LDA must be at least max( 1, m ), otherwise LDA must be at
95 * least max( 1, n ).
96 * Unchanged on exit.
97 *
98 * B - COMPLEX*16 array of DIMENSION ( LDB, n ).
99 * Before entry, the leading m by n part of the array B must
100 * contain the matrix B.
101 * Unchanged on exit.
102 *
103 * LDB - INTEGER.
104 * On entry, LDB specifies the first dimension of B as declared
105 * in the calling (sub) program. LDB must be at least
106 * max( 1, m ).
107 * Unchanged on exit.
108 *
109 * BETA - COMPLEX*16 .
110 * On entry, BETA specifies the scalar beta. When BETA is
111 * supplied as zero then C need not be set on input.
112 * Unchanged on exit.
113 *
114 * C - COMPLEX*16 array of DIMENSION ( LDC, n ).
115 * Before entry, the leading m by n part of the array C must
116 * contain the matrix C, except when beta is zero, in which
117 * case C need not be set on entry.
118 * On exit, the array C is overwritten by the m by n updated
119 * matrix.
120 *
121 * LDC - INTEGER.
122 * On entry, LDC specifies the first dimension of C as declared
123 * in the calling (sub) program. LDC must be at least
124 * max( 1, m ).
125 * Unchanged on exit.
126 *
127 * Further Details
128 * ===============
129 *
130 * Level 3 Blas routine.
131 *
132 * -- Written on 8-February-1989.
133 * Jack Dongarra, Argonne National Laboratory.
134 * Iain Duff, AERE Harwell.
135 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
136 * Sven Hammarling, Numerical Algorithms Group Ltd.
137 *
138 * =====================================================================
139 *
140 * .. External Functions ..
141 LOGICAL LSAME
142 EXTERNAL LSAME
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC DBLE,DCONJG,MAX
149 * ..
150 * .. Local Scalars ..
151 DOUBLE COMPLEX TEMP1,TEMP2
152 INTEGER I,INFO,J,K,NROWA
153 LOGICAL UPPER
154 * ..
155 * .. Parameters ..
156 DOUBLE COMPLEX ONE
157 PARAMETER (ONE= (1.0D+0,0.0D+0))
158 DOUBLE COMPLEX ZERO
159 PARAMETER (ZERO= (0.0D+0,0.0D+0))
160 * ..
161 *
162 * Set NROWA as the number of rows of A.
163 *
164 IF (LSAME(SIDE,'L')) THEN
165 NROWA = M
166 ELSE
167 NROWA = N
168 END IF
169 UPPER = LSAME(UPLO,'U')
170 *
171 * Test the input parameters.
172 *
173 INFO = 0
174 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
175 INFO = 1
176 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
177 INFO = 2
178 ELSE IF (M.LT.0) THEN
179 INFO = 3
180 ELSE IF (N.LT.0) THEN
181 INFO = 4
182 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
183 INFO = 7
184 ELSE IF (LDB.LT.MAX(1,M)) THEN
185 INFO = 9
186 ELSE IF (LDC.LT.MAX(1,M)) THEN
187 INFO = 12
188 END IF
189 IF (INFO.NE.0) THEN
190 CALL XERBLA('ZHEMM ',INFO)
191 RETURN
192 END IF
193 *
194 * Quick return if possible.
195 *
196 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
197 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
198 *
199 * And when alpha.eq.zero.
200 *
201 IF (ALPHA.EQ.ZERO) THEN
202 IF (BETA.EQ.ZERO) THEN
203 DO 20 J = 1,N
204 DO 10 I = 1,M
205 C(I,J) = ZERO
206 10 CONTINUE
207 20 CONTINUE
208 ELSE
209 DO 40 J = 1,N
210 DO 30 I = 1,M
211 C(I,J) = BETA*C(I,J)
212 30 CONTINUE
213 40 CONTINUE
214 END IF
215 RETURN
216 END IF
217 *
218 * Start the operations.
219 *
220 IF (LSAME(SIDE,'L')) THEN
221 *
222 * Form C := alpha*A*B + beta*C.
223 *
224 IF (UPPER) THEN
225 DO 70 J = 1,N
226 DO 60 I = 1,M
227 TEMP1 = ALPHA*B(I,J)
228 TEMP2 = ZERO
229 DO 50 K = 1,I - 1
230 C(K,J) = C(K,J) + TEMP1*A(K,I)
231 TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I))
232 50 CONTINUE
233 IF (BETA.EQ.ZERO) THEN
234 C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
235 ELSE
236 C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) +
237 + ALPHA*TEMP2
238 END IF
239 60 CONTINUE
240 70 CONTINUE
241 ELSE
242 DO 100 J = 1,N
243 DO 90 I = M,1,-1
244 TEMP1 = ALPHA*B(I,J)
245 TEMP2 = ZERO
246 DO 80 K = I + 1,M
247 C(K,J) = C(K,J) + TEMP1*A(K,I)
248 TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I))
249 80 CONTINUE
250 IF (BETA.EQ.ZERO) THEN
251 C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
252 ELSE
253 C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) +
254 + ALPHA*TEMP2
255 END IF
256 90 CONTINUE
257 100 CONTINUE
258 END IF
259 ELSE
260 *
261 * Form C := alpha*B*A + beta*C.
262 *
263 DO 170 J = 1,N
264 TEMP1 = ALPHA*DBLE(A(J,J))
265 IF (BETA.EQ.ZERO) THEN
266 DO 110 I = 1,M
267 C(I,J) = TEMP1*B(I,J)
268 110 CONTINUE
269 ELSE
270 DO 120 I = 1,M
271 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
272 120 CONTINUE
273 END IF
274 DO 140 K = 1,J - 1
275 IF (UPPER) THEN
276 TEMP1 = ALPHA*A(K,J)
277 ELSE
278 TEMP1 = ALPHA*DCONJG(A(J,K))
279 END IF
280 DO 130 I = 1,M
281 C(I,J) = C(I,J) + TEMP1*B(I,K)
282 130 CONTINUE
283 140 CONTINUE
284 DO 160 K = J + 1,N
285 IF (UPPER) THEN
286 TEMP1 = ALPHA*DCONJG(A(J,K))
287 ELSE
288 TEMP1 = ALPHA*A(K,J)
289 END IF
290 DO 150 I = 1,M
291 C(I,J) = C(I,J) + TEMP1*B(I,K)
292 150 CONTINUE
293 160 CONTINUE
294 170 CONTINUE
295 END IF
296 *
297 RETURN
298 *
299 * End of ZHEMM .
300 *
301 END