![[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|
Return to dsfrk.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
2: + C )
3: *
1.4 ! bertrand 4: * -- LAPACK routine (version 3.3.0) --
1.1 bertrand 5: *
6: * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
1.4 ! bertrand 7: * November 2010
1.1 bertrand 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' + beta*C,
30: *
31: * or
32: *
33: * C := alpha*A'*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: *
1.4 ! bertrand 42: * TRANSR (input) CHARACTER*1
1.1 bertrand 43: * = 'N': The Normal Form of RFP A is stored;
44: * = 'T': The Transpose Form of RFP A is stored.
45: *
1.4 ! bertrand 46: * UPLO (input) CHARACTER*1
1.1 bertrand 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: *
1.4 ! bertrand 59: * TRANS (input) CHARACTER*1
1.1 bertrand 60: * On entry, TRANS specifies the operation to be performed as
61: * follows:
62: *
63: * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
64: *
65: * TRANS = 'T' or 't' C := alpha*A'*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: * Arguments
110: * ==========
111: *
112: * ..
113: * .. Parameters ..
114: DOUBLE PRECISION ONE, ZERO
115: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
116: * ..
117: * .. Local Scalars ..
118: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
119: INTEGER INFO, NROWA, J, NK, N1, N2
120: * ..
121: * .. External Functions ..
122: LOGICAL LSAME
123: EXTERNAL LSAME
124: * ..
125: * .. External Subroutines ..
126: EXTERNAL XERBLA, DGEMM, DSYRK
127: * ..
128: * .. Intrinsic Functions ..
129: INTRINSIC MAX
130: * ..
131: * .. Executable Statements ..
132: *
133: * Test the input parameters.
134: *
135: INFO = 0
136: NORMALTRANSR = LSAME( TRANSR, 'N' )
137: LOWER = LSAME( UPLO, 'L' )
138: NOTRANS = LSAME( TRANS, 'N' )
139: *
140: IF( NOTRANS ) THEN
141: NROWA = N
142: ELSE
143: NROWA = K
144: END IF
145: *
146: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
147: INFO = -1
148: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
149: INFO = -2
150: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
151: INFO = -3
152: ELSE IF( N.LT.0 ) THEN
153: INFO = -4
154: ELSE IF( K.LT.0 ) THEN
155: INFO = -5
156: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
157: INFO = -8
158: END IF
159: IF( INFO.NE.0 ) THEN
160: CALL XERBLA( 'DSFRK ', -INFO )
161: RETURN
162: END IF
163: *
164: * Quick return if possible.
165: *
166: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
167: * done (it is in DSYRK for example) and left in the general case.
168: *
169: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
170: + ( BETA.EQ.ONE ) ) )RETURN
171: *
172: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
173: DO J = 1, ( ( N*( N+1 ) ) / 2 )
174: C( J ) = ZERO
175: END DO
176: RETURN
177: END IF
178: *
179: * C is N-by-N.
180: * If N is odd, set NISODD = .TRUE., and N1 and N2.
181: * If N is even, NISODD = .FALSE., and NK.
182: *
183: IF( MOD( N, 2 ).EQ.0 ) THEN
184: NISODD = .FALSE.
185: NK = N / 2
186: ELSE
187: NISODD = .TRUE.
188: IF( LOWER ) THEN
189: N2 = N / 2
190: N1 = N - N2
191: ELSE
192: N1 = N / 2
193: N2 = N - N1
194: END IF
195: END IF
196: *
197: IF( NISODD ) THEN
198: *
199: * N is odd
200: *
201: IF( NORMALTRANSR ) THEN
202: *
203: * N is odd and TRANSR = 'N'
204: *
205: IF( LOWER ) THEN
206: *
207: * N is odd, TRANSR = 'N', and UPLO = 'L'
208: *
209: IF( NOTRANS ) THEN
210: *
211: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
212: *
213: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
214: + BETA, C( 1 ), N )
215: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
216: + BETA, C( N+1 ), N )
217: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
218: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
219: *
220: ELSE
221: *
222: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
223: *
224: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
225: + BETA, C( 1 ), N )
226: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
227: + BETA, C( N+1 ), N )
228: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
229: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
230: *
231: END IF
232: *
233: ELSE
234: *
235: * N is odd, TRANSR = 'N', and UPLO = 'U'
236: *
237: IF( NOTRANS ) THEN
238: *
239: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
240: *
241: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
242: + BETA, C( N2+1 ), N )
243: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
244: + BETA, C( N1+1 ), N )
245: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
246: + LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
247: *
248: ELSE
249: *
250: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
251: *
252: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
253: + BETA, C( N2+1 ), N )
254: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
255: + BETA, C( N1+1 ), N )
256: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
257: + LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
258: *
259: END IF
260: *
261: END IF
262: *
263: ELSE
264: *
265: * N is odd, and TRANSR = 'T'
266: *
267: IF( LOWER ) THEN
268: *
269: * N is odd, TRANSR = 'T', and UPLO = 'L'
270: *
271: IF( NOTRANS ) THEN
272: *
273: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
274: *
275: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
276: + BETA, C( 1 ), N1 )
277: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
278: + BETA, C( 2 ), N1 )
279: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
280: + LDA, A( N1+1, 1 ), LDA, BETA,
281: + C( N1*N1+1 ), N1 )
282: *
283: ELSE
284: *
285: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
286: *
287: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
288: + BETA, C( 1 ), N1 )
289: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
290: + BETA, C( 2 ), N1 )
291: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
292: + LDA, A( 1, N1+1 ), LDA, BETA,
293: + C( N1*N1+1 ), N1 )
294: *
295: END IF
296: *
297: ELSE
298: *
299: * N is odd, TRANSR = 'T', and UPLO = 'U'
300: *
301: IF( NOTRANS ) THEN
302: *
303: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
304: *
305: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
306: + BETA, C( N2*N2+1 ), N2 )
307: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
308: + BETA, C( N1*N2+1 ), N2 )
309: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
310: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
311: *
312: ELSE
313: *
314: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
315: *
316: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
317: + BETA, C( N2*N2+1 ), N2 )
318: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
319: + BETA, C( N1*N2+1 ), N2 )
320: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
321: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
322: *
323: END IF
324: *
325: END IF
326: *
327: END IF
328: *
329: ELSE
330: *
331: * N is even
332: *
333: IF( NORMALTRANSR ) THEN
334: *
335: * N is even and TRANSR = 'N'
336: *
337: IF( LOWER ) THEN
338: *
339: * N is even, TRANSR = 'N', and UPLO = 'L'
340: *
341: IF( NOTRANS ) THEN
342: *
343: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
344: *
345: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
346: + BETA, C( 2 ), N+1 )
347: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
348: + BETA, C( 1 ), N+1 )
349: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
350: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
351: + N+1 )
352: *
353: ELSE
354: *
355: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
356: *
357: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
358: + BETA, C( 2 ), N+1 )
359: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
360: + BETA, C( 1 ), N+1 )
361: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
362: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
363: + N+1 )
364: *
365: END IF
366: *
367: ELSE
368: *
369: * N is even, TRANSR = 'N', and UPLO = 'U'
370: *
371: IF( NOTRANS ) THEN
372: *
373: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
374: *
375: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
376: + BETA, C( NK+2 ), N+1 )
377: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
378: + BETA, C( NK+1 ), N+1 )
379: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
380: + LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
381: + N+1 )
382: *
383: ELSE
384: *
385: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
386: *
387: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
388: + BETA, C( NK+2 ), N+1 )
389: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
390: + BETA, C( NK+1 ), N+1 )
391: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
392: + LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
393: + N+1 )
394: *
395: END IF
396: *
397: END IF
398: *
399: ELSE
400: *
401: * N is even, and TRANSR = 'T'
402: *
403: IF( LOWER ) THEN
404: *
405: * N is even, TRANSR = 'T', and UPLO = 'L'
406: *
407: IF( NOTRANS ) THEN
408: *
409: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
410: *
411: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
412: + BETA, C( NK+1 ), NK )
413: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
414: + BETA, C( 1 ), NK )
415: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
416: + LDA, A( NK+1, 1 ), LDA, BETA,
417: + C( ( ( NK+1 )*NK )+1 ), NK )
418: *
419: ELSE
420: *
421: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
422: *
423: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
424: + BETA, C( NK+1 ), NK )
425: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
426: + BETA, C( 1 ), NK )
427: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
428: + LDA, A( 1, NK+1 ), LDA, BETA,
429: + C( ( ( NK+1 )*NK )+1 ), NK )
430: *
431: END IF
432: *
433: ELSE
434: *
435: * N is even, TRANSR = 'T', and UPLO = 'U'
436: *
437: IF( NOTRANS ) THEN
438: *
439: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
440: *
441: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
442: + BETA, C( NK*( NK+1 )+1 ), NK )
443: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
444: + BETA, C( NK*NK+1 ), NK )
445: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
446: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
447: *
448: ELSE
449: *
450: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
451: *
452: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
453: + BETA, C( NK*( NK+1 )+1 ), NK )
454: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
455: + BETA, C( NK*NK+1 ), NK )
456: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
457: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
458: *
459: END IF
460: *
461: END IF
462: *
463: END IF
464: *
465: END IF
466: *
467: RETURN
468: *
469: * End of DSFRK
470: *
471: END