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