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