SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL AP( * ) * .. * * Purpose * ======= * * SSPTRF computes the factorization of a real symmetric matrix A stored * in packed format using the Bunch-Kaufman diagonal pivoting method: * * A = U*D*U**T or A = L*D*L**T * * where U (or L) is a product of permutation and unit upper (lower) * triangular matrices, and D is symmetric and block diagonal with * 1-by-1 and 2-by-2 diagonal blocks. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * = 'U': Upper triangle of A is stored; * = 'L': Lower triangle of A is stored. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * AP (input/output) REAL array, dimension (N*(N+1)/2) * On entry, the upper or lower triangle of the symmetric matrix * A, packed columnwise in a linear array. The j-th column of A * is stored in the array AP as follows: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. * * On exit, the block diagonal matrix D and the multipliers used * to obtain the factor U or L, stored as a packed triangular * matrix overwriting A (see below for further details). * * IPIV (output) INTEGER array, dimension (N) * Details of the interchanges and the block structure of D. * If IPIV(k) > 0, then rows and columns k and IPIV(k) were * interchanged and D(k,k) is a 1-by-1 diagonal block. * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, D(i,i) is exactly zero. The factorization * has been completed, but the block diagonal matrix D is * exactly singular, and division by zero will occur if it * is used to solve a system of equations. * * Further Details * =============== * * 5-96 - Based on modifications by J. Lewis, Boeing Computer Services * Company * * If UPLO = 'U', then A = U*D*U', where * U = P(n)*U(n)* ... *P(k)U(k)* ..., * i.e., U is a product of terms P(k)*U(k), where k decreases from n to * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such * that if the diagonal block D(k) is of order s (s = 1 or 2), then * * ( I v 0 ) k-s * U(k) = ( 0 I 0 ) s * ( 0 0 I ) n-k * k-s s n-k * * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), * and A(k,k), and v overwrites A(1:k-2,k-1:k). * * If UPLO = 'L', then A = L*D*L', where * L = P(1)*L(1)* ... *P(k)*L(k)* ..., * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such * that if the diagonal block D(k) is of order s (s = 1 or 2), then * * ( I 0 0 ) k-1 * L(k) = ( 0 I 0 ) s * ( 0 v I ) n-k-s+1 * k-1 s n-k-s+1 * * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) REAL EIGHT, SEVTEN PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC, $ KSTEP, KX, NPP REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1, $ ROWMAX, T, WK, WKM1, WKP1 * .. * .. External Functions .. LOGICAL LSAME INTEGER ISAMAX EXTERNAL LSAME, ISAMAX * .. * .. External Subroutines .. EXTERNAL SSCAL, SSPR, SSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SSPTRF', -INFO ) RETURN END IF * * Initialize ALPHA for use in choosing pivot block size. * ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT * IF( UPPER ) THEN * * Factorize A as U*D*U' using the upper triangle of A * * K is the main loop index, decreasing from N to 1 in steps of * 1 or 2 * K = N KC = ( N-1 )*N / 2 + 1 10 CONTINUE KNC = KC * * If K < 1, exit from loop * IF( K.LT.1 ) $ GO TO 110 KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( AP( KC+K-1 ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value * IF( K.GT.1 ) THEN IMAX = ISAMAX( K-1, AP( KC ), 1 ) COLMAX = ABS( AP( KC+IMAX-1 ) ) ELSE COLMAX = ZERO END IF * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN * * Column K is zero: set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * ROWMAX = ZERO JMAX = IMAX KX = IMAX*( IMAX+1 ) / 2 + IMAX DO 20 J = IMAX + 1, K IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN ROWMAX = ABS( AP( KX ) ) JMAX = J END IF KX = KX + J 20 CONTINUE KPC = ( IMAX-1 )*IMAX / 2 + 1 IF( IMAX.GT.1 ) THEN JMAX = ISAMAX( IMAX-1, AP( KPC ), 1 ) ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX ELSE * * interchange rows and columns K-1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * KK = K - KSTEP + 1 IF( KSTEP.EQ.2 ) $ KNC = KNC - K + 1 IF( KP.NE.KK ) THEN * * Interchange rows and columns KK and KP in the leading * submatrix A(1:k,1:k) * CALL SSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 ) KX = KPC + KP - 1 DO 30 J = KP + 1, KK - 1 KX = KX + J - 1 T = AP( KNC+J-1 ) AP( KNC+J-1 ) = AP( KX ) AP( KX ) = T 30 CONTINUE T = AP( KNC+KK-1 ) AP( KNC+KK-1 ) = AP( KPC+KP-1 ) AP( KPC+KP-1 ) = T IF( KSTEP.EQ.2 ) THEN T = AP( KC+K-2 ) AP( KC+K-2 ) = AP( KC+KP-1 ) AP( KC+KP-1 ) = T END IF END IF * * Update the leading submatrix * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column k now holds * * W(k) = U(k)*D(k) * * where U(k) is the k-th column of U * * Perform a rank-1 update of A(1:k-1,1:k-1) as * * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' * R1 = ONE / AP( KC+K-1 ) CALL SSPR( UPLO, K-1, -R1, AP( KC ), 1, AP ) * * Store U(k) in column k * CALL SSCAL( K-1, R1, AP( KC ), 1 ) ELSE * * 2-by-2 pivot block D(k): columns k and k-1 now hold * * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) * * where U(k) and U(k-1) are the k-th and (k-1)-th columns * of U * * Perform a rank-2 update of A(1:k-2,1:k-2) as * * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )' * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' * IF( K.GT.2 ) THEN * D12 = AP( K-1+( K-1 )*K / 2 ) D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 D11 = AP( K+( K-1 )*K / 2 ) / D12 T = ONE / ( D11*D22-ONE ) D12 = T / D12 * DO 50 J = K - 2, 1, -1 WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- $ AP( J+( K-1 )*K / 2 ) ) WK = D12*( D22*AP( J+( K-1 )*K / 2 )- $ AP( J+( K-2 )*( K-1 ) / 2 ) ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*WK - $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1 40 CONTINUE AP( J+( K-1 )*K / 2 ) = WK AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1 50 CONTINUE * END IF * END IF END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K-1 ) = -KP END IF * * Decrease K and return to the start of the main loop * K = K - KSTEP KC = KNC - K GO TO 10 * ELSE * * Factorize A as L*D*L' using the lower triangle of A * * K is the main loop index, increasing from 1 to N in steps of * 1 or 2 * K = 1 KC = 1 NPP = N*( N+1 ) / 2 60 CONTINUE KNC = KC * * If K > N, exit from loop * IF( K.GT.N ) $ GO TO 110 KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( AP( KC ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value * IF( K.LT.N ) THEN IMAX = K + ISAMAX( N-K, AP( KC+1 ), 1 ) COLMAX = ABS( AP( KC+IMAX-K ) ) ELSE COLMAX = ZERO END IF * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN * * Column K is zero: set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * ROWMAX = ZERO KX = KC + IMAX - K DO 70 J = K, IMAX - 1 IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN ROWMAX = ABS( AP( KX ) ) JMAX = J END IF KX = KX + N - J 70 CONTINUE KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1 IF( IMAX.LT.N ) THEN JMAX = IMAX + ISAMAX( N-IMAX, AP( KPC+1 ), 1 ) ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX ELSE * * interchange rows and columns K+1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * KK = K + KSTEP - 1 IF( KSTEP.EQ.2 ) $ KNC = KNC + N - K + 1 IF( KP.NE.KK ) THEN * * Interchange rows and columns KK and KP in the trailing * submatrix A(k:n,k:n) * IF( KP.LT.N ) $ CALL SSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 KX = KX + N - J + 1 T = AP( KNC+J-KK ) AP( KNC+J-KK ) = AP( KX ) AP( KX ) = T 80 CONTINUE T = AP( KNC ) AP( KNC ) = AP( KPC ) AP( KPC ) = T IF( KSTEP.EQ.2 ) THEN T = AP( KC+1 ) AP( KC+1 ) = AP( KC+KP-K ) AP( KC+KP-K ) = T END IF END IF * * Update the trailing submatrix * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column k now holds * * W(k) = L(k)*D(k) * * where L(k) is the k-th column of L * IF( K.LT.N ) THEN * * Perform a rank-1 update of A(k+1:n,k+1:n) as * * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' * R1 = ONE / AP( KC ) CALL SSPR( UPLO, N-K, -R1, AP( KC+1 ), 1, $ AP( KC+N-K+1 ) ) * * Store L(k) in column K * CALL SSCAL( N-K, R1, AP( KC+1 ), 1 ) END IF ELSE * * 2-by-2 pivot block D(k): columns K and K+1 now hold * * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) * * where L(k) and L(k+1) are the k-th and (k+1)-th columns * of L * IF( K.LT.N-1 ) THEN * * Perform a rank-2 update of A(k+2:n,k+2:n) as * * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )' * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )' * D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 T = ONE / ( D11*D22-ONE ) D21 = T / D21 * DO 100 J = K + 2, N WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- $ AP( J+K*( 2*N-K-1 ) / 2 ) ) WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) * DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1 90 CONTINUE * AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1 * 100 CONTINUE END IF END IF END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K+1 ) = -KP END IF * * Increase K and return to the start of the main loop * K = K + KSTEP KC = KNC + N - K + 2 GO TO 60 * END IF * 110 CONTINUE RETURN * * End of SSPTRF * END