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