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