159 $ IPIV2, WORK, LWORK, INFO )
169 INTEGER N, LDA, LTB, LWORK, INFO
172 INTEGER IPIV( * ), IPIV2( * )
173 DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
178 DOUBLE PRECISION ZERO, ONE
179 parameter( zero = 0.0d+0, one = 1.0d+0 )
182 LOGICAL UPPER, TQUERY, WQUERY
183 INTEGER I, J, K, I1, I2, TD
184 INTEGER LDTB, NB, KB, JB, NT, IINFO
190 EXTERNAL lsame, ilaenv
205 upper = lsame( uplo,
'U' )
206 wquery = ( lwork.EQ.-1 )
207 tquery = ( ltb.EQ.-1 )
208 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
210 ELSE IF( n.LT.0 )
THEN
212 ELSE IF( lda.LT.max( 1, n ) )
THEN
214 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery )
THEN
216 ELSE IF ( lwork .LT. n .AND. .NOT.wquery )
THEN
221 CALL xerbla(
'DSYTRF_AA_2STAGE', -info )
227 nb = ilaenv( 1,
'DSYTRF_AA_2STAGE', uplo, n, -1, -1, -1 )
236 IF( tquery .OR. wquery )
THEN
249 IF( ldtb .LT. 3*nb+1 )
THEN
252 IF( lwork .LT. nb*n )
THEN
286 IF( i .EQ. (j-1) )
THEN
291 CALL dgemm(
'NoTranspose',
'NoTranspose',
293 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
294 $ a( (i-1)*nb+1, j*nb+1 ), lda,
295 $ zero, work( i*nb+1 ), n )
303 CALL dgemm(
'NoTranspose',
'NoTranspose',
305 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
307 $ a( (i-2)*nb+1, j*nb+1 ), lda,
308 $ zero, work( i*nb+1 ), n )
314 CALL dlacpy(
'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
315 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
318 CALL dgemm(
'Transpose',
'NoTranspose',
320 $ -one, a( 1, j*nb+1 ), lda,
322 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
324 CALL dgemm(
'Transpose',
'NoTranspose',
326 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
327 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
328 $ zero, work( 1 ), n )
329 CALL dgemm(
'NoTranspose',
'NoTranspose',
331 $ -one, work( 1 ), n,
332 $ a( (j-2)*nb+1, j*nb+1 ), lda,
333 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
336 CALL dsygst( 1,
'Upper', kb,
337 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
338 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
345 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
346 $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
356 CALL dgemm(
'NoTranspose',
'NoTranspose',
358 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
359 $ a( (j-1)*nb+1, j*nb+1 ), lda,
360 $ zero, work( j*nb+1 ), n )
362 CALL dgemm(
'NoTranspose',
'NoTranspose',
364 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
366 $ a( (j-2)*nb+1, j*nb+1 ), lda,
367 $ zero, work( j*nb+1 ), n )
372 CALL dgemm(
'Transpose',
'NoTranspose',
373 $ nb, n-(j+1)*nb, j*nb,
374 $ -one, work( nb+1 ), n,
375 $ a( 1, (j+1)*nb+1 ), lda,
376 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
382 CALL dcopy( n-(j+1)*nb,
383 $ a( j*nb+k, (j+1)*nb+1 ), lda,
384 $ work( 1+(k-1)*n ), 1 )
389 CALL dgetrf( n-(j+1)*nb, nb,
391 $ ipiv( (j+1)*nb+1 ), iinfo )
399 CALL dcopy( n-(j+1)*nb,
400 $ work( 1+(k-1)*n ), 1,
401 $ a( j*nb+k, (j+1)*nb+1 ), lda )
406 kb = min(nb, n-(j+1)*nb)
407 CALL dlaset(
'Full', kb, nb, zero, zero,
408 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
409 CALL dlacpy(
'Upper', kb, nb,
411 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
413 CALL dtrsm(
'R',
'U',
'N',
'U', kb, nb, one,
414 $ a( (j-1)*nb+1, j*nb+1 ), lda,
415 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
423 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
424 $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
427 CALL dlaset(
'Lower', kb, nb, zero, one,
428 $ a( j*nb+1, (j+1)*nb+1), lda )
434 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
437 i2 = ipiv( (j+1)*nb+k )
440 CALL dswap( k-1, a( (j+1)*nb+1, i1 ), 1,
441 $ a( (j+1)*nb+1, i2 ), 1 )
444 $
CALL dswap( i2-i1-1, a( i1, i1+1 ), lda,
448 $
CALL dswap( n-i2, a( i1, i2+1 ), lda,
449 $ a( i2, i2+1 ), lda )
452 a( i1, i1 ) = a( i2, i2 )
456 CALL dswap( j*nb, a( 1, i1 ), 1,
482 CALL dgemm(
'NoTranspose',
'Transpose',
484 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
485 $ a( j*nb+1, (i-1)*nb+1 ), lda,
486 $ zero, work( i*nb+1 ), n )
494 CALL dgemm(
'NoTranspose',
'Transpose',
496 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
498 $ a( j*nb+1, (i-2)*nb+1 ), lda,
499 $ zero, work( i*nb+1 ), n )
505 CALL dlacpy(
'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
506 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
509 CALL dgemm(
'NoTranspose',
'NoTranspose',
511 $ -one, a( j*nb+1, 1 ), lda,
513 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
515 CALL dgemm(
'NoTranspose',
'NoTranspose',
517 $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
518 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
519 $ zero, work( 1 ), n )
520 CALL dgemm(
'NoTranspose',
'Transpose',
522 $ -one, work( 1 ), n,
523 $ a( j*nb+1, (j-2)*nb+1 ), lda,
524 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
527 CALL dsygst( 1,
'Lower', kb,
528 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
529 $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
536 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
537 $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
547 CALL dgemm(
'NoTranspose',
'Transpose',
549 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
550 $ a( j*nb+1, (j-1)*nb+1 ), lda,
551 $ zero, work( j*nb+1 ), n )
553 CALL dgemm(
'NoTranspose',
'Transpose',
555 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
557 $ a( j*nb+1, (j-2)*nb+1 ), lda,
558 $ zero, work( j*nb+1 ), n )
563 CALL dgemm(
'NoTranspose',
'NoTranspose',
564 $ n-(j+1)*nb, nb, j*nb,
565 $ -one, a( (j+1)*nb+1, 1 ), lda,
567 $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
572 CALL dgetrf( n-(j+1)*nb, nb,
573 $ a( (j+1)*nb+1, j*nb+1 ), lda,
574 $ ipiv( (j+1)*nb+1 ), iinfo )
581 kb = min(nb, n-(j+1)*nb)
582 CALL dlaset(
'Full', kb, nb, zero, zero,
583 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
584 CALL dlacpy(
'Upper', kb, nb,
585 $ a( (j+1)*nb+1, j*nb+1 ), lda,
586 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
588 CALL dtrsm(
'R',
'L',
'T',
'U', kb, nb, one,
589 $ a( j*nb+1, (j-1)*nb+1 ), lda,
590 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
598 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
599 $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
602 CALL dlaset(
'Upper', kb, nb, zero, one,
603 $ a( (j+1)*nb+1, j*nb+1), lda )
609 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
612 i2 = ipiv( (j+1)*nb+k )
615 CALL dswap( k-1, a( i1, (j+1)*nb+1 ), lda,
616 $ a( i2, (j+1)*nb+1 ), lda )
619 $
CALL dswap( i2-i1-1, a( i1+1, i1 ), 1,
620 $ a( i2, i1+1 ), lda )
623 $
CALL dswap( n-i2, a( i2+1, i1 ), 1,
627 a( i1, i1 ) = a( i2, i2 )
631 CALL dswap( j*nb, a( i1, 1 ), lda,
646 CALL dgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
subroutine xerbla(srname, info)
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
subroutine dgbtrf(m, n, kl, ku, ab, ldab, ipiv, info)
DGBTRF
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dgetrf(m, n, a, lda, ipiv, info)
DGETRF
subroutine dsygst(itype, uplo, n, a, lda, b, ldb, info)
DSYGST
subroutine dsytrf_aa_2stage(uplo, n, a, lda, tb, ltb, ipiv, ipiv2, work, lwork, info)
DSYTRF_AA_2STAGE
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM