!----------------------------------------------------------------- ! ! Purpose ! ======= ! ! LA_GESVX computes the solution to a either real or complex system of ! linear equations AX = B, ! where A is a square matrix and X and B are either rectangular ! matrices or vectors. ! LA_GESVX is an expert driver routine, which can also perform the ! following functions: ! - solves A^T X = B or A^H X = B, ! - estimates the condition number of A, ! - returns the pivot growth factor, ! - refines the solution and computes forward and backward error ! bounds, ! - equilibrates the system if A is poorly scaled. ! ! Arguments ! ========= ! SUBROUTINE LA_GESVX (A, B, X, AF, IPIV, FACT, TRANS, EQUED, ! R, C, FERR, BERR, RCOND, RPVGRW, INFO) ! (), INTENT(INOUT) :: A(:,:), ! (), INTENT(OUT) :: ! (), INTENT(INOUT), OPTIONAL :: AF(:,:) ! INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) ! CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED ! REAL(), INTENT(INOUT), OPTIONAL :: R(:), C(:) ! REAL(), INTENT(OUT), OPTIONAL :: , RCOND, RPVGRW ! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: FACT, TRANS ! INTEGER, INTENT(OUT), OPTIONAL :: INFO ! where ! ::= REAL | COMPLEX ! ::= KIND(1.0) | KIND(1.0D0) ! ::= B(:,:) | B(:) ! ::= X(:,:) | X(:) ! ::= FERR(:), BERR(:) | FERR, BERR ! ! Description ! =========== ! ! The following steps are performed: ! ! 1. If FACT is not present or FACT = 'N', and EQUED is present, ! real scaling factors are computed to equilibrate the system: ! TRANS = 'N': diag(R) A diag(C) inv(diag(C)) X = diag(R) B ! TRANS = 'T': (diag(R) A diag(C))^T inv(diag(R)) X = diag(C) B ! TRANS = 'C': (diag(R) A diag(C))^H inv(diag(R)) X = diag(C) B ! Whether or not the system will be equilibrated depends on the ! scaling of the matrix A, but if equilibration is used, A is ! overwritten by diag(R) A diag(C) and B by diag(R) B (if TRANS='N') ! or diag(C) B (if TRANS = 'T' or 'C'). ! ! 2. If FACT = 'N', the LU decomposition is used to factor the ! matrix A (after equilibration if EQUED is present) as A = PLU, ! where P is a permutation matrix, L is a unit lower triangular ! matrix, and U is upper triangular. ! ! 3. The factored form of A is used to estimate the condition number ! of the matrix A. If the reciprocal of the condition number is ! less than machine precision, steps 4-6 are skipped. ! ! 4. The system of equations is solved for X using the factored form ! of A. ! ! 5. Iterative refinement is applied to improve the computed solution ! matrix and calculate error bounds and backward error estimates ! for it. ! ! 6. If equilibration was used, the matrix X is premultiplied by ! diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so ! that it solves the original system before equilibration. ! ! Arguments ! ========= ! ! A (input/output) either REAL or COMPLEX square array, ! shape (:,:), size(A,1) == size(A,2). ! If FACT is not present or FACT = 'N', ! On entry, the matrix A. ! On exit, if EQUED is present, the matrix A may have been ! overwritten by the equilibrated matrix (see EQUED). ! If FACT is present and FACT = 'F', ! On entry, the matrix A, possibly equilibrated in a previous ! call to LA_GESVX (see EQUED). ! On exit, A is unchanged. ! ! B (input/output) either REAL or COMPLEX rectangular array, ! shape either (:,:) or (:), size(B,1) or size(B) == size(A,1). ! On entry, the right hand side vector(s) of matrix B for the ! system of equations AX = B. ! On exit, if EQUED is present, B may have been scaled in ! accordance with the equilibration of A (see EQUED); ! otherwise, B is unchanged. ! ! X (output) either REAL or COMPLEX rectangular array, ! shape either (:,:) or (:), size(X,1) or size(X) == size(A,1). ! If INFO = 0, the solution matrix X to the original ! system of equations. Note that X always returns the solution ! to the original system of equations; if equilibration has been ! performed (EQUED is present and EQUED /= 'N'), this does not ! correspond to the scaled A and B. ! ! AF Optional (input/output) either REAL or COMPLEX square array, ! shape (:,:), size(AF,1) == size(AF,2) == size(A,1). ! If FACT is not present or FACT = 'N', then AF is an output ! argument and returns the factors L and U from the factorization ! A = PLU of the original matrix A, possibly equilibrated if ! EQUED is present. ! If FACT is present and FACT = 'F', then AF is an input argument ! (and must be present); on entry, it must contain the factors ! L and U of A (possibly equilibrated if EQUED is present), ! returned by a previous call to LA_GESVX. ! ! IPIV Optional (input/output) INTEGER array, shape (:), ! size(IPIV) == size(A,1). ! If FACT is not present or FACT = 'N', then IPIV is an output ! argument and returns the pivot indices from the factorization ! A = PLU of the original matrix A, possibly equilibrated if ! EQUED is present. ! If FACT is present and FACT = 'F', then IPIV is an input argument ! (and must be present); on entry, it must contain the pivot ! indices from the factorization of A (possibly equilibrated ! if EQUED is present), returned by a previous call to LA_GESVX. ! ! TRANS Optional (input) CHARACTER*1 ! If TRANS is present, it specifies the form of the system ! of equations: ! = 'N': A X = B (No transpose) ! = 'T': A^T X = B (Transpose) ! = 'C': A^H X = B (Conjugate transpose = Transpose) ! otherwise TRANS = 'N' is assumed. ! ! FACT Optional (input) CHARACTER*1 ! Specifies whether or not the factored form of the matrix A is ! supplied on entry. ! If FACT is present then: ! = 'N': the matrix A will be equilibrated if EQUED is present, ! then copied to AF and factored. ! = 'F': on entry, AF and IPIV must contain the factored form ! of A (possibly equilibrated if EQUED is present). ! otherwise FACT = 'N' is assumed. ! ! EQUED Optional (input/output) CHARACTER*1 ! If FACT is not present or FACT = 'N', then EQUED is an output ! argument. If it is present, then the matrix is equilibrated, ! and on exit EQUED specifies the scaling of A which has ! actually been performed: ! = 'N': No equilibration. ! = 'R': Row equilibration, i.e., A has been premultiplied ! by diag(R); also B has been premultiplied by diag(R) ! if TRANS = 'N'. ! = 'C': Column equilibration, i.e., A has been postmultiplied ! by diag(C); also B has been postmultiplied by diag(C) ! if TRANS = 'T' or 'C'. ! = 'B': Both row and column equilibration: combines the ! effects of EQUED = 'R' and EQUED = 'C'. ! If FACT is present and FACT = 'F', then EQUED is an input ! argument; if it is present, it specifies the equilibration ! of A which was performed in a previous call to LA_GESVX with ! FACT not present or FACT = 'N'. ! ! R Optional (input/output) REAL array, shape (:), ! size(R) == size(A,1). ! R must be present if EQUED is present and EQUED = 'R' or 'B'; ! R is not referenced if EQUED = 'N' or 'C'. ! If FACT is not present or FACT = 'N', then R is an output ! argument. If EQUED = 'R' or 'B', R returns the row scale ! factors for equilibrating A. ! If FACT is present and FACT = 'F', then R is an input ! argument. If EQUED = 'R' or 'B', R must contain the row scale ! factors for equilibrating A, returned by a previous call to ! LA_GESVX; each element of R must be positive. ! ! C Optional (input/output) REAL array, shape (:), ! size(C) == size(A,1). ! C must be present if EQUED is present and EQUED = 'C' or 'B'; ! C is not referenced if EQUED = 'N' or 'R'. ! If FACT is not present or FACT = 'N', then C is an output ! argument. If EQUED = 'C' or 'B', C returns the column scale ! factors for equilibrating A. ! If FACT is present and FACT = 'F', then C is an input ! argument. If EQUED = 'C' or 'B', C must contain the column ! scale factors for equilibrating A, returned by a previous ! call to LA_GESVX; each element of C must be positive. ! ! FERR Optional (output) either REAL array of shape (:) or REAL ! scalar. If it is an array, size(FERR) == size(X,2). ! The estimated forward error bound for each solution vector ! X(j) (the j-th column of the solution matrix X). ! If XTRUE is the true solution corresponding to X(j), FERR(j) ! is an estimated upper bound for the magnitude of the largest ! element in (X(j) - XTRUE) divided by the magnitude of the ! largest element in X(j). The estimate is as reliable as ! the estimate for RCOND, and is almost always a slight ! overestimate of the true error. ! ! BERR Optional (output) either REAL array of shape (:) or REAL ! scalar. If it is an array, size(BERR) == size(X,2). ! The componentwise relative backward error of each solution ! vector X(j) (i.e., the smallest relative change in ! any element of A or B that makes X(j) an exact solution). ! ! RCOND Optional (output) REAL ! The estimate of the reciprocal condition number of the matrix ! A after equilibration (if done). If RCOND is less than the ! machine precision (in particular, if RCOND = 0), the matrix ! is singular to working precision. This condition is ! indicated by a return code of INFO > 0, and the solution and ! error bounds are not computed. ! ! RPVGRW Optional (output) REAL. ! The reciprocal pivot growth factor norm(A)/norm(U). ! If RPVGRW is much less than 1, then the stability ! of the LU factorization of the (equilibrated) matrix A ! could be poor. This also means that the solution X, condition ! estimator RCOND, and forward error bound FERR could be ! unreliable. If factorization fails with 0 0: if INFO = k, and k is ! <= N: U(k,k) is exactly zero. The factorization has ! been completed, but the factor U is exactly ! singular, so the solution and error bounds ! could not be computed. ! = N+1: RCOND is less than machine precision. The ! factorization has been completed, but the ! matrix is singular to working precision, and ! the solution and error bounds have not been ! computed. ! If INFO is not present and an error occurs, then the program is ! terminated with an error message. !-------------------------------------- .