1: *> \brief \b DSFRK
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: *> \date November 2011
162: *
163: *> \ingroup doubleOTHERcomputational
164: *
165: * =====================================================================
166: SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
167: $ C )
168: *
169: * -- LAPACK computational routine (version 3.4.0) --
170: * -- LAPACK is a software package provided by Univ. of Tennessee, --
171: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172: * November 2011
173: *
174: * .. Scalar Arguments ..
175: DOUBLE PRECISION ALPHA, BETA
176: INTEGER K, LDA, N
177: CHARACTER TRANS, TRANSR, UPLO
178: * ..
179: * .. Array Arguments ..
180: DOUBLE PRECISION A( LDA, * ), C( * )
181: * ..
182: *
183: * =====================================================================
184: *
185: * ..
186: * .. Parameters ..
187: DOUBLE PRECISION ONE, ZERO
188: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
189: * ..
190: * .. Local Scalars ..
191: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
192: INTEGER INFO, NROWA, J, NK, N1, N2
193: * ..
194: * .. External Functions ..
195: LOGICAL LSAME
196: EXTERNAL LSAME
197: * ..
198: * .. External Subroutines ..
199: EXTERNAL XERBLA, DGEMM, DSYRK
200: * ..
201: * .. Intrinsic Functions ..
202: INTRINSIC MAX
203: * ..
204: * .. Executable Statements ..
205: *
206: * Test the input parameters.
207: *
208: INFO = 0
209: NORMALTRANSR = LSAME( TRANSR, 'N' )
210: LOWER = LSAME( UPLO, 'L' )
211: NOTRANS = LSAME( TRANS, 'N' )
212: *
213: IF( NOTRANS ) THEN
214: NROWA = N
215: ELSE
216: NROWA = K
217: END IF
218: *
219: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
220: INFO = -1
221: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
222: INFO = -2
223: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
224: INFO = -3
225: ELSE IF( N.LT.0 ) THEN
226: INFO = -4
227: ELSE IF( K.LT.0 ) THEN
228: INFO = -5
229: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
230: INFO = -8
231: END IF
232: IF( INFO.NE.0 ) THEN
233: CALL XERBLA( 'DSFRK ', -INFO )
234: RETURN
235: END IF
236: *
237: * Quick return if possible.
238: *
239: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
240: * done (it is in DSYRK for example) and left in the general case.
241: *
242: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
243: $ ( BETA.EQ.ONE ) ) )RETURN
244: *
245: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
246: DO J = 1, ( ( N*( N+1 ) ) / 2 )
247: C( J ) = ZERO
248: END DO
249: RETURN
250: END IF
251: *
252: * C is N-by-N.
253: * If N is odd, set NISODD = .TRUE., and N1 and N2.
254: * If N is even, NISODD = .FALSE., and NK.
255: *
256: IF( MOD( N, 2 ).EQ.0 ) THEN
257: NISODD = .FALSE.
258: NK = N / 2
259: ELSE
260: NISODD = .TRUE.
261: IF( LOWER ) THEN
262: N2 = N / 2
263: N1 = N - N2
264: ELSE
265: N1 = N / 2
266: N2 = N - N1
267: END IF
268: END IF
269: *
270: IF( NISODD ) THEN
271: *
272: * N is odd
273: *
274: IF( NORMALTRANSR ) THEN
275: *
276: * N is odd and TRANSR = 'N'
277: *
278: IF( LOWER ) THEN
279: *
280: * N is odd, TRANSR = 'N', and UPLO = 'L'
281: *
282: IF( NOTRANS ) THEN
283: *
284: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
285: *
286: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
287: $ BETA, C( 1 ), N )
288: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
289: $ BETA, C( N+1 ), N )
290: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
291: $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
292: *
293: ELSE
294: *
295: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
296: *
297: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
298: $ BETA, C( 1 ), N )
299: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
300: $ BETA, C( N+1 ), N )
301: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
302: $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
303: *
304: END IF
305: *
306: ELSE
307: *
308: * N is odd, TRANSR = 'N', and UPLO = 'U'
309: *
310: IF( NOTRANS ) THEN
311: *
312: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
313: *
314: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
315: $ BETA, C( N2+1 ), N )
316: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
317: $ BETA, C( N1+1 ), N )
318: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
319: $ LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
320: *
321: ELSE
322: *
323: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
324: *
325: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
326: $ BETA, C( N2+1 ), N )
327: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
328: $ BETA, C( N1+1 ), N )
329: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
330: $ LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
331: *
332: END IF
333: *
334: END IF
335: *
336: ELSE
337: *
338: * N is odd, and TRANSR = 'T'
339: *
340: IF( LOWER ) THEN
341: *
342: * N is odd, TRANSR = 'T', and UPLO = 'L'
343: *
344: IF( NOTRANS ) THEN
345: *
346: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
347: *
348: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
349: $ BETA, C( 1 ), N1 )
350: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
351: $ BETA, C( 2 ), N1 )
352: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
353: $ LDA, A( N1+1, 1 ), LDA, BETA,
354: $ C( N1*N1+1 ), N1 )
355: *
356: ELSE
357: *
358: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
359: *
360: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
361: $ BETA, C( 1 ), N1 )
362: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
363: $ BETA, C( 2 ), N1 )
364: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
365: $ LDA, A( 1, N1+1 ), LDA, BETA,
366: $ C( N1*N1+1 ), N1 )
367: *
368: END IF
369: *
370: ELSE
371: *
372: * N is odd, TRANSR = 'T', and UPLO = 'U'
373: *
374: IF( NOTRANS ) THEN
375: *
376: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
377: *
378: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
379: $ BETA, C( N2*N2+1 ), N2 )
380: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
381: $ BETA, C( N1*N2+1 ), N2 )
382: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
383: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
384: *
385: ELSE
386: *
387: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
388: *
389: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
390: $ BETA, C( N2*N2+1 ), N2 )
391: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
392: $ BETA, C( N1*N2+1 ), N2 )
393: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
394: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
395: *
396: END IF
397: *
398: END IF
399: *
400: END IF
401: *
402: ELSE
403: *
404: * N is even
405: *
406: IF( NORMALTRANSR ) THEN
407: *
408: * N is even and TRANSR = 'N'
409: *
410: IF( LOWER ) THEN
411: *
412: * N is even, TRANSR = 'N', and UPLO = 'L'
413: *
414: IF( NOTRANS ) THEN
415: *
416: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
417: *
418: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
419: $ BETA, C( 2 ), N+1 )
420: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
421: $ BETA, C( 1 ), N+1 )
422: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
423: $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
424: $ N+1 )
425: *
426: ELSE
427: *
428: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
429: *
430: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
431: $ BETA, C( 2 ), N+1 )
432: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
433: $ BETA, C( 1 ), N+1 )
434: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
435: $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
436: $ N+1 )
437: *
438: END IF
439: *
440: ELSE
441: *
442: * N is even, TRANSR = 'N', and UPLO = 'U'
443: *
444: IF( NOTRANS ) THEN
445: *
446: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
447: *
448: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
449: $ BETA, C( NK+2 ), N+1 )
450: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
451: $ BETA, C( NK+1 ), N+1 )
452: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
453: $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
454: $ N+1 )
455: *
456: ELSE
457: *
458: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
459: *
460: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
461: $ BETA, C( NK+2 ), N+1 )
462: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
463: $ BETA, C( NK+1 ), N+1 )
464: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
465: $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
466: $ N+1 )
467: *
468: END IF
469: *
470: END IF
471: *
472: ELSE
473: *
474: * N is even, and TRANSR = 'T'
475: *
476: IF( LOWER ) THEN
477: *
478: * N is even, TRANSR = 'T', and UPLO = 'L'
479: *
480: IF( NOTRANS ) THEN
481: *
482: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
483: *
484: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
485: $ BETA, C( NK+1 ), NK )
486: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
487: $ BETA, C( 1 ), NK )
488: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
489: $ LDA, A( NK+1, 1 ), LDA, BETA,
490: $ C( ( ( NK+1 )*NK )+1 ), NK )
491: *
492: ELSE
493: *
494: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
495: *
496: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
497: $ BETA, C( NK+1 ), NK )
498: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
499: $ BETA, C( 1 ), NK )
500: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
501: $ LDA, A( 1, NK+1 ), LDA, BETA,
502: $ C( ( ( NK+1 )*NK )+1 ), NK )
503: *
504: END IF
505: *
506: ELSE
507: *
508: * N is even, TRANSR = 'T', and UPLO = 'U'
509: *
510: IF( NOTRANS ) THEN
511: *
512: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
513: *
514: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
515: $ BETA, C( NK*( NK+1 )+1 ), NK )
516: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
517: $ BETA, C( NK*NK+1 ), NK )
518: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
519: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
520: *
521: ELSE
522: *
523: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
524: *
525: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
526: $ BETA, C( NK*( NK+1 )+1 ), NK )
527: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
528: $ BETA, C( NK*NK+1 ), NK )
529: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
530: $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
531: *
532: END IF
533: *
534: END IF
535: *
536: END IF
537: *
538: END IF
539: *
540: RETURN
541: *
542: * End of DSFRK
543: *
544: 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