LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 realOTHERcomputational
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187