diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/dsysv.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/dsysv.f | 174 |
1 files changed, 174 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dsysv.f b/2.3-1/src/fortran/lapack/dsysv.f new file mode 100644 index 00000000..add53850 --- /dev/null +++ b/2.3-1/src/fortran/lapack/dsysv.f @@ -0,0 +1,174 @@ + SUBROUTINE DSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LDB, LWORK, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DSYSV computes the solution to a real system of linear equations +* A * X = B, +* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS +* matrices. +* +* The diagonal pivoting method is used to factor A as +* A = U * D * U**T, if UPLO = 'U', or +* A = L * D * L**T, if UPLO = 'L', +* 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. The factored form of A is then +* used to solve the system of equations A * X = B. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored; +* = 'L': Lower triangle of A is stored. +* +* N (input) INTEGER +* The number of linear equations, i.e., 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. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the symmetric matrix A. If UPLO = 'U', the leading +* N-by-N upper triangular part of A contains the upper +* triangular part of the matrix A, and the strictly lower +* triangular part of A is not referenced. If UPLO = 'L', the +* leading N-by-N lower triangular part of A contains the lower +* triangular part of the matrix A, and the strictly upper +* triangular part of A is not referenced. +* +* On exit, if INFO = 0, the block diagonal matrix D and the +* multipliers used to obtain the factor U or L from the +* factorization A = U*D*U**T or A = L*D*L**T as computed by +* DSYTRF. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* IPIV (output) INTEGER array, dimension (N) +* Details of the interchanges and the block structure of D, as +* determined by DSYTRF. 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. +* +* B (input/output) DOUBLE PRECISION 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). +* +* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The length of WORK. LWORK >= 1, and for best performance +* LWORK >= max(1,N*NB), where NB is the optimal blocksize for +* DSYTRF. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* 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, so the solution could not be computed. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER LWKOPT, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DSYTRF, DSYTRS, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF +* + IF( INFO.EQ.0 ) THEN + IF( N.EQ.0 ) THEN + LWKOPT = 1 + ELSE + NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) + LWKOPT = N*NB + END IF + WORK( 1 ) = LWKOPT + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYSV ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Compute the factorization A = U*D*U' or A = L*D*L'. +* + CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) + IF( INFO.EQ.0 ) THEN +* +* Solve the system A*X = B, overwriting B with X. +* + CALL DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) +* + END IF +* + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of DSYSV +* + END |