LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dsfrk()

subroutine dsfrk ( character transr,
character uplo,
character trans,
integer n,
integer k,
double precision alpha,
double precision, dimension( lda, * ) a,
integer lda,
double precision beta,
double precision, dimension( * ) c )

DSFRK performs a symmetric rank-k operation for matrix in RFP format.

Download DSFRK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> Level 3 BLAS like routine for C in RFP Format.
!>
!> DSFRK performs one of the symmetric rank--k operations
!>
!>    C := alpha*A*A**T + beta*C,
!>
!> or
!>
!>    C := alpha*A**T*A + beta*C,
!>
!> where alpha and beta are real scalars, C is an n--by--n symmetric
!> matrix and A is an n--by--k matrix in the first case and a k--by--n
!> matrix in the second case.
!>
Parameters
[in]TRANSR
!>          TRANSR is CHARACTER*1
!>          = 'N':  The Normal Form of RFP A is stored;
!>          = 'T':  The Transpose Form of RFP A is stored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>           On  entry, UPLO specifies whether the upper or lower
!>           triangular part of the array C is to be referenced as
!>           follows:
!>
!>              UPLO = 'U' or 'u'   Only the upper triangular part of C
!>                                  is to be referenced.
!>
!>              UPLO = 'L' or 'l'   Only the lower triangular part of C
!>                                  is to be referenced.
!>
!>           Unchanged on exit.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>           On entry, TRANS specifies the operation to be performed as
!>           follows:
!>
!>              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
!>
!>              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
!>
!>           Unchanged on exit.
!>
[in]N
!>          N is INTEGER
!>           On entry, N specifies the order of the matrix C. N must be
!>           at least zero.
!>           Unchanged on exit.
!>
[in]K
!>          K is INTEGER
!>           On entry with TRANS = 'N' or 'n', K specifies the number
!>           of  columns of the matrix A, and on entry with TRANS = 'T'
!>           or 't', K specifies the number of rows of the matrix A. K
!>           must be at least zero.
!>           Unchanged on exit.
!>
[in]ALPHA
!>          ALPHA is DOUBLE PRECISION
!>           On entry, ALPHA specifies the scalar alpha.
!>           Unchanged on exit.
!>
[in]A
!>          A is DOUBLE PRECISION array, dimension (LDA,ka)
!>           where KA
!>           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
!>           entry with TRANS = 'N' or 'n', the leading N--by--K part of
!>           the array A must contain the matrix A, otherwise the leading
!>           K--by--N part of the array A must contain the matrix A.
!>           Unchanged on exit.
!>
[in]LDA
!>          LDA is INTEGER
!>           On entry, LDA specifies the first dimension of A as declared
!>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
!>           be at least  max( 1, k ).
!>           Unchanged on exit.
!>
[in]BETA
!>          BETA is DOUBLE PRECISION
!>           On entry, BETA specifies the scalar beta.
!>           Unchanged on exit.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (NT)
!>           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
!>           Format. RFP Format is described by TRANSR, UPLO and N.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 162 of file dsfrk.f.

