Dgssv Dgssvx

Dgssvx() Sub Dgssvx ( Trans As  String, N As  Long, Val() As  Double, Ptr() As  Long, Ind() As  Long, Perm_c() As  Long, Perm_r() As  Long, Etree() As  Long, Equed As  String, R() As  Double, C() As  Double, B() As  Double, X() As  Double, Rpg As  Double, RCond As  Double, FErr() As  Double, BErr() As  Double, Optional Info As  Long, Optional Format As  Long = 0, Optional Base As  Long = -1, Optional Nrhs As  Long = 1, Optional ColPerm As  String = "A", Optional Refine As  String = "N", Optional Thresh As  Double = 1, Optional SymMode As  String = "N"  ) Solves the system of linear equations A*X = B or A^T*X = B (direct method) (sparse matrix) (SuperLU) (expert driver) Purpose This routine solves the system of linear equations A*X = B or A^T*X = B, using the LU factorization, where A is a general N x N sparse matrix. The system of equations is equilibrated. Error bounds on the solution and a condition estimate are also computed. It performs the following steps. If A is stored column-wise (CSC): 1.1. If Equed = "E", scaling factors are computed to equilibrate the system as follows. If Trans = "N": diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B If Trans = "T" or "C": (diag(R)*A*diag(C))^T * 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"). 1.2. Permute the columns of A, forming A*Pc, where Pc is a permutation matrix that usually preserves sparsity. 1.3. The LU decomposition is used to factor the matrix A (after equilibration if Equed = "E") as Pr*A*Pc = L*U, with Pr determined by partial pivoting. 1.4. Compute the reciprocal pivot growth factor. 1.5. If some U(i, i) = 0, so that U is exactly singular, then the routine returns with info = i. Otherwise, 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, Info = A.NCol + 1 is returned as a warning, but the routine still goes on to solve for X and computes error bounds as described below. 1.6. The system of equations is solved for X using the factored form of A. 1.7. If Refine <> "N", iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. 1.8. 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. If A is stored row-wise (CSR), apply the above algorithm to the transpose of A: 2.1. If option Equil = YES, scaling factors are computed to equilibrate the system as follows. If Trans = "N": diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B If Trans = "T" or "C": (diag(R)*A*diag(C))^T * 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^T is overwritten by diag(R)*A^T*diag(C) and B by diag(R)*B (if Trans = NOTRANS) or diag(C)*B (if Trans = TRANSPOSE or CONJ). 2.2. Permute columns of A^T (rows of A), forming A^T*Pc, where Pc is a permutation matrix that usually preserves sparsity. 2.3. The LU decomposition is used to factor the A^T (after equilibration if Fact = YES) as Pr*A^T*Pc = L*U with the permutation Pr determined by partial pivoting. 2.4. Compute the reciprocal pivot growth factor. 2.5. If some U(i,i) = 0, so that U is exactly singular, then the routine returns with info = i. Otherwise, the factored form of A^T is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, Info = N + 1 is returned as a warning, but the routine still goes on to solve for X and computes error bounds as described below. 2.6. The system of equations is solved for X using the factored form of A^T. 2.7. If Refine <> "N", iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. 2.8. 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. Parameters [in] Trans Specifies the form of the system of equations. = "N": A*X = B (No transpose). = "T" or "C": A^T*X = B (Transpose). [in] N Number of rows and columns of the matrix A. (N >= 0) (If N = 0, returns without computation) [in,out] Val() Array Val(LVal - 1) (LVal >= Nnz) (Nnz is the number of nonzero elements of the matrix) [in] Values of nonzero elements of the matrix A. [out] Not modified if Equed = "N" on entry, or if Equed = "E" but Equed = "N" on exit. Otherwise, if Equed = "E" and Equed <> "N", A is scaled as follows: If CSC: Equed = "R": A := diag(R)*A. Equed = "C": A := A*diag(C). Equed = "B": A := diag(R)*A*diag(C). If CSR: Equed = "R": A^T := diag(R)*A^T. Equed = "C": A^T := A^T*diag(C). Equed = "B": A^T := diag(R)*A^T*diag(C). [in] Ptr() Array Ptr(LPtr - 1) (LPtr >= N + 1) Column pointers (if CSC) or row pointers (if CSR) of the matrix A. [in] Ind() Array Ind(LInd - 1) (LInd >= Nnz) Row indices (if CSC) or column indices (if CSR) of the matrix A. [in,out] Perm_c() Array Perm_c(LPerm_c - 1) (LPerm_c >= N) [in] Input permutation vector if ColPerm = "U". Otherwise, an output argument. [out] If CSC, column permutation vector which defines the permutation matrix Pc. Perm_c(i) = j means column i of A is in position j in A*Pc. Similarly, if CSR, column permutation vector which describes permutation of columns of A^T (rows of A). Input permutation vector may be overwritten by the product of the input Perm_c() and a permutation that postorders the elimination tree of Pc^T*A^T*A*Pc. Perm_c() is not changed if the elimination tree is already in postorder. [out] Perm_r() Array Perm_r(LPerm_r - 1) (LPerm_r >= N) If CSC, row permutation vector which defines the permutation matrix Pr, and is determined by partial pivoting. Perm_r(i) = j means row i of A is in position j in Pr*A. Similarly, if CSR, permutation vector which determines permutation of rows of A^T (columns of A). [out] Etree() Array Etree(LEtree - 1) (LEtree >= N) Elimination tree of Pc^T*A^T*A*Pc. Note: Etree() is a vector of parent pointers for a forest whose vertices are the integers 0 to N - 1; Etree(root) = N. [in,out] Equed Specifies the form of equilibration (scaling). [in] Specify whether or not to try equilibration. = "N": Do not equilibrate. = "E": Try to equilibrate. [out] If input Equed = "E", the result of equilibration is returned as follows. = "N": No equilibration. = "R": Row equilibration, i.e., A was premultiplied by diag(R). = "C": Column equilibration, i.e., A was postmultiplied by diag(C). = "B": Both row and column equilibration, i.e., A was replaced by diag(R)*A*diag(C). [in,out] R() Array R(LR - 1) (LR >= N) The row scale factors for A or A^T. If Equed = "R" or "B", A (if CSC) or A^T (if CSR) is multiplied on the left by diag(R). If Equed = "N" or "C", R() is not accessed. [in,out] C() Array C(LC - 1) (LC >= N) The column scale factors for A or A^T. If Equed = "C" or "B", A (if CSC) or A^T (if CSR) is multiplied on the right by diag(C). If Equed = "N" or "R", C() is not accessed. [in,out] B() Array B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= Nrhs) (2D array) or B(LB - 1) (LB >= N, Nrhs = 1) (1D array) [in] N x Nrhs right hand side matrix B. [out] If Equed = "N" on exit, B is not modified. Otherwise, B is modified as follows. If A is CSC: If Trans = "N" and Equed = "R" or "B", B is overwritten by diag(R)*B. If Trans = "T" or "C" and Equed = "C" or "B", B is overwritten by diag(C)*B. If A is CSR: If Trans = "N" and Equed = "C" or "B", B is overwritten by diag(C)*B. If Trans = "T" or "C" and Equed = "R" or "B", B is overwritten by diag(R)*B. [out] X() Array X(LX1 - 1, LX2 - 1) (LX1 >= N, LX2 >= Nrhs) (2D array) or X(LX - 1) (LX >= N, Nrhs = 1) (1D array) If Info = 0 or Info = N + 1, X contains the solution matrix to the original system of equations (before equilibration). Note that A and B are modified on exit if Equed <> "N", and the solution to the equilibrated system is inv(diag(C))*X if Trans = "N" and Equed = "C" or "B", or inv(diag(R))*X if Trans = "T" or "C" and Equed = "R" or "B". [out] Rpg The reciprocal pivot growth factor max_j(norm(A_j)/norm(U_j)). The infinity norm is used. If Rpg is much less than 1, the stability of the LU factorization could be poor. [out] RCond 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. [out] FErr() Array FErr(LFErr - 1) (LFErr >= Nrhs) 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. If Refine = "N", FErr(j) = 1. [out] BErr() Array BErr(LBErr - 1) (LBErr >= Nrhs) 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). If Refine = "N", BErr(j) = 1. [out] Info (Optional) = 0: Successful exit. < 0: The (-Info)-th argument is invalid. = i > 0 and <= N: U(i, i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed. = N + 1: U is nonsingular, but RCond is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCond would suggest. > N + 1: Number of bytes allocated when memory allocation failure occurred, plus N + 1. [in] Format (Optional) Sparse matrix format. (default = 0) = 0: CSR format. = 1: CSC format. [in] Base (Optional) Indexing of Rowptr() and Colind(). = 0: Zero-based (C style) indexing: Starting index is 0. = 1: One-based (Fortran style) indexing: Starting index is 1. (default: Assumes 1 if Ptr(0) = 1, otherwise 0) [in] Nrhs (Optional) Number of right hand sides, i.e., number of columns of the matrix B. (Nrhs >= 0) (If Nrhs = 0, returns without computation) (default = 1) [in] ColPerm (Optional) Specifies the type of column ordering to reduce fill-in. = "A": Approximate minimum degree (AMD) ordering. = "N": Natural ordering (No ordering) (Pc = I). = "M": Multiple minimum degree (MMD) ordering on A^T*A. = "P": Multiple minimum degree (MMD) ordering on A^T + A. = "U": User specified ordering. (default = "A") [in] Refine (Optional) Specifies whether to perform iterative refinement. = "N": No iterative refinement. = "S": Perform iterative refinement in single precision. = "D": Perform iterative refinement in double precision. = "X": Perform iterative refinement in extra precision. (default = "N") [in] Thresh (Optional) Specifies the threshold used for a diagonal entry to be an acceptable pivot (0 <= Thresh <= 1). (default = 1) [in] SymMode (Optional) Specifies whether to use symmetric mode. Symmetric mode gives preference to diagonal pivots, and uses an (A^T + A)-based column permutation algorithm. = "N": Do not use symmetric mode. = "S": Use symmetric mode. (default = "N") Reference SuperLU Example Program Solve the system of linear equations Ax = B and estimate the reciprocal of the condition number (RCond) of A, where ( 0.2 -0.11 -0.93 ) ( -0.3727 ) A = ( -0.32 0.81 0.37 ), B = ( 0.4319 ) ( -0.8 -0.92 -0.29 ) ( -1.4247 ) Sub Ex_Dgssvx() Const N = 3, Nnz = N * N, Nrhs = 1 Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Double, X(N - 1) As Double Dim R(N - 1) As Double, C(N - 1) As Double Dim Perm_c(N - 1) As Long, Perm_r(N - 1) As Long, Etree(N - 1) As Long Dim FErr(Nrhs - 1) As Double, BErr(Nrhs - 1) As Double Dim Equed As String, Rpg As Double, RCond As Double, Info As Long A(0) = 0.2: A(1) = -0.11: A(2) = -0.93: A(3) = -0.32: A(4) = 0.81: A(5) = 0.37: A(6) = -0.8: A(7) = -0.92: A(8) = -0.29 Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9 Ja(0) = 0: Ja(1) = 1: Ja(2) = 2: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2: Ja(6) = 0: Ja(7) = 1: Ja(8) = 2 B(0) = -0.3727: B(1) = 0.4319: B(2) = -1.4247 Equed = "E" Call Dgssvx("N", N, A(), Ia(), Ja(), Perm_c(), Perm_r(), Etree(), Equed, R(), C(), B(), X(), Rpg, RCond, FErr(), BErr(), Info) Debug.Print "X =", X(0), X(1), X(2) Debug.Print "Equed = " + Equed + ", RCond = " + CStr(RCond) + ", Info = " + CStr(Info); "" End Sub Example Results X = 0.86 0.64 0.51 Equed = N, RCond = 0.232708473186076, Info = 0