160 SUBROUTINE dsptrf( UPLO, N, AP, IPIV, INFO )
173 DOUBLE PRECISION AP( * )
179 DOUBLE PRECISION ZERO, ONE
180 parameter ( zero = 0.0d+0, one = 1.0d+0 )
181 DOUBLE PRECISION EIGHT, SEVTEN
182 parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
186 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
188 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
189 $ rowmax, t, wk, wkm1, wkp1
194 EXTERNAL lsame, idamax
200 INTRINSIC abs, max, sqrt
207 upper = lsame( uplo,
'U' )
208 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
210 ELSE IF( n.LT.0 )
THEN
214 CALL xerbla(
'DSPTRF', -info )
220 alpha = ( one+sqrt( sevten ) ) / eight
230 kc = ( n-1 )*n / 2 + 1
243 absakk = abs( ap( kc+k-1 ) )
249 imax = idamax( k-1, ap( kc ), 1 )
250 colmax = abs( ap( kc+imax-1 ) )
255 IF( max( absakk, colmax ).EQ.zero )
THEN
263 IF( absakk.GE.alpha*colmax )
THEN
272 kx = imax*( imax+1 ) / 2 + imax
273 DO 20 j = imax + 1, k
274 IF( abs( ap( kx ) ).GT.rowmax )
THEN
275 rowmax = abs( ap( kx ) )
280 kpc = ( imax-1 )*imax / 2 + 1
282 jmax = idamax( imax-1, ap( kpc ), 1 )
283 rowmax = max( rowmax, abs( ap( kpc+jmax-1 ) ) )
286 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
291 ELSE IF( abs( ap( kpc+imax-1 ) ).GE.alpha*rowmax )
THEN
315 CALL dswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
317 DO 30 j = kp + 1, kk - 1
320 ap( knc+j-1 ) = ap( kx )
324 ap( knc+kk-1 ) = ap( kpc+kp-1 )
326 IF( kstep.EQ.2 )
THEN
328 ap( kc+k-2 ) = ap( kc+kp-1 )
335 IF( kstep.EQ.1 )
THEN
347 r1 = one / ap( kc+k-1 )
348 CALL dspr( uplo, k-1, -r1, ap( kc ), 1, ap )
352 CALL dscal( k-1, r1, ap( kc ), 1 )
369 d12 = ap( k-1+( k-1 )*k / 2 )
370 d22 = ap( k-1+( k-2 )*( k-1 ) / 2 ) / d12
371 d11 = ap( k+( k-1 )*k / 2 ) / d12
372 t = one / ( d11*d22-one )
375 DO 50 j = k - 2, 1, -1
376 wkm1 = d12*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
377 $ ap( j+( k-1 )*k / 2 ) )
378 wk = d12*( d22*ap( j+( k-1 )*k / 2 )-
379 $ ap( j+( k-2 )*( k-1 ) / 2 ) )
381 ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
382 $ ap( i+( k-1 )*k / 2 )*wk -
383 $ ap( i+( k-2 )*( k-1 ) / 2 )*wkm1
385 ap( j+( k-1 )*k / 2 ) = wk
386 ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
396 IF( kstep.EQ.1 )
THEN
431 absakk = abs( ap( kc ) )
437 imax = k + idamax( n-k, ap( kc+1 ), 1 )
438 colmax = abs( ap( kc+imax-k ) )
443 IF( max( absakk, colmax ).EQ.zero )
THEN
451 IF( absakk.GE.alpha*colmax )
THEN
463 DO 70 j = k, imax - 1
464 IF( abs( ap( kx ) ).GT.rowmax )
THEN
465 rowmax = abs( ap( kx ) )
470 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
472 jmax = imax + idamax( n-imax, ap( kpc+1 ), 1 )
473 rowmax = max( rowmax, abs( ap( kpc+jmax-imax ) ) )
476 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
481 ELSE IF( abs( ap( kpc ) ).GE.alpha*rowmax )
THEN
499 $ knc = knc + n - k + 1
506 $
CALL dswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
509 DO 80 j = kk + 1, kp - 1
512 ap( knc+j-kk ) = ap( kx )
516 ap( knc ) = ap( kpc )
518 IF( kstep.EQ.2 )
THEN
520 ap( kc+1 ) = ap( kc+kp-k )
527 IF( kstep.EQ.1 )
THEN
542 CALL dspr( uplo, n-k, -r1, ap( kc+1 ), 1,
547 CALL dscal( n-k, r1, ap( kc+1 ), 1 )
568 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 )
569 d11 = ap( k+1+k*( 2*n-k-1 ) / 2 ) / d21
570 d22 = ap( k+( k-1 )*( 2*n-k ) / 2 ) / d21
571 t = one / ( d11*d22-one )
575 wk = d21*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-
576 $ ap( j+k*( 2*n-k-1 ) / 2 ) )
577 wkp1 = d21*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
578 $ ap( j+( k-1 )*( 2*n-k ) / 2 ) )
581 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
582 $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
583 $ 2 )*wk - ap( i+k*( 2*n-k-1 ) / 2 )*wkp1
586 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
587 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
596 IF( kstep.EQ.1 )
THEN
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dspr(UPLO, N, ALPHA, X, INCX, AP)
DSPR
subroutine dsptrf(UPLO, N, AP, IPIV, INFO)
DSPTRF