165*
166* -- LAPACK computational routine --
167* -- LAPACK is a software package provided by Univ. of Tennessee, --
168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*
170* .. Scalar Arguments ..
171 DOUBLE PRECISION ALPHA, BETA
172 INTEGER K, LDA, N
173 CHARACTER TRANS, TRANSR, UPLO
174* ..
175* .. Array Arguments ..
176 DOUBLE PRECISION A( LDA, * ), C( * )
177* ..
178*
179* =====================================================================
180*
181* ..
182* .. Parameters ..
183 DOUBLE PRECISION ONE, ZERO
184 parameter( one = 1.0d+0, zero = 0.0d+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 xerbla, dgemm, dsyrk
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( 'DSFRK ', -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 DSYRK 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 dsyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
283 $ beta, c( 1 ), n )
284 CALL dsyrk( 'U', 'N', n2, k, alpha, a( n1+1, 1 ),
285 $ lda,
286 $ beta, c( n+1 ), n )
287 CALL dgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1,
288 $ 1 ),
289 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
290*
291 ELSE
292*
293* N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
294*
295 CALL dsyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
296 $ beta, c( 1 ), n )
297 CALL dsyrk( 'U', 'T', n2, k, alpha, a( 1, n1+1 ),
298 $ lda,
299 $ beta, c( n+1 ), n )
300 CALL dgemm( 'T', 'N', n2, n1, k, alpha, a( 1,
301 $ 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 ),
317 $ lda,
318 $ beta, c( n1+1 ), n )
319 CALL dgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
320 $ lda, a( n2, 1 ), lda, beta, c( 1 ), n )
321*
322 ELSE
323*
324* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
325*
326 CALL dsyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
327 $ beta, c( n2+1 ), n )
328 CALL dsyrk( 'U', 'T', n2, k, alpha, a( 1, n2 ),
329 $ lda,
330 $ beta, c( n1+1 ), n )
331 CALL dgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
332 $ lda, a( 1, n2 ), lda, beta, c( 1 ), n )
333*
334 END IF
335*
336 END IF
337*
338 ELSE
339*
340* N is odd, and TRANSR = 'T'
341*
342 IF( lower ) THEN
343*
344* N is odd, TRANSR = 'T', and UPLO = 'L'
345*
346 IF( notrans ) THEN
347*
348* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
349*
350 CALL dsyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
351 $ beta, c( 1 ), n1 )
352 CALL dsyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ),
353 $ lda,
354 $ beta, c( 2 ), n1 )
355 CALL dgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
356 $ lda, a( n1+1, 1 ), lda, beta,
357 $ c( n1*n1+1 ), n1 )
358*
359 ELSE
360*
361* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
362*
363 CALL dsyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
364 $ beta, c( 1 ), n1 )
365 CALL dsyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ),
366 $ lda,
367 $ beta, c( 2 ), n1 )
368 CALL dgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
369 $ lda, a( 1, n1+1 ), lda, beta,
370 $ c( n1*n1+1 ), n1 )
371*
372 END IF
373*
374 ELSE
375*
376* N is odd, TRANSR = 'T', and UPLO = 'U'
377*
378 IF( notrans ) THEN
379*
380* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
381*
382 CALL dsyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
383 $ beta, c( n2*n2+1 ), n2 )
384 CALL dsyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ),
385 $ lda,
386 $ beta, c( n1*n2+1 ), n2 )
387 CALL dgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1,
388 $ 1 ),
389 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
390*
391 ELSE
392*
393* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
394*
395 CALL dsyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
396 $ beta, c( n2*n2+1 ), n2 )
397 CALL dsyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ),
398 $ lda,
399 $ beta, c( n1*n2+1 ), n2 )
400 CALL dgemm( 'T', 'N', n2, n1, k, alpha, a( 1,
401 $ n1+1 ),
402 $ lda, a( 1, 1 ), lda, beta, c( 1 ), n2 )
403*
404 END IF
405*
406 END IF
407*
408 END IF
409*
410 ELSE
411*
412* N is even
413*
414 IF( normaltransr ) THEN
415*
416* N is even and TRANSR = 'N'
417*
418 IF( lower ) THEN
419*
420* N is even, TRANSR = 'N', and UPLO = 'L'
421*
422 IF( notrans ) THEN
423*
424* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
425*
426 CALL dsyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
427 $ beta, c( 2 ), n+1 )
428 CALL dsyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ),
429 $ lda,
430 $ beta, c( 1 ), n+1 )
431 CALL dgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1,
432 $ 1 ),
433 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
434 $ n+1 )
435*
436 ELSE
437*
438* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
439*
440 CALL dsyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
441 $ beta, c( 2 ), n+1 )
442 CALL dsyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ),
443 $ lda,
444 $ beta, c( 1 ), n+1 )
445 CALL dgemm( 'T', 'N', nk, nk, k, alpha, a( 1,
446 $ nk+1 ),
447 $ lda, a( 1, 1 ), lda, beta, c( nk+2 ),
448 $ n+1 )
449*
450 END IF
451*
452 ELSE
453*
454* N is even, TRANSR = 'N', and UPLO = 'U'
455*
456 IF( notrans ) THEN
457*
458* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
459*
460 CALL dsyrk( 'L', 'N', nk, k, alpha, a( 1, 1 ), lda,
461 $ beta, c( nk+2 ), n+1 )
462 CALL dsyrk( 'U', 'N', nk, k, alpha, a( nk+1, 1 ),
463 $ lda,
464 $ beta, c( nk+1 ), n+1 )
465 CALL dgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
466 $ lda, a( nk+1, 1 ), lda, beta, c( 1 ),
467 $ n+1 )
468*
469 ELSE
470*
471* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
472*
473 CALL dsyrk( 'L', 'T', nk, k, alpha, a( 1, 1 ), lda,
474 $ beta, c( nk+2 ), n+1 )
475 CALL dsyrk( 'U', 'T', nk, k, alpha, a( 1, nk+1 ),
476 $ lda,
477 $ beta, c( nk+1 ), n+1 )
478 CALL dgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
479 $ lda, a( 1, nk+1 ), lda, beta, c( 1 ),
480 $ n+1 )
481*
482 END IF
483*
484 END IF
485*
486 ELSE
487*
488* N is even, and TRANSR = 'T'
489*
490 IF( lower ) THEN
491*
492* N is even, TRANSR = 'T', and UPLO = 'L'
493*
494 IF( notrans ) THEN
495*
496* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
497*
498 CALL dsyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
499 $ beta, c( nk+1 ), nk )
500 CALL dsyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ),
501 $ lda,
502 $ beta, c( 1 ), nk )
503 CALL dgemm( 'N', 'T', nk, nk, k, alpha, a( 1, 1 ),
504 $ lda, a( nk+1, 1 ), lda, beta,
505 $ c( ( ( nk+1 )*nk )+1 ), nk )
506*
507 ELSE
508*
509* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
510*
511 CALL dsyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
512 $ beta, c( nk+1 ), nk )
513 CALL dsyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ),
514 $ lda,
515 $ beta, c( 1 ), nk )
516 CALL dgemm( 'T', 'N', nk, nk, k, alpha, a( 1, 1 ),
517 $ lda, a( 1, nk+1 ), lda, beta,
518 $ c( ( ( nk+1 )*nk )+1 ), nk )
519*
520 END IF
521*
522 ELSE
523*
524* N is even, TRANSR = 'T', and UPLO = 'U'
525*
526 IF( notrans ) THEN
527*
528* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
529*
530 CALL dsyrk( 'U', 'N', nk, k, alpha, a( 1, 1 ), lda,
531 $ beta, c( nk*( nk+1 )+1 ), nk )
532 CALL dsyrk( 'L', 'N', nk, k, alpha, a( nk+1, 1 ),
533 $ lda,
534 $ beta, c( nk*nk+1 ), nk )
535 CALL dgemm( 'N', 'T', nk, nk, k, alpha, a( nk+1,
536 $ 1 ),
537 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
538*
539 ELSE
540*
541* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
542*
543 CALL dsyrk( 'U', 'T', nk, k, alpha, a( 1, 1 ), lda,
544 $ beta, c( nk*( nk+1 )+1 ), nk )
545 CALL dsyrk( 'L', 'T', nk, k, alpha, a( 1, nk+1 ),
546 $ lda,
547 $ beta, c( nk*nk+1 ), nk )
548 CALL dgemm( 'T', 'N', nk, nk, k, alpha, a( 1,
549 $ nk+1 ),
550 $ lda, a( 1, 1 ), lda, beta, c( 1 ), nk )
551*
552 END IF
553*
554 END IF
555*
556 END IF
557*
558 END IF
559*
560 RETURN
561*
562* End of DSFRK
563*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dsyrk
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: