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