158 SUBROUTINE dsptrf( UPLO, N, AP, IPIV, INFO )
170 DOUBLE PRECISION AP( * )
176 DOUBLE PRECISION ZERO, ONE
177 parameter( zero = 0.0d+0, one = 1.0d+0 )
178 DOUBLE PRECISION EIGHT, SEVTEN
179 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
183 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
185 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
186 $ ROWMAX, T, WK, WKM1, WKP1
191 EXTERNAL lsame, idamax
197 INTRINSIC abs, max, sqrt
204 upper = lsame( uplo,
'U' )
205 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
207 ELSE IF( n.LT.0 )
THEN
211 CALL xerbla(
'DSPTRF', -info )
217 alpha = ( one+sqrt( sevten ) ) / eight
227 kc = ( n-1 )*n / 2 + 1
240 absakk = abs( ap( kc+k-1 ) )
246 imax = idamax( k-1, ap( kc ), 1 )
247 colmax = abs( ap( kc+imax-1 ) )
252 IF( max( absakk, colmax ).EQ.zero )
THEN
260 IF( absakk.GE.alpha*colmax )
THEN
269 kx = imax*( imax+1 ) / 2 + imax
270 DO 20 j = imax + 1, k
271 IF( abs( ap( kx ) ).GT.rowmax )
THEN
272 rowmax = abs( ap( kx ) )
277 kpc = ( imax-1 )*imax / 2 + 1
279 jmax = idamax( imax-1, ap( kpc ), 1 )
280 rowmax = max( rowmax, abs( ap( kpc+jmax-1 ) ) )
283 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
288 ELSE IF( abs( ap( kpc+imax-1 ) ).GE.alpha*rowmax )
THEN
312 CALL dswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
314 DO 30 j = kp + 1, kk - 1
317 ap( knc+j-1 ) = ap( kx )
321 ap( knc+kk-1 ) = ap( kpc+kp-1 )
323 IF( kstep.EQ.2 )
THEN
325 ap( kc+k-2 ) = ap( kc+kp-1 )
332 IF( kstep.EQ.1 )
THEN
344 r1 = one / ap( kc+k-1 )
345 CALL dspr( uplo, k-1, -r1, ap( kc ), 1, ap )
349 CALL dscal( k-1, r1, ap( kc ), 1 )
366 d12 = ap( k-1+( k-1 )*k / 2 )
367 d22 = ap( k-1+( k-2 )*( k-1 ) / 2 ) / d12
368 d11 = ap( k+( k-1 )*k / 2 ) / d12
369 t = one / ( d11*d22-one )
372 DO 50 j = k - 2, 1, -1
373 wkm1 = d12*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
374 $ ap( j+( k-1 )*k / 2 ) )
375 wk = d12*( d22*ap( j+( k-1 )*k / 2 )-
376 $ ap( j+( k-2 )*( k-1 ) / 2 ) )
378 ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
379 $ ap( i+( k-1 )*k / 2 )*wk -
380 $ ap( i+( k-2 )*( k-1 ) / 2 )*wkm1
382 ap( j+( k-1 )*k / 2 ) = wk
383 ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
393 IF( kstep.EQ.1 )
THEN
428 absakk = abs( ap( kc ) )
434 imax = k + idamax( n-k, ap( kc+1 ), 1 )
435 colmax = abs( ap( kc+imax-k ) )
440 IF( max( absakk, colmax ).EQ.zero )
THEN
448 IF( absakk.GE.alpha*colmax )
THEN
460 DO 70 j = k, imax - 1
461 IF( abs( ap( kx ) ).GT.rowmax )
THEN
462 rowmax = abs( ap( kx ) )
467 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
469 jmax = imax + idamax( n-imax, ap( kpc+1 ), 1 )
470 rowmax = max( rowmax, abs( ap( kpc+jmax-imax ) ) )
473 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
478 ELSE IF( abs( ap( kpc ) ).GE.alpha*rowmax )
THEN
496 $ knc = knc + n - k + 1
503 $
CALL dswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
506 DO 80 j = kk + 1, kp - 1
509 ap( knc+j-kk ) = ap( kx )
513 ap( knc ) = ap( kpc )
515 IF( kstep.EQ.2 )
THEN
517 ap( kc+1 ) = ap( kc+kp-k )
524 IF( kstep.EQ.1 )
THEN
539 CALL dspr( uplo, n-k, -r1, ap( kc+1 ), 1,
544 CALL dscal( n-k, r1, ap( kc+1 ), 1 )
565 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 )
566 d11 = ap( k+1+k*( 2*n-k-1 ) / 2 ) / d21
567 d22 = ap( k+( k-1 )*( 2*n-k ) / 2 ) / d21
568 t = one / ( d11*d22-one )
572 wk = d21*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-
573 $ ap( j+k*( 2*n-k-1 ) / 2 ) )
574 wkp1 = d21*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
575 $ ap( j+( k-1 )*( 2*n-k ) / 2 ) )
578 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
579 $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
580 $ 2 )*wk - ap( i+k*( 2*n-k-1 ) / 2 )*wkp1
583 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
584 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
593 IF( kstep.EQ.1 )
THEN
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dspr(UPLO, N, ALPHA, X, INCX, AP)
DSPR
subroutine dsptrf(UPLO, N, AP, IPIV, INFO)
DSPTRF