1SUBROUTINE ZGBSVX_F95(A, B, X, KL, AFB, IPIV, FACT, TRANS, & 2 EQUED, R, C, FERR, BERR, RCOND, RPVGRW, 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: GBSVX_F77 => LA_GBSVX 12! .. IMPLICIT STATEMENT .. 13 IMPLICIT NONE 14! .. SCALAR ARGUMENTS .. 15 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT 16 CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED 17 INTEGER, INTENT(IN), OPTIONAL :: KL 18 INTEGER, INTENT(OUT), OPTIONAL :: INFO 19 REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, RPVGRW 20! .. ARRAY ARGUMENTS .. 21 COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:,:) 22 COMPLEX(WP), INTENT(OUT) :: X(:,:) 23 INTEGER, INTENT(INOUT), OPTIONAL, TARGET :: IPIV(:) 24 COMPLEX(WP), INTENT(INOUT), OPTIONAL, TARGET :: AFB(:,:) 25 REAL(WP), INTENT(INOUT), OPTIONAL, TARGET :: C(:), R(:) 26 REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: FERR(:), BERR(:) 27!---------------------------------------------------------------------- 28! 29! Purpose 30! ======= 31! 32! LA_GBSVX computes the solution to a real or complex linear system of 33! equations of the form A*X = B, A^T*X = B or A^H*X = B, where A is a 34! square band matrix and X and B are rectangular matrices or vectors. 35! LA_GBSVX can also optionally equilibrate the system if A is poorly 36! scaled, estimate the condition number of (the equilibrated) A, return 37! the pivot growth factor, and compute error bounds. 38! 39! ========= 40! 41! SUBROUTINE LA_GBSVX( AB, B, X, KL=kl, AFB=afb, IPIV=ipiv, & 42! FACT=fact, TRANS=trans, EQUED=equed, R=r, C=c, & 43! FERR=ferr, BERR=berr, RCOND=rcond, & 44! RPVGRW=rpvgrw, INFO=info ) 45! <type>(<wp>), INTENT(INOUT) :: AB(:,:), <rhs> 46! <type>(<wp>), INTENT(OUT) :: <sol> 47! INTEGER, INTENT(IN), OPTIONAL :: KL 48! <type>(<wp>), INTENT(INOUT), OPTIONAL :: AFB(:,:) 49! INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) 50! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT 51! CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED 52! REAL(<wp>), INTENT(INOUT), OPTIONAL :: C(:), R(:) 53! REAL(<wp>), INTENT(OUT), OPTIONAL :: <err>, RCOND, RPVGRW 54! INTEGER, INTENT(OUT), OPTIONAL :: INFO 55! where 56! <type> ::= REAL | COMPLEX 57! <wp> ::= KIND(1.0) | KIND(1.0D0) 58! <rhs> ::= B(:,:) | B(:) 59! <sol> ::= X(:,:) | X(:) 60! <err> ::= FERR(:), BERR(:) | FERR, BERR 61! 62! Arguments 63! ========= 64! 65! AB (input/output) REAL or COMPLEX rectangular array, shape (:,:) 66! with size(AB,1) = kl + ku + 1 and size(AB,2) = n, where kl 67! and ku are, respectively, the numbers of subdiagonals and 68! superdiagonals in the band of A, and n is the order of A. 69! On entry, the matrix A or its equilibration in band storage. 70! The (kl + ku + 1) diagonals of A are stored in rows 1 to 71! (kl + ku + 1) of AB, so that the j-th column of A is 72! stored in the j-th column of AB as follows: 73! AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl) 74! 1<=j<=n. 75! The remaining elements in AB need not be set. 76! If FACT = 'F' and EQUED /= 'N' then A has been equilibrated 77! by the scaling factors in R and/or C during the previous call 78! to LA GBSVX. 79! On exit, if FACT = 'E', the equilibrated version of A is 80! stored in AB; otherwise, AB is unchanged. 81! B (input/output) REAL or COMPLEX array, shape (:,:) with 82! size(B,1) = n or shape (:) with size(B) = n. 83! On entry, the matrix B. 84! On exit, the scaled version of B if the system has been 85! equilibrated; otherwise, B is unchanged. 86! X (output) REAL or COMPLEX array, shape (:,:) with size(X,1)=n 87! and size(X,2) = size(B,2), or shape (:) with size(X) = n. 88! The solution matrix X . 89! KL Optional (input) INTEGER. 90! The number of subdiagonals in the band of A (KL = kl). 91! The number of superdiagonals in the band is given by 92! ku = size(AB,1) - kl - 1. 93! Default value: (size(AB,1) - 1) / 2. 94! AFB Optional (input or output) REAL or COMPLEX rectangular array, 95! shape (:,:) with size(AFB,1) = 2*kl+ku+1 and size(AFB,2)=n 96! If FACT = 'F' then AFB is an input argument that contains the 97! details of the factorization of (the equilibrated) A returned 98! by a previous call to LA_GBSVX. 99! If FACT /= 'F' then AFB is an output argument that contains 100! the details of the factorization of (the equilibrated) A. U is 101! an upper triangular band matrix with (kl + ku + 1) diagonals. 102! These are stored in the first (kl + ku + 1) rows of AFB. The 103! multipliers that arise during the factorization are stored in 104! the remaining rows. 105! IPIV Optional (input or output) INTEGER array, shape (:) with 106! size(IPIV) = n. 107! If FACT = 'F' then IPIV is an input argument that contains 108! the pivot indices from the factorization of (the equilibrated) 109! A, returned by a previous call to LA_GBSVX. 110! If FACT /= 'F' then IPIV is an output argument that contains 111! the pivot indices from the factorization of (the equilibrated) 112! A. 113! FACT Optional (input) CHARACTER(LEN=1). 114! Specifies whether the factored form of the matrix A is 115! supplied on entry, and, if not, whether the matrix A should 116! be equilibrated before it is factored. 117! = 'N': The matrix A will be copied to AFB and factored (no 118! equilibration). 119! = 'E': The matrix A will be equilibrated, then copied to 120! AFB and factored. 121! = 'F': AFB and IPIV contain the factored form of (the 122! equilibrated) A. 123! Default value: 'N'. 124! TRANS Optional (input) CHARACTER(LEN=1). 125! Specifies the form of the system of equations: 126! = 'N': A*X = B (No transpose) 127! = 'T': A^T*X = B (Transpose) 128! = 'C': A^H*X = B (Conjugate transpose) 129! Default value: 'N'. 130! EQUED Optional (input or output) CHARACTER(LEN=1). 131! Specifies the form of equilibration that was done. 132! EQUED is an input argument if FACT = 'F', otherwise it is an 133! output argument: 134! = 'N': No equilibration (always true if FACT = 'N'). 135! = 'R': Row equilibration, i.e., A has been premultiplied 136! by diag(R). 137! = 'C': Column equilibration, i.e., A has been postmultiplied 138! by diag(C). 139! = 'B': Both row and column equilibration. 140! Default value: 'N'. 141! R Optional (input or output) REAL array, shape (:) with 142! size(R) = size(A,1). 143! The row scale factors for A. 144! R is an input argument if FACT = 'F' and EQUED = 'R' or 'B'. 145! R is an output argument if FACT = 'E' and EQUED = 'R' or 'B'. 146! C Optional (input or output) REAL array, shape (:) with 147! size(C) = size(A,1). 148! The column scale factors for A. 149! C is an input argument if FACT = 'F' and EQUED = 'C' or 'B'. 150! C is an output argument if FACT = 'E' and EQUED = 'C' or 'B'. 151! FERR Optional (output) REAL array of shape (:), with size(FERR) = 152! size(X,2), or REAL scalar. 153! The estimated forward error bound for each solution vector 154! X(j) (the j-th column of the solution matrix X). If XTRUE is 155! the true solution corresponding to X(j), FERR(j) is an 156! estimated upper bound for the magnitude of the largest element 157! in (X(j)-XTRUE) divided by the magnitude of the largest 158! element in X(j). The estimate is as reliable as the estimate 159! for RCOND and is almost always a slight overestimate of the 160! true error. 161! BERR Optional (output) REAL array of shape (:), with size(BERR) = 162! size(X,2), or REAL scalar. 163! The componentwise relative backward error of each solution 164! vector X(j) (i.e., the smallest relative change in any element 165! of A or B that makes X(j) an exact solution). 166! RCOND Optional (output) REAL. 167! The estimate of the reciprocal condition number of (the 168! equilibrated) A. If RCOND is less than the machine precision, 169! the matrix is singular to working precision. This condition is 170! indicated by a return code of INFO > 0. 171! RPVGRW Optional (output) REAL. 172! The reciprocal pivot growth factor ||A||inf = ||U||inf. If 173! RPVGRW is much less than 1, then the stability of the LU 174! factorization of the (equilibrated) matrix A could be poor. 175! This also means that the solution X , condition estimator 176! RCOND, and forward error bound FERR could be unreliable. If 177! the factorization fails with 0 < INFO <= size(A,1), then 178! RPVGRW contains the reciprocal pivot growth factor for the 179! leading INFO columns of A. 180! INFO Optional (output) INTEGER 181! = 0: successful exit. 182! < 0: if INFO = -i, the i-th argument had an illegal value. 183! > 0: if INFO = i, and i is 184! <= n: U(i,i) = 0. The factorization has been completed, 185! but the factor U is singular, so the solution could 186! not be computed. 187! = n+1: U is nonsingular, but RCOND is less than machine 188! precision, so the matrix is singular to working 189! precision. Nevertheless, the solution and error 190! bounds are computed because the computed solution can 191! be more accurate than the value of RCOND would suggest. 192! If INFO is not present and an error occurs, then the program 193! is terminated with an error message. 194!---------------------------------------------------------------------- 195! .. PARAMETERS .. 196 CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GBSVX' 197! .. LOCAL SCALARS .. 198 CHARACTER(LEN=1) :: LFACT, LTRANS, LEQUED 199 INTEGER :: LINFO, NRHS, N, ISTAT, ISTAT1, SIPIV, S1AFB, S2AFB, & 200 SC, SR, SFERR, SBERR, LD, LKL, LKU, LDA 201 REAL(WP) :: LRCOND, MVR, MVC 202! .. LOCAL POINTERS .. 203 INTEGER, POINTER :: LPIV(:) 204 REAL(WP), POINTER :: RWORK(:), LC(:), LR(:), LFERR(:), LBERR(:) 205 COMPLEX(WP), POINTER :: WORK(:), LAFB(:, :) 206! .. INTRINSIC FUNCTIONS .. 207 INTRINSIC MAX, PRESENT, SIZE, MINVAL, TINY 208! .. EXECUTABLE STATEMENTS .. 209 LINFO = 0; ISTAT = 0 210 LDA = SIZE(A,1); N = SIZE(A, 2); NRHS = SIZE(B, 2); LD = MAX(1,N) 211 IF( PRESENT(KL) ) THEN; LKL = KL; ELSE; LKL = (LDA-1)/2; ENDIF 212 LKU = LDA -LKL -1 213 IF( PRESENT(RCOND) ) RCOND = 1.0_WP 214 IF( PRESENT(RPVGRW) ) RPVGRW = 1.0_WP 215 IF( PRESENT(FACT) )THEN; LFACT = FACT; ELSE; LFACT='N'; END IF 216 IF( PRESENT(EQUED) .AND. LSAME(LFACT,'F') )THEN; LEQUED = EQUED 217 ELSE; LEQUED='N'; END IF 218 IF( PRESENT(IPIV) )THEN; SIPIV = SIZE(IPIV); ELSE; SIPIV = N; END IF 219 IF( PRESENT(AFB) )THEN; S1AFB = SIZE(AFB,1); S2AFB = SIZE(AFB,2) 220 ELSE; S1AFB = 2*LKL+LKU+1; S2AFB = N; END IF 221 IF( ( PRESENT(C) ) )THEN; SC = SIZE(C); ELSE; SC = N; END IF 222 IF( ( PRESENT(C) .AND. LSAME(LFACT,'F') ) .AND. & 223 ( LSAME(LEQUED,'C') .OR. LSAME(LEQUED,'B') ) )THEN; MVC = MINVAL(C) 224 ELSE; MVC = TINY(1.0_WP); END IF 225 IF( PRESENT(R) )THEN; SR = SIZE(R); ELSE; SR = N; END IF 226 IF( ( PRESENT(R) .AND. LSAME(LFACT,'F') ) .AND. & 227 ( LSAME(LEQUED,'R') .OR. LSAME(LEQUED,'B') ) )THEN; MVR = MINVAL(R) 228 ELSE; MVR = TINY(1.0_WP); END IF 229 IF( PRESENT(FERR) )THEN; SFERR = SIZE(FERR); ELSE; SFERR = NRHS; END IF 230 IF( PRESENT(BERR) )THEN; SBERR = SIZE(BERR); ELSE; SBERR = NRHS; END IF 231 IF(PRESENT(TRANS))THEN; LTRANS = TRANS; ELSE; LTRANS='N'; END IF 232! .. TEST THE ARGUMENTS 233 IF( LDA < 0 .OR. N < 0 ) THEN; LINFO = -1 234 ELSE IF( SIZE(B, 1) /= N .OR. NRHS < 0 )THEN; LINFO = -2 235 ELSE IF( SIZE(X, 1) /= N .OR. SIZE(X, 2) /= NRHS )THEN; LINFO = -3 236 ELSE IF( LKL < 0 .OR. LKU < 0 ) THEN; LINFO = -4 237 ELSE IF( S1AFB /= 2*LKL+LKU+1 .OR. S2AFB /= N ) THEN; LINFO = -5 238 ELSE IF( SIPIV /= N )THEN; LINFO = -6 239 ELSE IF( SR /= N .OR. MVR <= 0.0_WP )THEN; LINFO = -10 240 ELSE IF( SC /= N .OR. MVC <= 0.0_WP )THEN; LINFO = -11 241 ELSE IF( SFERR /= NRHS )THEN; LINFO = -12 242 ELSE IF( SBERR /= NRHS )THEN; LINFO = -13 243 ELSE IF( ( .NOT. ( LSAME(LFACT,'F') .OR. LSAME(LFACT,'N') .OR. & 244 LSAME(LFACT,'E') ) ) .OR. & 245 ( LSAME(LFACT,'F') .AND. .NOT.( PRESENT(AFB) .AND. PRESENT(IPIV) ) ) )THEN 246 LINFO = -7 247 ELSE IF( .NOT.( LSAME(LTRANS,'N') .OR. LSAME(LTRANS,'T') .OR. & 248 LSAME(LTRANS,'C') ) )THEN; LINFO = -8 249 ELSE IF( ( .NOT.( LSAME(LEQUED,'N') .OR. LSAME(LEQUED,'R') .OR. & 250 LSAME(LEQUED,'C') .OR. LSAME(LEQUED,'B') ) & 251 .AND. LSAME(LFACT,'F') ) .OR. & 252 ( ( LSAME(LEQUED,'R') .OR. LSAME(LEQUED,'B') ) .AND. & 253 .NOT.PRESENT(R) ) .OR. & 254 ( ( LSAME(LEQUED,'C') .OR. LSAME(LEQUED,'B') ) .AND. & 255 .NOT.PRESENT(C) ) )THEN; LINFO = -9 256 ELSE IF ( N > 0 )THEN 257 IF( .NOT.PRESENT(AFB) ) THEN; ALLOCATE( LAFB(S1AFB,N), STAT=ISTAT ) 258 ELSE; LAFB => AFB; END IF 259 IF( ISTAT == 0 )THEN 260 IF( .NOT.PRESENT(IPIV) )THEN; ALLOCATE( LPIV(N), STAT=ISTAT ) 261 ELSE; LPIV => IPIV; END IF 262 END IF 263 IF( ISTAT == 0 )THEN 264 IF( .NOT.PRESENT(R) )THEN; ALLOCATE( LR(N), STAT=ISTAT ) 265 ELSE; LR => R; END IF 266 END IF 267 IF( ISTAT == 0 )THEN 268 IF( .NOT.PRESENT(C) )THEN; ALLOCATE( LC(N), STAT=ISTAT ) 269 ELSE; LC => C; END IF 270 END IF 271 IF( ISTAT == 0 )THEN 272 IF( .NOT.PRESENT(FERR) )THEN; ALLOCATE( LFERR(NRHS), STAT=ISTAT ) 273 ELSE; LFERR => FERR; END IF 274 END IF 275 IF( ISTAT == 0 )THEN 276 IF( .NOT.PRESENT(BERR) )THEN; ALLOCATE( LBERR(NRHS), STAT=ISTAT ) 277 ELSE; LBERR => BERR; END IF 278 END IF 279 IF( ISTAT == 0 ) ALLOCATE(WORK(2*N), RWORK(N), STAT=ISTAT ) 280 IF( ISTAT == 0 )THEN 281 CALL GBSVX_F77( LFACT, LTRANS, N, LKL, LKU, NRHS, A, LDA, LAFB, S1AFB, & 282 LPIV, LEQUED, LR, LC, B, LD, X, LD, LRCOND, & 283 LFERR, LBERR, WORK, RWORK, LINFO ) 284 ELSE; LINFO = -100; END IF 285 IF( .NOT.PRESENT(R) ) DEALLOCATE( LR, STAT=ISTAT1 ) 286 IF( .NOT.PRESENT(C) ) DEALLOCATE( LC, STAT=ISTAT1 ) 287 IF( .NOT.PRESENT(AFB) ) DEALLOCATE( LAFB, STAT=ISTAT1 ) 288 IF( .NOT.PRESENT(IPIV) ) DEALLOCATE( LPIV, STAT=ISTAT1 ) 289 IF( .NOT.PRESENT(FERR) ) DEALLOCATE( LFERR, STAT=ISTAT1 ) 290 IF( .NOT.PRESENT(BERR) ) DEALLOCATE( LBERR, STAT=ISTAT1 ) 291 IF( PRESENT(RCOND) ) RCOND=LRCOND 292 IF( PRESENT(EQUED) .AND. .NOT.LSAME(LFACT,'F') ) EQUED=LEQUED 293 IF( PRESENT(RPVGRW) ) RPVGRW=RWORK(1) 294 DEALLOCATE( WORK, RWORK, STAT=ISTAT1 ) 295 END IF 296 CALL ERINFO( LINFO, SRNAME, INFO, ISTAT ) 297END SUBROUTINE ZGBSVX_F95 298