1SUBROUTINE ZPBSVX_F95(A, B, X, UPLO, AF, FACT, EQUED, & 2 S, FERR, BERR, RCOND, INFO) 3! 4! -- LAPACK95 interface driver routine (version 3.0) -- 5! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK 6! September, 2000 7! 8! .. USE STATEMENTS .. 9 USE LA_PRECISION, ONLY: WP => DP 10 USE LA_AUXMOD, ONLY: LSAME, ERINFO 11 USE F77_LAPACK, ONLY: PBSVX_F77 => LA_PBSVX 12! .. IMPLICIT STATEMENT .. 13 IMPLICIT NONE 14! .. SCALAR ARGUMENTS .. 15 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO, FACT 16 CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED 17 INTEGER, INTENT(OUT), OPTIONAL :: INFO 18 REAL(WP), INTENT(OUT), OPTIONAL :: RCOND 19! .. ARRAY ARGUMENTS .. 20 COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:,:) 21 COMPLEX(WP), INTENT(OUT) :: X(:,:) 22 COMPLEX(WP), INTENT(INOUT), OPTIONAL, TARGET :: AF(:,:) 23 REAL(WP), INTENT(INOUT), OPTIONAL, TARGET :: S(:) 24 REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: FERR(:), BERR(:) 25!---------------------------------------------------------------------- 26! 27! Purpose 28! ======= 29! 30! LA_PBSVX computes the solution to a linear system of equations 31! A*X = B, where A has band form and is real symmetric or complex 32! Hermitian and, in either case, positive definite, and where X and B 33! are rectangular matrices or vectors. 34! LA_PBSVX can also optionally equilibrate the system if A is poorly 35! scaled, estimate the condition number of (the equilibrated) A, and 36! compute error bounds. 37! 38! ========= 39! 40! SUBROUTINE LA_PBSVX( AB, B, X, UPLO=uplo, AFB=afb, FACT=fact, & 41! EQUED=equed, S=s, FERR=ferr, BERR=berr, & 42! RCOND=rcond, INFO=info ) 43! <type>(<wp>), INTENT(INOUT) :: AB(:,:), <rhs> 44! <type>(<wp>), INTENT(OUT) :: <sol> 45! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO 46! <type>(<wp>), INTENT(INOUT), OPTIONAL :: AFB(:,:) 47! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: FACT 48! CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED 49! REAL(<wp>), INTENT(INOUT), OPTIONAL :: S(:) 50! REAL(<wp>), INTENT(OUT), OPTIONAL :: <err> 51! REAL(<wp>), INTENT(OUT), OPTIONAL :: RCOND 52! INTEGER, INTENT(OUT), OPTIONAL :: INFO 53! where 54! <type> ::= REAL | COMPLEX 55! <wp> ::= KIND(1.0) | KIND(1.0D0) 56! <rhs> ::= B(:,:) | B(:) 57! <sol> ::= X(:,:) | X(:) 58! <err> ::= FERR(:), BERR(:) | FERR, BERR 59! 60! Arguments 61! ========= 62! 63! AB (input/output) REAL or COMPLEX array, shape (:,:) with 64! size(AB,1) = kd + 1 and size(AB,2) = n, where kd is the number 65! of superdiagonals or subdiagonals and n is the order of A. 66! On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') 67! triangle of matrix A, or its equilibration, in band storage. 68! The (kd + 1) diagonals of A are stored in the rows of AB so 69! that the j-th column of A is stored in the j-th column of AB 70! as follows: 71! if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j 72! 1<=j<=n 73! if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd) 74! 1<=j<=n. 75! On exit, if FACT = 'E', then the equilibrated version of A is 76! stored in AB; otherwise, AB is unchanged. 77! B (input/output) REAL or COMPLEX array, shape (:,:) with 78! size(B,1) = n or shape (:) with size(B) = n. 79! On entry, the matrix B. 80! On exit, the scaled version of B if the system has been 81! equilibrated; otherwise, B is unchanged. 82! X (output) REAL or COMPLEX array, shape (:,:) with size(X,1)=n 83! and size(X,2) = size(B,2), or shape (:) with size(X) = n. 84! The solution matrix X . 85! UPLO Optional (input) CHARACTER(LEN=1). 86! = 'U': Upper triangle of A is stored; 87! = 'L': Lower triangle of A is stored. 88! Default value: 'U'. 89! AFB Optional (input or output) REAL or COMPLEX array, shape (:) 90! with the same size as AB. 91! If FACT = 'F' then AFB is an input argument that contains the 92! factor U or L from the Cholesky factorization of (the 93! equilibrated) A, in the same storage format as A, returned by 94! a previous call to LA_PBSVX. 95! If FACT /= 'F' then AFB is an output argument that contains 96! the factor U or L from the Cholesky factorization of (the 97! equilibrated) A in the same storage format as A. 98! FACT Optional (input) CHARACTER(LEN=1). 99! Specifies whether the factored form of the matrix A is supplied 100! on entry, and if not, whether A should be equilibrated before 101! it is factored. 102! = 'N': The matrix A will be copied to AFB and factored (no 103! equilibration). 104! = 'E': The matrix A will be equilibrated, then copied to 105! AFB and factored. 106! = 'F': AFB contains the factored form of (the equilibrated) 107! A. 108! Default value: 'N'. 109! EQUED Optional (input or output) CHARACTER(LEN=1). 110! Specifies the form of equilibration that was done. 111! EQUED is an input argument if FACT = 'F', otherwise it is an 112! output argument 113! = 'N': No equilibration (always true if FACT = 'N'). 114! = 'Y': Equilibration, i.e., A has been premultiplied and 115! postmultiplied by diag(S). 116! Default value: 'N'. 117! S Optional (input or output) REAL array, shape (:) with size(S) 118! = size(A,1). 119! The scaling factors for A. 120! S is an input argument if FACT = 'F' and EQUED = 'Y'. 121! S is an output argument if FACT = 'E' and EQUED = 'Y'. 122! FERR Optional (output) REAL array of shape (:), with size(FERR) = 123! size(X,2), or REAL scalar. 124! The estimated forward error bound for each solution vector 125! X(j) (the j-th column of the solution matrix X). If XTRUE is 126! the true solution corresponding to X(j) , FERR(j) is an 127! estimated upper bound for the magnitude of the largest element 128! in (X(j)-XTRUE) divided by the magnitude of the largest 129! element in X(j). The estimate is as reliable as the estimate 130! for RCOND, and is almost always a slight overestimate of the 131! true error. 132! BERR Optional (output) REAL array of shape (:), with size(BERR) = 133! size(X,2), or REAL scalar. 134! The componentwise relative backward error of each solution 135! vector X(j) (i.e., the smallest relative change in any element 136! of A or B that makes X(j) an exact solution). 137! RCOND Optional (output) REAL 138! The estimate of the reciprocal condition number of (the 139! equilibrated) A. If RCOND is less than the machine precision, 140! the matrix is singular to working precision. This condition is 141! indicated by a return code of INFO > 0. 142! INFO Optional (output) INTEGER 143! = 0: successful exit. 144! < 0: if INFO = -i, the i-th argument had an illegal value. 145! > 0: if INFO = i, and i is 146! <= n: the leading minor of order i of (the equilibrated) A 147! is not positive definite, so the factorization could 148! not be completed and the solution and error bounds 149! could not be computed. RCOND= 0 is returned. 150! = n+1: U or L is nonsingular, but RCOND is less than 151! machine precision, so the matrix is singular to 152! working precision. Nevertheless, the solution and 153! error bounds are computed because the computed solution 154! can be more accurate than the value of RCOND would 155! suggest. 156! If INFO is not present and an error occurs, then the program 157! is terminated with an error message. 158!---------------------------------------------------------------------- 159! .. PARAMETERS .. 160 CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_PBSVX' 161! .. LOCAL SCALARS .. 162 CHARACTER(LEN=1) :: LFACT, LUPLO, LEQUED 163 INTEGER :: LINFO, NRHS, N, ISTAT, ISTAT1, S1AF, S2AF, & 164 SS, SFERR, SBERR, KD 165 REAL(WP) :: LRCOND, MVS 166! .. LOCAL POINTERS .. 167 REAL(WP), POINTER :: RWORK(:), LS(:), LFERR(:), LBERR(:) 168 COMPLEX(WP), POINTER :: WORK(:), LAF(:, :) 169! .. INTRINSIC FUNCTIONS .. 170 INTRINSIC PRESENT, SIZE, MINVAL, TINY 171! .. EXECUTABLE STATEMENTS .. 172 LINFO = 0; ISTAT = 0 173 KD = SIZE(A,1)-1; N = SIZE(A, 2); NRHS = SIZE(B, 2) 174 IF( PRESENT(RCOND) ) RCOND = 1.0_WP 175 IF( PRESENT(FACT) )THEN; LFACT = FACT; ELSE; LFACT='N'; END IF 176 IF( PRESENT(EQUED) .AND. LSAME(LFACT,'F') )THEN; LEQUED = EQUED 177 ELSE; LEQUED='N'; END IF 178 IF( PRESENT(AF) )THEN; S1AF = SIZE(AF,1); S2AF = SIZE(AF,2) 179 ELSE; S1AF = KD+1; S2AF = N; END IF 180 IF( ( PRESENT(S) ) )THEN; SS = SIZE(S); ELSE; SS = N; END IF 181 IF( PRESENT(S) .AND. LSAME(LFACT,'F') .AND. LSAME(LEQUED,'Y') ) THEN 182 MVS = MINVAL(S); ELSE; MVS = TINY(1.0_WP); ENDIF 183 IF( PRESENT(FERR) )THEN; SFERR = SIZE(FERR); ELSE; SFERR = NRHS; END IF 184 IF( PRESENT(BERR) )THEN; SBERR = SIZE(BERR); ELSE; SBERR = NRHS; END IF 185 IF(PRESENT(UPLO))THEN; LUPLO = UPLO; ELSE; LUPLO='U'; END IF 186! .. TEST THE ARGUMENTS 187 IF( SIZE(A,1) < 0 .OR. N < 0 ) THEN; LINFO = -1 188 ELSE IF( SIZE(B, 1) /= N .OR. NRHS < 0 )THEN; LINFO = -2 189 ELSE IF( SIZE(X, 1) /= N .OR. SIZE(X, 2) /= NRHS )THEN; LINFO = -3 190 ELSE IF( .NOT.LSAME(LUPLO,'U') .AND. .NOT.LSAME(LUPLO,'L') )THEN; LINFO = -4 191 ELSE IF( S1AF /= KD+1 .OR. S2AF /= N ) THEN; LINFO = -5 192 ELSE IF( ( .NOT. ( LSAME(LFACT,'F') .OR. LSAME(LFACT,'N') .OR. & 193 LSAME(LFACT,'E') ) ) .OR. & 194 ( LSAME(LFACT,'F') .AND. .NOT.PRESENT(AF) ) )THEN; LINFO = -6 195 ELSE IF( .NOT.LSAME(LEQUED,'N') .AND. .NOT.LSAME(LEQUED,'Y') )THEN; LINFO = -7 196 ELSE IF( SS /= N .OR. LSAME(LFACT,'F') .AND. LSAME(LEQUED,'Y') & 197 .AND. MVS <= 0.0_WP )THEN; LINFO = -8 198 ELSE IF( SFERR /= NRHS )THEN; LINFO = -9 199 ELSE IF( SBERR /= NRHS )THEN; LINFO = -10 200 ELSE IF ( N > 0 )THEN 201 IF( .NOT.PRESENT(AF) ) THEN; ALLOCATE( LAF(S1AF,N), STAT=ISTAT ) 202 ELSE; LAF => AF; END IF 203 IF( ISTAT == 0 )THEN 204 IF( .NOT.PRESENT(S) )THEN; ALLOCATE( LS(N), STAT=ISTAT ) 205 ELSE; LS => S; END IF 206 END IF 207 IF( ISTAT == 0 )THEN 208 IF( .NOT.PRESENT(FERR) )THEN; ALLOCATE( LFERR(NRHS), STAT=ISTAT ) 209 ELSE; LFERR => FERR; END IF 210 END IF 211 IF( ISTAT == 0 )THEN 212 IF( .NOT.PRESENT(BERR) )THEN; ALLOCATE( LBERR(NRHS), STAT=ISTAT ) 213 ELSE; LBERR => BERR; END IF 214 END IF 215 IF( ISTAT == 0 ) ALLOCATE(WORK(2*N), RWORK(N), STAT=ISTAT ) 216 IF( ISTAT == 0 )THEN 217 CALL PBSVX_F77( LFACT, LUPLO, N, KD, NRHS, A, KD+1, LAF, S1AF, & 218 LEQUED, LS, B, N, X, N, LRCOND, & 219 LFERR, LBERR, WORK, RWORK, LINFO ) 220 ELSE; LINFO = -100; END IF 221 IF( .NOT.PRESENT(S) ) DEALLOCATE( LS, STAT=ISTAT1 ) 222 IF( .NOT.PRESENT(AF) ) DEALLOCATE( LAF, STAT=ISTAT1 ) 223 IF( .NOT.PRESENT(FERR) ) DEALLOCATE( LFERR, STAT=ISTAT1 ) 224 IF( .NOT.PRESENT(BERR) ) DEALLOCATE( LBERR, STAT=ISTAT1 ) 225 IF( PRESENT(RCOND) ) RCOND=LRCOND 226 IF( PRESENT(EQUED) .AND. .NOT.LSAME(LFACT,'F') ) EQUED=LEQUED 227 DEALLOCATE( WORK, RWORK, STAT=ISTAT1 ) 228 END IF 229 CALL ERINFO( LINFO, SRNAME, INFO, ISTAT ) 230END SUBROUTINE ZPBSVX_F95 231