diff options
Diffstat (limited to 'SRC/zgtsv.f')
-rw-r--r-- | SRC/zgtsv.f | 173 |
1 files changed, 173 insertions, 0 deletions
diff --git a/SRC/zgtsv.f b/SRC/zgtsv.f new file mode 100644 index 00000000..ea466b38 --- /dev/null +++ b/SRC/zgtsv.f @@ -0,0 +1,173 @@ + SUBROUTINE ZGTSV( N, NRHS, DL, D, DU, B, LDB, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDB, N, NRHS +* .. +* .. Array Arguments .. + COMPLEX*16 B( LDB, * ), D( * ), DL( * ), DU( * ) +* .. +* +* Purpose +* ======= +* +* ZGTSV solves the equation +* +* A*X = B, +* +* where A is an N-by-N tridiagonal matrix, by Gaussian elimination with +* partial pivoting. +* +* Note that the equation A'*X = B may be solved by interchanging the +* order of the arguments DU and DL. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* NRHS (input) INTEGER +* The number of right hand sides, i.e., the number of columns +* of the matrix B. NRHS >= 0. +* +* DL (input/output) COMPLEX*16 array, dimension (N-1) +* On entry, DL must contain the (n-1) subdiagonal elements of +* A. +* On exit, DL is overwritten by the (n-2) elements of the +* second superdiagonal of the upper triangular matrix U from +* the LU factorization of A, in DL(1), ..., DL(n-2). +* +* D (input/output) COMPLEX*16 array, dimension (N) +* On entry, D must contain the diagonal elements of A. +* On exit, D is overwritten by the n diagonal elements of U. +* +* DU (input/output) COMPLEX*16 array, dimension (N-1) +* On entry, DU must contain the (n-1) superdiagonal elements +* of A. +* On exit, DU is overwritten by the (n-1) elements of the first +* superdiagonal of U. +* +* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) +* On entry, the N-by-NRHS right hand side matrix B. +* On exit, if INFO = 0, the N-by-NRHS solution matrix X. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, U(i,i) is exactly zero, and the solution +* has not been computed. The factorization has not been +* completed unless i = N. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER J, K + COMPLEX*16 MULT, TEMP, ZDUM +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DIMAG, MAX +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Statement Functions .. + DOUBLE PRECISION CABS1 +* .. +* .. Statement Function definitions .. + CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -2 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGTSV ', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* + DO 30 K = 1, N - 1 + IF( DL( K ).EQ.ZERO ) THEN +* +* Subdiagonal is zero, no elimination is required. +* + IF( D( K ).EQ.ZERO ) THEN +* +* Diagonal is zero: set INFO = K and return; a unique +* solution can not be found. +* + INFO = K + RETURN + END IF + ELSE IF( CABS1( D( K ) ).GE.CABS1( DL( K ) ) ) THEN +* +* No row interchange required +* + MULT = DL( K ) / D( K ) + D( K+1 ) = D( K+1 ) - MULT*DU( K ) + DO 10 J = 1, NRHS + B( K+1, J ) = B( K+1, J ) - MULT*B( K, J ) + 10 CONTINUE + IF( K.LT.( N-1 ) ) + $ DL( K ) = ZERO + ELSE +* +* Interchange rows K and K+1 +* + MULT = D( K ) / DL( K ) + D( K ) = DL( K ) + TEMP = D( K+1 ) + D( K+1 ) = DU( K ) - MULT*TEMP + IF( K.LT.( N-1 ) ) THEN + DL( K ) = DU( K+1 ) + DU( K+1 ) = -MULT*DL( K ) + END IF + DU( K ) = TEMP + DO 20 J = 1, NRHS + TEMP = B( K, J ) + B( K, J ) = B( K+1, J ) + B( K+1, J ) = TEMP - MULT*B( K+1, J ) + 20 CONTINUE + END IF + 30 CONTINUE + IF( D( N ).EQ.ZERO ) THEN + INFO = N + RETURN + END IF +* +* Back solve with the matrix U from the factorization. +* + DO 50 J = 1, NRHS + B( N, J ) = B( N, J ) / D( N ) + IF( N.GT.1 ) + $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 ) + DO 40 K = N - 2, 1, -1 + B( K, J ) = ( B( K, J )-DU( K )*B( K+1, J )-DL( K )* + $ B( K+2, J ) ) / D( K ) + 40 CONTINUE + 50 CONTINUE +* + RETURN +* +* End of ZGTSV +* + END |