1: SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
2: $ C )
3: *
4: * -- LAPACK routine (version 3.3.1) --
5: *
6: * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
7: * -- April 2011 --
8: *
9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11: *
12: * ..
13: * .. Scalar Arguments ..
14: DOUBLE PRECISION ALPHA, BETA
15: INTEGER K, LDA, N
16: CHARACTER TRANS, TRANSR, UPLO
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION A( LDA, * ), C( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * Level 3 BLAS like routine for C in RFP Format.
26: *
27: * DSFRK performs one of the symmetric rank--k operations
28: *
29: * C := alpha*A*A**T + beta*C,
30: *
31: * or
32: *
33: * C := alpha*A**T*A + beta*C,
34: *
35: * where alpha and beta are real scalars, C is an n--by--n symmetric
36: * matrix and A is an n--by--k matrix in the first case and a k--by--n
37: * matrix in the second case.
38: *
39: * Arguments
40: * ==========
41: *
42: * TRANSR (input) CHARACTER*1
43: * = 'N': The Normal Form of RFP A is stored;
44: * = 'T': The Transpose Form of RFP A is stored.
45: *
46: * UPLO (input) 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: *
57: * Unchanged on exit.
58: *
59: * TRANS (input) CHARACTER*1
60: * On entry, TRANS specifies the operation to be performed as
61: * follows:
62: *
63: * TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C.
64: *
65: * TRANS = 'T' or 't' C := alpha*A**T*A + beta*C.
66: *
67: * Unchanged on exit.
68: *
69: * N (input) INTEGER
70: * On entry, N specifies the order of the matrix C. N must be
71: * at least zero.
72: * Unchanged on exit.
73: *
74: * K (input) INTEGER
75: * On entry with TRANS = 'N' or 'n', K specifies the number
76: * of columns of the matrix A, and on entry with TRANS = 'T'
77: * or 't', K specifies the number of rows of the matrix A. K
78: * must be at least zero.
79: * Unchanged on exit.
80: *
81: * ALPHA (input) DOUBLE PRECISION
82: * On entry, ALPHA specifies the scalar alpha.
83: * Unchanged on exit.
84: *
85: * A (input) DOUBLE PRECISION array, dimension (LDA,ka)
86: * where KA
87: * is K when TRANS = 'N' or 'n', and is N otherwise. Before
88: * entry with TRANS = 'N' or 'n', the leading N--by--K part of
89: * the array A must contain the matrix A, otherwise the leading
90: * K--by--N part of the array A must contain the matrix A.
91: * Unchanged on exit.
92: *
93: * LDA (input) INTEGER
94: * On entry, LDA specifies the first dimension of A as declared
95: * in the calling (sub) program. When TRANS = 'N' or 'n'
96: * then LDA must be at least max( 1, n ), otherwise LDA must
97: * be at least max( 1, k ).
98: * Unchanged on exit.
99: *
100: * BETA (input) DOUBLE PRECISION
101: * On entry, BETA specifies the scalar beta.
102: * Unchanged on exit.
103: *
104: *
105: * C (input/output) DOUBLE PRECISION array, dimension (NT)
106: * NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
107: * Format. RFP Format is described by TRANSR, UPLO and N.
108: *
109: * =====================================================================
110: *
111: * ..
112: * .. Parameters ..
113: DOUBLE PRECISION ONE, ZERO
114: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
115: * ..
116: * .. Local Scalars ..
117: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
118: INTEGER INFO, NROWA, J, NK, N1, N2
119: * ..
120: * .. External Functions ..
121: LOGICAL LSAME
122: EXTERNAL LSAME
123: * ..
124: * .. External Subroutines ..
125: EXTERNAL XERBLA, DGEMM, DSYRK
126: * ..
127: * .. Intrinsic Functions ..
128: INTRINSIC MAX
129: * ..
130: * .. Executable Statements ..
131: *
132: * Test the input parameters.
133: *
134: INFO = 0
135: NORMALTRANSR = LSAME( TRANSR, 'N' )
136: LOWER = LSAME( UPLO, 'L' )
137: NOTRANS = LSAME( TRANS, 'N' )
138: *
139: IF( NOTRANS ) THEN
140: NROWA = N
141: ELSE
142: NROWA = K
143: END IF
144: *
145: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
146: INFO = -1
147: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
148: INFO = -2
149: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
150: INFO = -3
151: ELSE IF( N.LT.0 ) THEN
152: INFO = -4
153: ELSE IF( K.LT.0 ) THEN
154: INFO = -5
155: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
156: INFO = -8
157: END IF
158: IF( INFO.NE.0 ) THEN
159: CALL XERBLA( 'DSFRK ', -INFO )
160: RETURN
161: END IF
162: *
163: * Quick return if possible.
164: *
165: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
166: * done (it is in DSYRK for example) and left in the general case.
167: *
168: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
169: $ ( BETA.EQ.ONE ) ) )RETURN
170: *
171: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
172: DO J = 1, ( ( N*( N+1 ) ) / 2 )
173: C( J ) = ZERO
174: END DO
175: RETURN
176: END IF
177: *
178: * C is N-by-N.
179: * If N is odd, set NISODD = .TRUE., and N1 and N2.
180: * If N is even, NISODD = .FALSE., and NK.
181: *
182: IF( MOD( N, 2 ).EQ.0 ) THEN
183: NISODD = .FALSE.
184: NK = N / 2
185: ELSE
186: NISODD = .TRUE.
187: IF( LOWER ) THEN
188: N2 = N / 2
189: N1 = N - N2
190: ELSE
191: N1 = N / 2
192: N2 = N - N1
193: END IF
194: END IF
195: *
196: IF( NISODD ) THEN
197: *
198: * N is odd
199: *
200: IF( NORMALTRANSR ) THEN
201: *
202: * N is odd and TRANSR = 'N'
203: *
204: IF( LOWER ) THEN
205: *
206: * N is odd, TRANSR = 'N', and UPLO = 'L'
207: *
208: IF( NOTRANS ) THEN
209: *
210: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
211: *
212: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
213: $ BETA, C( 1 ), N )
214: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
215: $ BETA, C( N+1 ), N )
216: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
217: $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
218: *
219: ELSE
220: *
221: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
222: *
223: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
224: $ BETA, C( 1 ), N )
225: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
226: $ BETA, C( N+1 ), N )
227: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
228: $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
229: *
230: END IF
231: *
232: ELSE
233: *
234: * N is odd, TRANSR = 'N', and UPLO = 'U'
235: *
236: IF( NOTRANS ) THEN
237: *
238: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
239: *
240: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
241: $ BETA, C( N2+1 ), N )
242: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
243: $ BETA, C( N1+1 ), N )
244: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
245: $ LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
246: *
247: ELSE
248: *
249: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
250: *
251: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
252: $ BETA, C( N2+1 ), N )
253: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
254: $ BETA, C( N1+1 ), N )
255: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
256: $ LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
257: *
258: END IF
259: *
260: END IF
261: *
262: ELSE
263: *
264: * N is odd, and TRANSR = 'T'
265: *
266: IF( LOWER ) THEN
267: *
268: * N is odd, TRANSR = 'T', and UPLO = 'L'
269: *
270: IF( NOTRANS ) THEN
271: *
272: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
273: *
274: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
275: $ BETA, C( 1 ), N1 )
276: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
277: $ BETA, C( 2 ), N1 )
278: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
279: $ LDA, A( N1+1, 1 ), LDA, BETA,
280: $ C( N1*N1+1 ), N1 )
281: *
282: ELSE
283: *
284: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
285: *
286: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
287: $ BETA, C( 1 ), N1 )
288: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
289: $ BETA, C( 2 ), N1 )
290: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
291: $ LDA, A( 1, N1+1 ), LDA, BETA,
292: $ C( N1*N1+1 ), N1 )
293: *
294: END IF
295: *
296: ELSE
297: *
298: * N is odd, TRANSR = 'T', and UPLO = 'U'
299: *
300: IF( NOTRANS ) THEN
301: *
302: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
303: *
304: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
305: $ BETA, C( N2*N2+1 ), N2 )
306: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
307: $ BETA, C( N1*N2+1 ), N2 )
308: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
309: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
310: *
311: ELSE
312: *
313: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
314: *
315: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
316: $ BETA, C( N2*N2+1 ), N2 )
317: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
318: $ BETA, C( N1*N2+1 ), N2 )
319: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
320: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
321: *
322: END IF
323: *
324: END IF
325: *
326: END IF
327: *
328: ELSE
329: *
330: * N is even
331: *
332: IF( NORMALTRANSR ) THEN
333: *
334: * N is even and TRANSR = 'N'
335: *
336: IF( LOWER ) THEN
337: *
338: * N is even, TRANSR = 'N', and UPLO = 'L'
339: *
340: IF( NOTRANS ) THEN
341: *
342: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
343: *
344: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
345: $ BETA, C( 2 ), N+1 )
346: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
347: $ BETA, C( 1 ), N+1 )
348: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
349: $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
350: $ N+1 )
351: *
352: ELSE
353: *
354: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
355: *
356: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
357: $ BETA, C( 2 ), N+1 )
358: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
359: $ BETA, C( 1 ), N+1 )
360: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
361: $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
362: $ N+1 )
363: *
364: END IF
365: *
366: ELSE
367: *
368: * N is even, TRANSR = 'N', and UPLO = 'U'
369: *
370: IF( NOTRANS ) THEN
371: *
372: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
373: *
374: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
375: $ BETA, C( NK+2 ), N+1 )
376: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
377: $ BETA, C( NK+1 ), N+1 )
378: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
379: $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
380: $ N+1 )
381: *
382: ELSE
383: *
384: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
385: *
386: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
387: $ BETA, C( NK+2 ), N+1 )
388: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
389: $ BETA, C( NK+1 ), N+1 )
390: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
391: $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
392: $ N+1 )
393: *
394: END IF
395: *
396: END IF
397: *
398: ELSE
399: *
400: * N is even, and TRANSR = 'T'
401: *
402: IF( LOWER ) THEN
403: *
404: * N is even, TRANSR = 'T', and UPLO = 'L'
405: *
406: IF( NOTRANS ) THEN
407: *
408: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
409: *
410: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
411: $ BETA, C( NK+1 ), NK )
412: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
413: $ BETA, C( 1 ), NK )
414: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
415: $ LDA, A( NK+1, 1 ), LDA, BETA,
416: $ C( ( ( NK+1 )*NK )+1 ), NK )
417: *
418: ELSE
419: *
420: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
421: *
422: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
423: $ BETA, C( NK+1 ), NK )
424: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
425: $ BETA, C( 1 ), NK )
426: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
427: $ LDA, A( 1, NK+1 ), LDA, BETA,
428: $ C( ( ( NK+1 )*NK )+1 ), NK )
429: *
430: END IF
431: *
432: ELSE
433: *
434: * N is even, TRANSR = 'T', and UPLO = 'U'
435: *
436: IF( NOTRANS ) THEN
437: *
438: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
439: *
440: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
441: $ BETA, C( NK*( NK+1 )+1 ), NK )
442: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
443: $ BETA, C( NK*NK+1 ), NK )
444: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
445: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
446: *
447: ELSE
448: *
449: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
450: *
451: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
452: $ BETA, C( NK*( NK+1 )+1 ), NK )
453: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
454: $ BETA, C( NK*NK+1 ), NK )
455: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
456: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
457: *
458: END IF
459: *
460: END IF
461: *
462: END IF
463: *
464: END IF
465: *
466: RETURN
467: *
468: * End of DSFRK
469: *
470: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsfrk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack