![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsyr2k.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas|
Return to dsyr2k.f CVS log | Up to [local] / rpl / lapack / blas |
1.8 bertrand 1: *> \brief \b DSYR2K
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.13 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: * Definition:
9: * ===========
10: *
11: * SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
1.13 bertrand 12: *
1.8 bertrand 13: * .. Scalar Arguments ..
14: * DOUBLE PRECISION ALPHA,BETA
15: * INTEGER K,LDA,LDB,LDC,N
16: * CHARACTER TRANS,UPLO
17: * ..
18: * .. Array Arguments ..
19: * DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
20: * ..
1.13 bertrand 21: *
1.8 bertrand 22: *
23: *> \par Purpose:
24: * =============
25: *>
26: *> \verbatim
27: *>
28: *> DSYR2K performs one of the symmetric rank 2k operations
29: *>
30: *> C := alpha*A*B**T + alpha*B*A**T + beta*C,
31: *>
32: *> or
33: *>
34: *> C := alpha*A**T*B + alpha*B**T*A + beta*C,
35: *>
36: *> where alpha and beta are scalars, C is an n by n symmetric matrix
37: *> and A and B are n by k matrices in the first case and k by n
38: *> matrices 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*B**T + alpha*B*A**T +
65: *> beta*C.
66: *>
67: *> TRANS = 'T' or 't' C := alpha*A**T*B + alpha*B**T*A +
68: *> beta*C.
69: *>
70: *> TRANS = 'C' or 'c' C := alpha*A**T*B + alpha*B**T*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 = 'T' or 't' or 'C' or 'c', K specifies the number
87: *> of rows of the matrices A and B. K must be at least zero.
88: *> \endverbatim
89: *>
90: *> \param[in] ALPHA
91: *> \verbatim
92: *> ALPHA is DOUBLE PRECISION.
93: *> On entry, ALPHA specifies the scalar alpha.
94: *> \endverbatim
95: *>
96: *> \param[in] A
97: *> \verbatim
1.14 bertrand 98: *> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
1.8 bertrand 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
1.14 bertrand 117: *> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
1.8 bertrand 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 DOUBLE PRECISION.
137: *> On entry, BETA specifies the scalar beta.
138: *> \endverbatim
139: *>
140: *> \param[in,out] C
141: *> \verbatim
1.14 bertrand 142: *> C is DOUBLE PRECISION array, dimension ( LDC, N )
1.8 bertrand 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 symmetric 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 symmetric 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: *> \endverbatim
156: *>
157: *> \param[in] LDC
158: *> \verbatim
159: *> LDC is INTEGER
160: *> On entry, LDC specifies the first dimension of C as declared
161: *> in the calling (sub) program. LDC must be at least
162: *> max( 1, n ).
163: *> \endverbatim
164: *
165: * Authors:
166: * ========
167: *
1.13 bertrand 168: *> \author Univ. of Tennessee
169: *> \author Univ. of California Berkeley
170: *> \author Univ. of Colorado Denver
171: *> \author NAG Ltd.
1.8 bertrand 172: *
173: *> \ingroup double_blas_level3
174: *
175: *> \par Further Details:
176: * =====================
177: *>
178: *> \verbatim
179: *>
180: *> Level 3 Blas routine.
181: *>
182: *>
183: *> -- Written on 8-February-1989.
184: *> Jack Dongarra, Argonne National Laboratory.
185: *> Iain Duff, AERE Harwell.
186: *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
187: *> Sven Hammarling, Numerical Algorithms Group Ltd.
188: *> \endverbatim
189: *>
190: * =====================================================================
1.1 bertrand 191: SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
1.8 bertrand 192: *
1.16 ! bertrand 193: * -- Reference BLAS level3 routine --
1.8 bertrand 194: * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
195: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196: *
1.1 bertrand 197: * .. Scalar Arguments ..
198: DOUBLE PRECISION ALPHA,BETA
199: INTEGER K,LDA,LDB,LDC,N
200: CHARACTER TRANS,UPLO
201: * ..
202: * .. Array Arguments ..
203: DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
204: * ..
205: *
206: * =====================================================================
207: *
208: * .. External Functions ..
209: LOGICAL LSAME
210: EXTERNAL LSAME
211: * ..
212: * .. External Subroutines ..
213: EXTERNAL XERBLA
214: * ..
215: * .. Intrinsic Functions ..
216: INTRINSIC MAX
217: * ..
218: * .. Local Scalars ..
219: DOUBLE PRECISION TEMP1,TEMP2
220: INTEGER I,INFO,J,L,NROWA
221: LOGICAL UPPER
222: * ..
223: * .. Parameters ..
224: DOUBLE PRECISION ONE,ZERO
225: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
226: * ..
227: *
228: * Test the input parameters.
229: *
230: IF (LSAME(TRANS,'N')) THEN
231: NROWA = N
232: ELSE
233: NROWA = K
234: END IF
235: UPPER = LSAME(UPLO,'U')
236: *
237: INFO = 0
238: IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
239: INFO = 1
240: ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
241: + (.NOT.LSAME(TRANS,'T')) .AND.
242: + (.NOT.LSAME(TRANS,'C'))) THEN
243: INFO = 2
244: ELSE IF (N.LT.0) THEN
245: INFO = 3
246: ELSE IF (K.LT.0) THEN
247: INFO = 4
248: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
249: INFO = 7
250: ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
251: INFO = 9
252: ELSE IF (LDC.LT.MAX(1,N)) THEN
253: INFO = 12
254: END IF
255: IF (INFO.NE.0) THEN
256: CALL XERBLA('DSYR2K',INFO)
257: RETURN
258: END IF
259: *
260: * Quick return if possible.
261: *
262: IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
263: + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
264: *
265: * And when alpha.eq.zero.
266: *
267: IF (ALPHA.EQ.ZERO) THEN
268: IF (UPPER) THEN
269: IF (BETA.EQ.ZERO) THEN
270: DO 20 J = 1,N
271: DO 10 I = 1,J
272: C(I,J) = ZERO
273: 10 CONTINUE
274: 20 CONTINUE
275: ELSE
276: DO 40 J = 1,N
277: DO 30 I = 1,J
278: C(I,J) = BETA*C(I,J)
279: 30 CONTINUE
280: 40 CONTINUE
281: END IF
282: ELSE
283: IF (BETA.EQ.ZERO) THEN
284: DO 60 J = 1,N
285: DO 50 I = J,N
286: C(I,J) = ZERO
287: 50 CONTINUE
288: 60 CONTINUE
289: ELSE
290: DO 80 J = 1,N
291: DO 70 I = J,N
292: C(I,J) = BETA*C(I,J)
293: 70 CONTINUE
294: 80 CONTINUE
295: END IF
296: END IF
297: RETURN
298: END IF
299: *
300: * Start the operations.
301: *
302: IF (LSAME(TRANS,'N')) THEN
303: *
1.7 bertrand 304: * Form C := alpha*A*B**T + alpha*B*A**T + C.
1.1 bertrand 305: *
306: IF (UPPER) THEN
307: DO 130 J = 1,N
308: IF (BETA.EQ.ZERO) THEN
309: DO 90 I = 1,J
310: C(I,J) = ZERO
311: 90 CONTINUE
312: ELSE IF (BETA.NE.ONE) THEN
313: DO 100 I = 1,J
314: C(I,J) = BETA*C(I,J)
315: 100 CONTINUE
316: END IF
317: DO 120 L = 1,K
318: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
319: TEMP1 = ALPHA*B(J,L)
320: TEMP2 = ALPHA*A(J,L)
321: DO 110 I = 1,J
322: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
323: + B(I,L)*TEMP2
324: 110 CONTINUE
325: END IF
326: 120 CONTINUE
327: 130 CONTINUE
328: ELSE
329: DO 180 J = 1,N
330: IF (BETA.EQ.ZERO) THEN
331: DO 140 I = J,N
332: C(I,J) = ZERO
333: 140 CONTINUE
334: ELSE IF (BETA.NE.ONE) THEN
335: DO 150 I = J,N
336: C(I,J) = BETA*C(I,J)
337: 150 CONTINUE
338: END IF
339: DO 170 L = 1,K
340: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
341: TEMP1 = ALPHA*B(J,L)
342: TEMP2 = ALPHA*A(J,L)
343: DO 160 I = J,N
344: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
345: + B(I,L)*TEMP2
346: 160 CONTINUE
347: END IF
348: 170 CONTINUE
349: 180 CONTINUE
350: END IF
351: ELSE
352: *
1.7 bertrand 353: * Form C := alpha*A**T*B + alpha*B**T*A + C.
1.1 bertrand 354: *
355: IF (UPPER) THEN
356: DO 210 J = 1,N
357: DO 200 I = 1,J
358: TEMP1 = ZERO
359: TEMP2 = ZERO
360: DO 190 L = 1,K
361: TEMP1 = TEMP1 + A(L,I)*B(L,J)
362: TEMP2 = TEMP2 + B(L,I)*A(L,J)
363: 190 CONTINUE
364: IF (BETA.EQ.ZERO) THEN
365: C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
366: ELSE
367: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
368: + ALPHA*TEMP2
369: END IF
370: 200 CONTINUE
371: 210 CONTINUE
372: ELSE
373: DO 240 J = 1,N
374: DO 230 I = J,N
375: TEMP1 = ZERO
376: TEMP2 = ZERO
377: DO 220 L = 1,K
378: TEMP1 = TEMP1 + A(L,I)*B(L,J)
379: TEMP2 = TEMP2 + B(L,I)*A(L,J)
380: 220 CONTINUE
381: IF (BETA.EQ.ZERO) THEN
382: C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
383: ELSE
384: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
385: + ALPHA*TEMP2
386: END IF
387: 230 CONTINUE
388: 240 CONTINUE
389: END IF
390: END IF
391: *
392: RETURN
393: *
1.16 ! bertrand 394: * End of DSYR2K
1.1 bertrand 395: *
396: END