1 SUBROUTINE DSYSV_F95( A, B, UPLO, IPIV, INFO ) 2! 3! -- LAPACK95 interface driver routine (version 3.0) -- 4! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK 5! September, 2000 6! 7! .. USE STATEMENTS .. 8 USE LA_PRECISION, ONLY: WP => DP 9 USE LA_AUXMOD, ONLY: ERINFO, LSAME 10 USE F77_LAPACK, ONLY: SYSV_F77 => LA_SYSV, ILAENV_F77 => ILAENV 11! .. IMPLICIT STATEMENT .. 12 IMPLICIT NONE 13! .. SCALAR ARGUMENTS .. 14 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO 15 INTEGER, INTENT(OUT), OPTIONAL :: INFO 16! .. ARRAY ARGUMENTS .. 17 INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IPIV(:) 18 REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) 19!---------------------------------------------------------------------- 20! 21! Purpose 22! ======= 23! 24! LA_SYSV computes the solution to a linear system of equations 25! A*X = B, where A is a real or complex symmetric matrix and X and B are 26! rectangular matrices or vectors. A diagonal pivoting method is used to 27! factor A as 28! A = U*D*U^T if UPLO = 'U', or A = L*D*L^T if UPLO = 'L' 29! where U (or L) is a product of permutation and unit upper (or lower) 30! triangular matrices, and D is a symmetric block diagonal matrix with 31! 1 by 1 and 2 by 2 diagonal blocks. The factored form of A is then used 32! to solve the above system. 33! LA_HESV computes the solution to a linear system of equations 34! A*X = B, where A is a complex Hermitian matrix and X and B are 35! rectangular matrices or vectors. A diagonal pivoting method is used to 36! factor A as 37! A = U*D*U^H if UPLO = 'U', or A = L*D*L^H if UPLO = 'L' 38! where U (or L) is a product of permutation and unit upper (or lower) 39! triangular matrices, and D is a complex Hermitian block diagonal 40! matrix with 1 by 1 and 2 by 2 diagonal blocks. The factored form of A 41! is then used to solve the above system. 42! 43! ========= 44! 45! SUBROUTINE LA_SYSV / LA_HESV( A, B, UPLO=uplo, & 46! IPIV=ipiv, INFO=info ) 47! <type>(<wp>), INTENT(INOUT) :: A(:,:), <rhs> 48! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO 49! INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) 50! INTEGER, INTENT(OUT), OPTIONAL :: INFO 51! where 52! <type> ::= REAL | COMPLEX 53! <wp> ::= KIND(1.0) | KIND(1.0D0) 54! <rhs> ::= B(:,:) | B(:) 55! 56! Arguments 57! ========= 58! 59! A (input/output) REAL or COMPLEX square array, shape (:,:). 60! On entry, the matrix A. 61! If UPLO = 'U', the upper triangular part of A contains the upper 62! triangular part of the matrix A, and the strictly lower 63! triangular part of A is not referenced. 64! If UPLO = 'L', the lower triangular part of A contains the lower 65! triangular part of the matrix A, and the strictly upper 66! triangular part of A is not referenced. 67! On exit, the block diagonal matrix D and the multipliers used to 68! obtain the factor U or L from the factorization of A. 69! B (input/output) REAL or COMPLEX array, shape (:,:) with 70! size(B,1) = size(A,1) or shape (:) with size(B) = size(A,1). 71! On entry, the matrix B. 72! On exit, the solution matrix X. 73! UPLO Optional (input) CHARACTER(LEN=1) 74! = 'U': Upper triangle of A is stored; 75! = 'L': Lower triangle of A is stored. 76! Default value: 'U'. 77! IPIV Optional (output) INTEGER array, shape (:) with size(IPIV) = 78! size(A,1). 79! Details of the row and column interchanges and the block 80! structure of D. 81! If IPIV(k) > 0, then rows and columns k and IPIV(k) were 82! interchanged, and D(k,k) is a 1 by 1 diagonal block. 83! If IPIV k < 0, then there are two cases: 84! 1. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 85! columns (k-1) and -IPIV(k) were interchanged and 86! D(k-1:k,k-1:k) is a 2 by 2 diagonal block. 87! 2. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0, then rows and 88! columns (k + 1) and -IPIV(k) were interchanged and 89! D(k:k+1,k:k+1) is a 2 by 2 diagonal block. 90! INFO Optional (output) INTEGER 91! = 0: successful exit. 92! < 0: if INFO = -i, the i-th argument had an illegal value. 93! > 0: if INFO = i, D(i,i) = 0. The factorization has been 94! completed, but the block diagonal matrix D is singular, so 95! the solution could not be computed. 96! If INFO is not present and an error occurs, then the program is 97! terminated with an error message. 98!----------------------------------------------------------------------- 99! .. PARAMETERS .. 100 CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_SYSV' 101 CHARACTER(LEN=6), PARAMETER :: BSNAME = 'DSYTRF' 102! .. LOCAL SCALARS .. 103 CHARACTER(LEN=1) :: LUPLO 104 INTEGER :: LINFO, ISTAT, ISTAT1, SIPIV, N, NRHS, LWORK, NB 105! .. LOCAL POINTERS .. 106 INTEGER, POINTER :: LPIV(:) 107 REAL(WP), POINTER :: WORK(:) 108! .. INTRINSIC FUNCTIONS .. 109 INTRINSIC SIZE, PRESENT 110! .. EXECUTABLE STATEMENTS .. 111 LINFO = 0; ISTAT = 0; N = SIZE(A,1); NRHS = SIZE(B,2) 112 IF( PRESENT(UPLO) ) THEN; LUPLO = UPLO; ELSE; LUPLO = 'U'; END IF 113 IF( PRESENT(IPIV) )THEN; SIPIV = SIZE(IPIV); ELSE; SIPIV = SIZE(A,1); END IF 114! .. TEST THE ARGUMENTS 115 IF( SIZE( A, 2 ) /= N .OR. N < 0 ) THEN; LINFO = -1 116 ELSE IF( SIZE( B, 1 ) /= N .OR. NRHS < 0 ) THEN; LINFO = -2 117 ELSE IF( .NOT.LSAME(LUPLO,'U') .AND. .NOT.LSAME(LUPLO,'L') )THEN; LINFO = -3 118 ELSE IF( SIPIV /= N )THEN; LINFO = -4 119 ELSE IF ( N > 0 ) THEN 120! .. DETERMINE THE WORKSPACE 121 IF( PRESENT(IPIV) )THEN; LPIV => IPIV 122 ELSE; ALLOCATE( LPIV(N), STAT = ISTAT ); END IF 123 IF( ISTAT == 0 )THEN 124 NB = ILAENV_F77( 1, BSNAME, LUPLO, N, -1, -1, -1 ) 125 IF( NB <= 1 .OR. NB >= N ) NB = 1; LWORK = N*NB 126 ALLOCATE(WORK(LWORK), STAT=ISTAT) 127 IF( ISTAT /= 0 )THEN 128 DEALLOCATE(WORK, STAT=ISTAT1); LWORK = 3*N 129 ALLOCATE(WORK(LWORK), STAT=ISTAT) 130 IF( ISTAT /= 0 ) THEN; LINFO = - 100 131 ELSE; CALL ERINFO( -200, SRNAME, LINFO ); ENDIF 132 ENDIF 133 IF ( ISTAT == 0 ) & 134! .. CALL LAPACK77 ROUTINE 135 CALL SYSV_F77( LUPLO, N, NRHS, A, N, LPIV, B, N, WORK, LWORK, LINFO ) 136 ELSE; LINFO = -100; END IF 137 IF( .NOT.PRESENT(IPIV) )DEALLOCATE(LPIV, STAT = ISTAT1 ) 138 END IF 139 CALL ERINFO( LINFO, SRNAME, INFO, ISTAT ) 140 END SUBROUTINE DSYSV_F95 141