summaryrefslogtreecommitdiff
path: root/src/fortran/lapack/dsptrf.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/fortran/lapack/dsptrf.f')
-rw-r--r--src/fortran/lapack/dsptrf.f547
1 files changed, 547 insertions, 0 deletions
diff --git a/src/fortran/lapack/dsptrf.f b/src/fortran/lapack/dsptrf.f
new file mode 100644
index 0000000..8b8a918
--- /dev/null
+++ b/src/fortran/lapack/dsptrf.f
@@ -0,0 +1,547 @@
+ SUBROUTINE DSPTRF( 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( * )
+ DOUBLE PRECISION AP( * )
+* ..
+*
+* Purpose
+* =======
+*
+* DSPTRF 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) DOUBLE PRECISION 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 ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ DOUBLE PRECISION EIGHT, SEVTEN
+ PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
+ $ KSTEP, KX, NPP
+ DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
+ $ ROWMAX, T, WK, WKM1, WKP1
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ EXTERNAL LSAME, IDAMAX
+* ..
+* .. External Subroutines ..
+ EXTERNAL DSCAL, DSPR, DSWAP, 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( 'DSPTRF', -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 = IDAMAX( 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 = IDAMAX( 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 DSWAP( 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 DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
+*
+* Store U(k) in column k
+*
+ CALL DSCAL( 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 + IDAMAX( 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 + IDAMAX( 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 DSWAP( 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 DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
+ $ AP( KC+N-K+1 ) )
+*
+* Store L(k) in column K
+*
+ CALL DSCAL( 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 DSPTRF
+*
+ END
generated by cgit (git 2.34.1) at 2024-12-20 17:51:43 +0000