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