LAPACK 3.12.0
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*> \htmlonly
9*> Download SSFRK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssfrk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssfrk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssfrk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
22* C )
23*
24* .. Scalar Arguments ..
25* REAL ALPHA, BETA
26* INTEGER K, LDA, N
27* CHARACTER TRANS, TRANSR, UPLO
28* ..
29* .. Array Arguments ..
30* REAL 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*> SSFRK 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 REAL
114*> On entry, ALPHA specifies the scalar alpha.
115*> Unchanged on exit.
116*> \endverbatim
117*>
118*> \param[in] A
119*> \verbatim
120*> A is REAL 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 REAL
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 REAL 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 hfrk
162*
163* =====================================================================
164 SUBROUTINE ssfrk( 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 REAL ALPHA, BETA
173 INTEGER K, LDA, N
174 CHARACTER TRANS, TRANSR, UPLO
175* ..
176* .. Array Arguments ..
177 REAL A( LDA, * ), C( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 REAL ONE, ZERO
184 parameter( one = 1.0e+0, zero = 0.0e+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 sgemm, ssyrk, xerbla
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( 'SSFRK ', -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 SSYRK 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 ssyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
283 $ beta, c( 1 ), n )
284 CALL ssyrk( 'U', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
285 $ beta, c( n+1 ), n )
286 CALL sgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1, 1 ),
287 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
288*
289 ELSE
290*
291* N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
292*
293 CALL ssyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
294 $ beta, c( 1 ), n )
295 CALL ssyrk( 'U', 'T', n2, k, alpha, a( 1, n1+1 ), lda,
296 $ beta, c( n+1 ), n )
297 CALL sgemm( 'T', 'N', n2, n1, k, alpha, a( 1, n1+1 ),
298 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
299*
300 END IF
301*
302 ELSE
303*
304* N is odd, TRANSR = 'N', and UPLO = 'U'
305*
306 IF( notrans ) THEN
307*
308* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
309*
310 CALL ssyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
311 $ beta, c( n2+1 ), n )
312 CALL ssyrk( 'U', 'N', n2, k, alpha, a( n2, 1 ), lda,
313 $ beta, c( n1+1 ), n )
314 CALL sgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
315 $ lda, a( n2, 1 ), lda, beta, c( 1 ), n )
316*
317 ELSE
318*
319* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
320*
321 CALL ssyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
322 $ beta, c( n2+1 ), n )
323 CALL ssyrk( 'U', 'T', n2, k, alpha, a( 1, n2 ), lda,
324 $ beta, c( n1+1 ), n )
325 CALL sgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
326 $ lda, a( 1, n2 ), lda, beta, c( 1 ), n )
327*
328 END IF
329*
330 END IF
331*
332 ELSE
333*
334* N is odd, and TRANSR = 'T'
335*
336 IF( lower ) THEN
337*
338* N is odd, TRANSR = 'T', and UPLO = 'L'
339*
340 IF( notrans ) THEN
341*
342* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
343*
344 CALL ssyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
345 $ beta, c( 1 ), n1 )
346 CALL ssyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
347 $ beta, c( 2 ), n1 )
348 CALL sgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
349 $ lda, a( n1+1, 1 ), lda, beta,
350 $ c( n1*n1+1 ), n1 )
351*
352 ELSE
353*
354* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
355*
356 CALL ssyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
357 $ beta, c( 1 ), n1 )
358 CALL ssyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ), lda,
359 $ beta, c( 2 ), n1 )
360 CALL sgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
361 $ lda, a( 1, n1+1 ), lda, beta,
362 $ c( n1*n1+1 ), n1 )
363*
364 END IF
365*
366 ELSE
367*
368* N is odd, TRANSR = 'T', and UPLO = 'U'
369*
370 IF( notrans ) THEN
371*
372* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
373*
374 CALL ssyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
375 $ beta, c( n2*n2+1 ), n2 )
376 CALL ssyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
377 $ beta, c( n1*n2+1 ), n2 )
378 CALL sgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1, 1 ),
379 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
380*
381 ELSE
382*
383* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
384*
385 CALL ssyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
386 $ beta, c( n2*n2+1 ), n2 )
387 CALL ssyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ), lda,
388 $ beta, c( n1*n2+1 ), n2 )
389 CALL sgemm( 'T', 'N', n2, n1, k, alpha, a( 1, n1+1 ),
390 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
391*
392 END IF
393*
394 END IF
395*
396 END IF
397*
398 ELSE
399*
400* N is even
401*
402 IF( normaltransr ) THEN
403*
404* N is even and TRANSR = 'N'
405*
406 IF( lower ) THEN
407*
408* N is even, TRANSR = 'N', and UPLO = 'L'
409*
410 IF( notrans ) THEN
411*
412* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
413*
414 CALL ssyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
415 $ beta, c( 2 ), n+1 )
416 CALL ssyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
417 $ beta, c( 1 ), n+1 )
418 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1, 1 ),
419 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
420 $ n+1 )
421*
422 ELSE
423*
424* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
425*
426 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
427 $ beta, c( 2 ), n+1 )
428 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ), lda,
429 $ beta, c( 1 ), n+1 )
430 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, nk+1 ),
431 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
432 $ n+1 )
433*
434 END IF
435*
436 ELSE
437*
438* N is even, TRANSR = 'N', and UPLO = 'U'
439*
440 IF( notrans ) THEN
441*
442* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
443*
444 CALL ssyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
445 $ beta, c( nk+2 ), n+1 )
446 CALL ssyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
447 $ beta, c( nk+1 ), n+1 )
448 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
449 $ lda, a( nk+1, 1 ), lda, beta, c( 1 ),
450 $ n+1 )
451*
452 ELSE
453*
454* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
455*
456 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
457 $ beta, c( nk+2 ), n+1 )
458 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ), lda,
459 $ beta, c( nk+1 ), n+1 )
460 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
461 $ lda, a( 1, nk+1 ), lda, beta, c( 1 ),
462 $ n+1 )
463*
464 END IF
465*
466 END IF
467*
468 ELSE
469*
470* N is even, and TRANSR = 'T'
471*
472 IF( lower ) THEN
473*
474* N is even, TRANSR = 'T', and UPLO = 'L'
475*
476 IF( notrans ) THEN
477*
478* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
479*
480 CALL ssyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
481 $ beta, c( nk+1 ), nk )
482 CALL ssyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
483 $ beta, c( 1 ), nk )
484 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
485 $ lda, a( nk+1, 1 ), lda, beta,
486 $ c( ( ( nk+1 )*nk )+1 ), nk )
487*
488 ELSE
489*
490* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
491*
492 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
493 $ beta, c( nk+1 ), nk )
494 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ), lda,
495 $ beta, c( 1 ), nk )
496 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
497 $ lda, a( 1, nk+1 ), lda, beta,
498 $ c( ( ( nk+1 )*nk )+1 ), nk )
499*
500 END IF
501*
502 ELSE
503*
504* N is even, TRANSR = 'T', and UPLO = 'U'
505*
506 IF( notrans ) THEN
507*
508* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
509*
510 CALL ssyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
511 $ beta, c( nk*( nk+1 )+1 ), nk )
512 CALL ssyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ), lda,
513 $ beta, c( nk*nk+1 ), nk )
514 CALL sgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1, 1 ),
515 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
516*
517 ELSE
518*
519* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
520*
521 CALL ssyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
522 $ beta, c( nk*( nk+1 )+1 ), nk )
523 CALL ssyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ), lda,
524 $ beta, c( nk*nk+1 ), nk )
525 CALL sgemm( 'T', 'N', nk, nk, k, alpha, a( 1, nk+1 ),
526 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
527*
528 END IF
529*
530 END IF
531*
532 END IF
533*
534 END IF
535*
536 RETURN
537*
538* End of SSFRK
539*
540 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:166