LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssfrk.f
Go to the documentation of this file.
1*> \brief \b SSFRK 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 SSFRK + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssfrk.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssfrk.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssfrk.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
20* C )
21*
22* .. Scalar Arguments ..
23* REAL ALPHA, BETA
24* INTEGER K, LDA, N
25* CHARACTER TRANS, TRANSR, UPLO
26* ..
27* .. Array Arguments ..
28* REAL 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*> SSFRK 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 REAL
112*> On entry, ALPHA specifies the scalar alpha.
113*> Unchanged on exit.
114*> \endverbatim
115*>
116*> \param[in] A
117*> \verbatim
118*> A is REAL 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 REAL
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 REAL 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 ssfrk( 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 REAL ALPHA, BETA
172 INTEGER K, LDA, N
173 CHARACTER TRANS, TRANSR, UPLO
174* ..
175* .. Array Arguments ..
176 REAL A( LDA, * ), C( * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 REAL ONE, ZERO
183 PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
184* ..
185* .. Local Scalars ..
186 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
187 INTEGER INFO, NROWA, J, NK, N1, N2
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 EXTERNAL LSAME
192* ..
193* .. External Subroutines ..
194 EXTERNAL sgemm, ssyrk, xerbla
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC max
198* ..
199* .. Executable Statements ..
200*
201* Test the input parameters.
202*
203 info = 0
204 normaltransr = lsame( transr, 'N' )
205 lower = lsame( uplo, 'L' )
206 notrans = lsame( trans, 'N' )
207*
208 IF( notrans ) THEN
209 nrowa = n
210 ELSE
211 nrowa = k
212 END IF
213*
214 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
215 info = -1
216 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
217 info = -2
218 ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'T' ) ) THEN
219 info = -3
220 ELSE IF( n.LT.0 ) THEN
221 info = -4
222 ELSE IF( k.LT.0 ) THEN
223 info = -5
224 ELSE IF( lda.LT.max( 1, nrowa ) ) THEN
225 info = -8
226 END IF
227 IF( info.NE.0 ) THEN
228 CALL xerbla( 'SSFRK ', -info )
229 RETURN
230 END IF
231*
232* Quick return if possible.
233*
234* The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
235* done (it is in SSYRK for example) and left in the general case.
236*
237 IF( ( n.EQ.0 ) .OR. ( ( ( alpha.EQ.zero ) .OR. ( k.EQ.0 ) ) .AND.
238 $ ( beta.EQ.one ) ) )RETURN
239*
240 IF( ( alpha.EQ.zero ) .AND. ( beta.EQ.zero ) ) THEN
241 DO j = 1, ( ( n*( n+1 ) ) / 2 )
242 c( j ) = zero
243 END DO
244 RETURN
245 END IF
246*
247* C is N-by-N.
248* If N is odd, set NISODD = .TRUE., and N1 and N2.
249* If N is even, NISODD = .FALSE., and NK.
250*
251 IF( mod( n, 2 ).EQ.0 ) THEN
252 nisodd = .false.
253 nk = n / 2
254 ELSE
255 nisodd = .true.
256 IF( lower ) THEN
257 n2 = n / 2
258 n1 = n - n2
259 ELSE
260 n1 = n / 2
261 n2 = n - n1
262 END IF
263 END IF
264*
265 IF( nisodd ) THEN
266*
267* N is odd
268*
269 IF( normaltransr ) THEN
270*
271* N is odd and TRANSR = 'N'
272*
273 IF( lower ) THEN
274*
275* N is odd, TRANSR = 'N', and UPLO = 'L'
276*
277 IF( notrans ) THEN
278*
279* N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
280*
281 CALL ssyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
282 $ beta, c( 1 ), n )
283 CALL ssyrk( 'U', 'N', n2, k, alpha, a( n1+1, 1 ),
284 $ lda,
285 $ beta, c( n+1 ), n )
286 CALL sgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1,
287 $ 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 ssyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
295 $ beta, c( 1 ), n )
296 CALL ssyrk( 'U', 'T', n2, k, alpha, a( 1, n1+1 ),
297 $ lda,
298 $ beta, c( n+1 ), n )
299 CALL sgemm( 'T', 'N', n2, n1, k, alpha, a( 1,
300 $ n1+1 ),
301 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
302*
303 END IF
304*
305 ELSE
306*
307* N is odd, TRANSR = 'N', and UPLO = 'U'
308*
309 IF( notrans ) THEN
310*
311* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
312*
313 CALL ssyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
314 $ beta, c( n2+1 ), n )
315 CALL ssyrk( 'U', 'N', n2, k, alpha, a( n2, 1 ),
316 $ lda,
317 $ beta, c( n1+1 ), n )
318 CALL sgemm( '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 ssyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
326 $ beta, c( n2+1 ), n )
327 CALL ssyrk( 'U', 'T', n2, k, alpha, a( 1, n2 ),
328 $ lda,
329 $ beta, c( n1+1 ), n )
330 CALL sgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
331 $ lda, a( 1, n2 ), lda, beta, c( 1 ), n )
332*
333 END IF
334*
335 END IF
336*
337 ELSE
338*
339* N is odd, and TRANSR = 'T'
340*
341 IF( lower ) THEN
342*
343* N is odd, TRANSR = 'T', and UPLO = 'L'
344*
345 IF( notrans ) THEN
346*
347* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
348*
349 CALL ssyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
350 $ beta, c( 1 ), n1 )
351 CALL ssyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ),
352 $ lda,
353 $ beta, c( 2 ), n1 )
354 CALL sgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
355 $ lda, a( n1+1, 1 ), lda, beta,
356 $ c( n1*n1+1 ), n1 )
357*
358 ELSE
359*
360* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
361*
362 CALL ssyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
363 $ beta, c( 1 ), n1 )
364 CALL ssyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ),
365 $ lda,
366 $ beta, c( 2 ), n1 )
367 CALL sgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
368 $ lda, a( 1, n1+1 ), lda, beta,
369 $ c( n1*n1+1 ), n1 )
370*
371 END IF
372*
373 ELSE
374*
375* N is odd, TRANSR = 'T', and UPLO = 'U'
376*
377 IF( notrans ) THEN
378*
379* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
380*
381 CALL ssyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
382 $ beta, c( n2*n2+1 ), n2 )
383 CALL ssyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ),
384 $ lda,
385 $ beta, c( n1*n2+1 ), n2 )
386 CALL sgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1,
387 $ 1 ),
388 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
389*
390 ELSE
391*
392* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
393*
394 CALL ssyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
395 $ beta, c( n2*n2+1 ), n2 )
396 CALL ssyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ),
397 $ lda,
398 $ beta, c( n1*n2+1 ), n2 )
399 CALL sgemm( 'T', 'N', n2, n1, k, alpha, a( 1,
400 $ n1+1 ),
401 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
402*
403 END IF
404*
405 END IF
406*
407 END IF
408*
409 ELSE
410*
411* N is even
412*
413 IF( normaltransr ) THEN
414*
415* N is even and TRANSR = 'N'
416*
417 IF( lower ) THEN
418*
419* N is even, TRANSR = 'N', and UPLO = 'L'
420*
421 IF( notrans ) THEN
422*
423* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
424*
425 CALL ssyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
426 $ beta, c( 2 ), n+1 )
427 CALL ssyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ),
428 $ lda,
429 $ beta, c( 1 ), n+1 )
430 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1,
431 $ 1 ),
432 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
433 $ n+1 )
434*
435 ELSE
436*
437* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
438*
439 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
440 $ beta, c( 2 ), n+1 )
441 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ),
442 $ lda,
443 $ beta, c( 1 ), n+1 )
444 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1,
445 $ nk+1 ),
446 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
447 $ n+1 )
448*
449 END IF
450*
451 ELSE
452*
453* N is even, TRANSR = 'N', and UPLO = 'U'
454*
455 IF( notrans ) THEN
456*
457* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
458*
459 CALL ssyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
460 $ beta, c( nk+2 ), n+1 )
461 CALL ssyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ),
462 $ lda,
463 $ beta, c( nk+1 ), n+1 )
464 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
465 $ lda, a( nk+1, 1 ), lda, beta, c( 1 ),
466 $ n+1 )
467*
468 ELSE
469*
470* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
471*
472 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
473 $ beta, c( nk+2 ), n+1 )
474 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ),
475 $ lda,
476 $ beta, c( nk+1 ), n+1 )
477 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
478 $ lda, a( 1, nk+1 ), lda, beta, c( 1 ),
479 $ n+1 )
480*
481 END IF
482*
483 END IF
484*
485 ELSE
486*
487* N is even, and TRANSR = 'T'
488*
489 IF( lower ) THEN
490*
491* N is even, TRANSR = 'T', and UPLO = 'L'
492*
493 IF( notrans ) THEN
494*
495* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
496*
497 CALL ssyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
498 $ beta, c( nk+1 ), nk )
499 CALL ssyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ),
500 $ lda,
501 $ beta, c( 1 ), nk )
502 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
503 $ lda, a( nk+1, 1 ), lda, beta,
504 $ c( ( ( nk+1 )*nk )+1 ), nk )
505*
506 ELSE
507*
508* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
509*
510 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
511 $ beta, c( nk+1 ), nk )
512 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ),
513 $ lda,
514 $ beta, c( 1 ), nk )
515 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
516 $ lda, a( 1, nk+1 ), lda, beta,
517 $ c( ( ( nk+1 )*nk )+1 ), nk )
518*
519 END IF
520*
521 ELSE
522*
523* N is even, TRANSR = 'T', and UPLO = 'U'
524*
525 IF( notrans ) THEN
526*
527* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
528*
529 CALL ssyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
530 $ beta, c( nk*( nk+1 )+1 ), nk )
531 CALL ssyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ),
532 $ lda,
533 $ beta, c( nk*nk+1 ), nk )
534 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1,
535 $ 1 ),
536 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
537*
538 ELSE
539*
540* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
541*
542 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
543 $ beta, c( nk*( nk+1 )+1 ), nk )
544 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ),
545 $ lda,
546 $ beta, c( nk*nk+1 ), nk )
547 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1,
548 $ nk+1 ),
549 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
550*
551 END IF
552*
553 END IF
554*
555 END IF
556*
557 END IF
558*
559 RETURN
560*
561* End of SSFRK
562*
563 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
ssfrk
subroutine ssfrk(transr, uplo, trans, n, k, alpha, a, lda, beta, c)
SSFRK performs a symmetric rank-k operation for matrix in RFP format.
Definition ssfrk.f:165