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