COMPLEX*16 or DOUBLE COMPLEX routines for symmetric matrix
zsycon
USAGE:
rcond, info = NumRu::Lapack.zsycon( uplo, a, ipiv, anorm, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYCON( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, INFO )
* Purpose
* =======
*
* ZSYCON estimates the reciprocal of the condition number (in the
* 1-norm) of a complex symmetric matrix A using the factorization
* A = U*D*U**T or A = L*D*L**T computed by ZSYTRF.
*
* An estimate is obtained for norm(inv(A)), and the reciprocal of the
* condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by ZSYTRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* ANORM (input) DOUBLE PRECISION
* The 1-norm of the original matrix A.
*
* RCOND (output) DOUBLE PRECISION
* The reciprocal of the condition number of the matrix A,
* computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
* estimate of the 1-norm of inv(A) computed in this routine.
*
* WORK (workspace) COMPLEX*16 array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
go to the page top
zsyconv
USAGE:
info = NumRu::Lapack.zsyconv( uplo, way, a, ipiv, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
* Purpose
* =======
*
* ZSYCONV converts A given by ZHETRF into L and D or vice-versa.
* Get nondiagonal elements of D (returned in workspace) and
* apply or reverse permutation done in TRF.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* WAY (input) CHARACTER*1
* = 'C': Convert
* = 'R': Revert
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) DOUBLE COMPLEX array, dimension (LDA,N)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by ZSYTRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* WORK (workspace) DOUBLE COMPLEX array, dimension (N)
*
* LWORK (input) INTEGER
* The length of WORK. LWORK >=1.
* LWORK = N
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
go to the page top
zsyequb
USAGE:
s, scond, amax, info = NumRu::Lapack.zsyequb( uplo, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
* Purpose
* =======
*
* ZSYEQUB computes row and column scalings intended to equilibrate a
* symmetric matrix A and reduce its condition number
* (with respect to the two-norm). S contains the scale factors,
* S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
* elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This
* choice of S puts the condition number of B within a factor N of the
* smallest possible condition number over all possible diagonal
* scalings.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The N-by-N symmetric matrix whose scaling
* factors are to be computed. Only the diagonal elements of A
* are referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* S (output) DOUBLE PRECISION array, dimension (N)
* If INFO = 0, S contains the scale factors for A.
*
* SCOND (output) DOUBLE PRECISION
* If INFO = 0, S contains the ratio of the smallest S(i) to
* the largest S(i). If SCOND >= 0.1 and AMAX is neither too
* large nor too small, it is not worth scaling by S.
*
* AMAX (output) DOUBLE PRECISION
* Absolute value of largest matrix element. If AMAX is very
* close to overflow or very close to underflow, the matrix
* should be scaled.
*
* WORK (workspace) COMPLEX*16 array, dimension (3*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, the i-th diagonal element is nonpositive.
*
* Further Details
* ======= =======
*
* Reference: Livne, O.E. and Golub, G.H., "Scaling by Binormalization",
* Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004.
* DOI 10.1023/B:NUMA.0000016606.32820.69
* Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf
*
* =====================================================================
*
go to the page top
zsymv
USAGE:
y = NumRu::Lapack.zsymv( uplo, alpha, a, x, incx, beta, y, incy, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY )
* Purpose
* =======
*
* ZSYMV performs the matrix-vector operation
*
* y := alpha*A*x + beta*y,
*
* where alpha and beta are scalars, x and y are n element vectors and
* A is an n by n symmetric matrix.
*
* Arguments
* ==========
*
* UPLO (input) CHARACTER*1
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array A is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of A
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of A
* is to be referenced.
*
* Unchanged on exit.
*
* N (input) INTEGER
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA (input) COMPLEX*16
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A (input) COMPLEX*16 array, dimension ( LDA, N )
* Before entry, with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular part of the symmetric matrix and the strictly
* lower triangular part of A is not referenced.
* Before entry, with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular part of the symmetric matrix and the strictly
* upper triangular part of A is not referenced.
* Unchanged on exit.
*
* LDA (input) INTEGER
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, N ).
* Unchanged on exit.
*
* X (input) COMPLEX*16 array, dimension at least
* ( 1 + ( N - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the N-
* element vector x.
* Unchanged on exit.
*
* INCX (input) INTEGER
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA (input) COMPLEX*16
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y (input/output) COMPLEX*16 array, dimension at least
* ( 1 + ( N - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y. On exit, Y is overwritten by the updated
* vector y.
*
* INCY (input) INTEGER
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* =====================================================================
*
go to the page top
zsyr
USAGE:
a = NumRu::Lapack.zsyr( uplo, alpha, x, incx, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYR( UPLO, N, ALPHA, X, INCX, A, LDA )
* Purpose
* =======
*
* ZSYR performs the symmetric rank 1 operation
*
* A := alpha*x*( x' ) + A,
*
* where alpha is a complex scalar, x is an n element vector and A is an
* n by n symmetric matrix.
*
* Arguments
* ==========
*
* UPLO (input) CHARACTER*1
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array A is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of A
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of A
* is to be referenced.
*
* Unchanged on exit.
*
* N (input) INTEGER
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA (input) COMPLEX*16
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X (input) COMPLEX*16 array, dimension at least
* ( 1 + ( N - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the N-
* element vector x.
* Unchanged on exit.
*
* INCX (input) INTEGER
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* A (input/output) COMPLEX*16 array, dimension ( LDA, N )
* Before entry, with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular part of the symmetric matrix and the strictly
* lower triangular part of A is not referenced. On exit, the
* upper triangular part of the array A is overwritten by the
* upper triangular part of the updated matrix.
* Before entry, with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular part of the symmetric matrix and the strictly
* upper triangular part of A is not referenced. On exit, the
* lower triangular part of the array A is overwritten by the
* lower triangular part of the updated matrix.
*
* LDA (input) INTEGER
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, N ).
* Unchanged on exit.
*
* =====================================================================
*
go to the page top
zsyrfs
USAGE:
ferr, berr, info, x = NumRu::Lapack.zsyrfs( uplo, a, af, ipiv, b, x, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
* Purpose
* =======
*
* ZSYRFS improves the computed solution to a system of linear
* equations when the coefficient matrix is symmetric indefinite, and
* provides error bounds and backward error estimates for the solution.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The symmetric matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of A contains the upper triangular part
* of the matrix A, and the strictly lower triangular part of A
* is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of A contains the lower triangular part of
* the matrix A, and the strictly upper triangular part of A is
* not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input) COMPLEX*16 array, dimension (LDAF,N)
* The factored form of the matrix A. AF contains the block
* diagonal matrix D and the multipliers used to obtain the
* factor U or L from the factorization A = U*D*U**T or
* A = L*D*L**T as computed by ZSYTRF.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* B (input) COMPLEX*16 array, dimension (LDB,NRHS)
* The right hand side matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
* On entry, the solution matrix X, as computed by ZSYTRS.
* On exit, the improved solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* FERR (output) DOUBLE PRECISION array, dimension (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.
*
* BERR (output) DOUBLE PRECISION array, dimension (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).
*
* WORK (workspace) COMPLEX*16 array, dimension (2*N)
*
* RWORK (workspace) DOUBLE PRECISION array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Internal Parameters
* ===================
*
* ITMAX is the maximum number of steps of iterative refinement.
*
* =====================================================================
*
go to the page top
zsyrfsx
USAGE:
rcond, berr, err_bnds_norm, err_bnds_comp, info, s, x, params = NumRu::Lapack.zsyrfsx( uplo, equed, a, af, ipiv, s, b, x, params, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO )
* Purpose
* =======
*
* ZSYRFSX improves the computed solution to a system of linear
* equations when the coefficient matrix is symmetric indefinite, and
* provides error bounds and backward error estimates for the
* solution. In addition to normwise error bound, the code provides
* maximum componentwise error bound if possible. See comments for
* ERR_BNDS_NORM and ERR_BNDS_COMP for details of the error bounds.
*
* The original system of linear equations may have been equilibrated
* before calling this routine, as described by arguments EQUED and S
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
* Arguments
* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* defaults are acceptable. If the defaults are acceptable, users
* can pass NPARAMS = 0 which prevents the source code from accessing
* the PARAMS argument.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* EQUED (input) CHARACTER*1
* Specifies the form of equilibration that was done to A
* before calling this routine. This is needed to compute
* the solution and error bounds correctly.
* = 'N': No equilibration
* = 'Y': Both row and column equilibration, i.e., A has been
* replaced by diag(S) * A * diag(S).
* The right hand side B has been changed accordingly.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The symmetric matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of A contains the upper triangular
* part of the matrix A, and the strictly lower triangular
* part of A is not referenced. If UPLO = 'L', the leading
* N-by-N lower triangular part of A contains the lower
* triangular part of the matrix A, and the strictly upper
* triangular part of A is not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input) COMPLEX*16 array, dimension (LDAF,N)
* The factored form of the matrix A. AF contains the block
* diagonal matrix D and the multipliers used to obtain the
* factor U or L from the factorization A = U*D*U**T or A =
* L*D*L**T as computed by DSYTRF.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by DSYTRF.
*
* S (input or output) DOUBLE PRECISION array, dimension (N)
* The scale factors for A. If EQUED = 'Y', A is multiplied on
* the left and right by diag(S). S is an input argument if FACT =
* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
* = 'Y', each element of S must be positive. If S is output, each
* element of S is a power of the radix. If S is input, each element
* of S should be a power of the radix to ensure a reliable solution
* and error estimates. Scaling by powers of the radix does not cause
* rounding errors unless the result underflows or overflows.
* Rounding errors during scaling lead to refining with a matrix that
* is not equivalent to the input matrix, producing error estimates
* that may not be reliable.
*
* B (input) COMPLEX*16 array, dimension (LDB,NRHS)
* The right hand side matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
* On entry, the solution matrix X, as computed by DGETRS.
* On exit, the improved solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* RCOND (output) DOUBLE PRECISION
* Reciprocal scaled condition number. This is an estimate of the
* reciprocal Skeel condition number of the matrix A after
* equilibration (if done). If this is less than the machine
* precision (in particular, if it is zero), the matrix is singular
* to working precision. Note that the error may still be small even
* if this number is very small and the matrix appears ill-
* conditioned.
*
* BERR (output) DOUBLE PRECISION array, dimension (NRHS)
* Componentwise relative backward error. This is 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).
*
* N_ERR_BNDS (input) INTEGER
* Number of error bounds to return for each right hand side
* and each type (normwise or componentwise). See ERR_BNDS_NORM and
* ERR_BNDS_COMP below.
*
* ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
* For each right-hand side, this array contains information about
* various error bounds and condition numbers corresponding to the
* normwise relative error, which is defined as follows:
*
* Normwise relative error in the ith solution vector:
* max_j (abs(XTRUE(j,i) - X(j,i)))
* ------------------------------
* max_j abs(X(j,i))
*
* The array is indexed by the type of error information as described
* below. There currently are up to three pieces of information
* returned.
*
* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
* right-hand side.
*
* The second index in ERR_BNDS_NORM(:,err) contains the following
* three fields:
* err = 1 "Trust/don't trust" boolean. Trust the answer if the
* reciprocal condition number is less than the threshold
* sqrt(n) * dlamch('Epsilon').
*
* err = 2 "Guaranteed" error bound: The estimated forward error,
* almost certainly within a factor of 10 of the true error
* so long as the next entry is greater than the threshold
* sqrt(n) * dlamch('Epsilon'). This error bound should only
* be trusted if the previous boolean is true.
*
* err = 3 Reciprocal condition number: Estimated normwise
* reciprocal condition number. Compared with the threshold
* sqrt(n) * dlamch('Epsilon') to determine if the error
* estimate is "guaranteed". These reciprocal condition
* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
* appropriately scaled matrix Z.
* Let Z = S*A, where S scales each row by a power of the
* radix so all absolute row sums of Z are approximately 1.
*
* See Lapack Working Note 165 for further details and extra
* cautions.
*
* ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
* For each right-hand side, this array contains information about
* various error bounds and condition numbers corresponding to the
* componentwise relative error, which is defined as follows:
*
* Componentwise relative error in the ith solution vector:
* abs(XTRUE(j,i) - X(j,i))
* max_j ----------------------
* abs(X(j,i))
*
* The array is indexed by the right-hand side i (on which the
* componentwise relative error depends), and the type of error
* information as described below. There currently are up to three
* pieces of information returned for each right-hand side. If
* componentwise accuracy is not requested (PARAMS(3) = 0.0), then
* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
* the first (:,N_ERR_BNDS) entries are returned.
*
* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
* right-hand side.
*
* The second index in ERR_BNDS_COMP(:,err) contains the following
* three fields:
* err = 1 "Trust/don't trust" boolean. Trust the answer if the
* reciprocal condition number is less than the threshold
* sqrt(n) * dlamch('Epsilon').
*
* err = 2 "Guaranteed" error bound: The estimated forward error,
* almost certainly within a factor of 10 of the true error
* so long as the next entry is greater than the threshold
* sqrt(n) * dlamch('Epsilon'). This error bound should only
* be trusted if the previous boolean is true.
*
* err = 3 Reciprocal condition number: Estimated componentwise
* reciprocal condition number. Compared with the threshold
* sqrt(n) * dlamch('Epsilon') to determine if the error
* estimate is "guaranteed". These reciprocal condition
* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
* appropriately scaled matrix Z.
* Let Z = S*(A*diag(x)), where x is the solution for the
* current right-hand side and S scales each row of
* A*diag(x) by a power of the radix so all absolute row
* sums of Z are approximately 1.
*
* See Lapack Working Note 165 for further details and extra
* cautions.
*
* NPARAMS (input) INTEGER
* Specifies the number of parameters set in PARAMS. If .LE. 0, the
* PARAMS array is never referenced and default values are used.
*
* PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
* Specifies algorithm parameters. If an entry is .LT. 0.0, then
* that entry will be filled with default value used for that
* parameter. Only positions up to NPARAMS are accessed; defaults
* are used for higher-numbered parameters.
*
* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
* refinement or not.
* Default: 1.0D+0
* = 0.0 : No refinement is performed, and no error bounds are
* computed.
* = 1.0 : Use the double-precision refinement algorithm,
* possibly with doubled-single computations if the
* compilation environment does not support DOUBLE
* PRECISION.
* (other values are reserved for future use)
*
* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
* computations allowed for refinement.
* Default: 10
* Aggressive: Set to 100 to permit convergence using approximate
* factorizations or factorizations other than LU. If
* the factorization uses a technique other than
* Gaussian elimination, the guarantees in
* err_bnds_norm and err_bnds_comp may no longer be
* trustworthy.
*
* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
* will attempt to find a solution with small componentwise
* relative error in the double-precision algorithm. Positive
* is true, 0.0 is false.
* Default: 1.0 (attempt componentwise convergence)
*
* WORK (workspace) COMPLEX*16 array, dimension (2*N)
*
* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: Successful exit. The solution to every right-hand side is
* guaranteed.
* < 0: If INFO = -i, the i-th argument had an illegal value
* > 0 and <= N: U(INFO,INFO) 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. RCOND = 0
* is returned.
* = N+J: The solution corresponding to the Jth right-hand side is
* not guaranteed. The solutions corresponding to other right-
* hand sides K with K > J may not be guaranteed as well, but
* only the first such right-hand side is reported. If a small
* componentwise error is not requested (PARAMS(3) = 0.0) then
* the Jth right-hand side is the first with a normwise error
* bound that is not guaranteed (the smallest J such
* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
* the Jth right-hand side is the first with either a normwise or
* componentwise error bound that is not guaranteed (the smallest
* J such that either ERR_BNDS_NORM(J,1) = 0.0 or
* ERR_BNDS_COMP(J,1) = 0.0). See the definition of
* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
* ==================================================================
*
go to the page top
zsysv
USAGE:
ipiv, work, info, a, b = NumRu::Lapack.zsysv( uplo, a, b, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO )
* Purpose
* =======
*
* ZSYSV computes the solution to a complex system of linear equations
* A * X = B,
* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
* matrices.
*
* The diagonal pivoting method is used to factor A as
* A = U * D * U**T, if UPLO = 'U', or
* A = L * D * L**T, if UPLO = 'L',
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, and D is symmetric and block diagonal with
* 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then
* used to solve the system of equations A * X = B.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
* N-by-N upper triangular part of A contains the upper
* triangular part of the matrix A, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading N-by-N lower triangular part of A contains the lower
* triangular part of the matrix A, and the strictly upper
* triangular part of A is not referenced.
*
* On exit, if INFO = 0, the block diagonal matrix D and the
* multipliers used to obtain the factor U or L from the
* factorization A = U*D*U**T or A = L*D*L**T as computed by
* ZSYTRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (output) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D, as
* determined by ZSYTRF. If IPIV(k) > 0, then rows and columns
* k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
* diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
* then rows and columns k-1 and -IPIV(k) were interchanged and
* D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and
* IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
* -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
* diagonal block.
*
* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
* On entry, the N-by-NRHS right hand side matrix B.
* On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The length of WORK. LWORK >= 1, and for best performance
* LWORK >= max(1,N*NB), where NB is the optimal blocksize for
* ZSYTRF.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, D(i,i) is exactly zero. The factorization
* has been completed, but the block diagonal matrix D is
* exactly singular, so the solution could not be computed.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL LQUERY
INTEGER LWKOPT, NB
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
EXTERNAL XERBLA, ZSYTRF, ZSYTRS2
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
go to the page top
zsysvx
USAGE:
x, rcond, ferr, berr, work, info, af, ipiv = NumRu::Lapack.zsysvx( fact, uplo, a, af, ipiv, b, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, RWORK, INFO )
* Purpose
* =======
*
* ZSYSVX uses the diagonal pivoting factorization to compute the
* solution to a complex system of linear equations A * X = B,
* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
* matrices.
*
* Error bounds on the solution and a condition estimate are also
* provided.
*
* Description
* ===========
*
* The following steps are performed:
*
* 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
* The form of the factorization is
* A = U * D * U**T, if UPLO = 'U', or
* A = L * D * L**T, if UPLO = 'L',
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, and D is symmetric and block diagonal with
* 1-by-1 and 2-by-2 diagonal blocks.
*
* 2. If some D(i,i)=0, so that D 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 = N+1 is returned as a warning, but the routine still goes on
* to solve for X and compute error bounds as described below.
*
* 3. The system of equations is solved for X using the factored form
* of A.
*
* 4. Iterative refinement is applied to improve the computed solution
* matrix and calculate error bounds and backward error estimates
* for it.
*
* Arguments
* =========
*
* FACT (input) CHARACTER*1
* Specifies whether or not the factored form of A has been
* supplied on entry.
* = 'F': On entry, AF and IPIV contain the factored form
* of A. A, AF and IPIV will not be modified.
* = 'N': The matrix A will be copied to AF and factored.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The symmetric matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of A contains the upper triangular part
* of the matrix A, and the strictly lower triangular part of A
* is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of A contains the lower triangular part of
* the matrix A, and the strictly upper triangular part of A is
* not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
* If FACT = 'F', then AF is an input argument and on entry
* contains the block diagonal matrix D and the multipliers used
* to obtain the factor U or L from the factorization
* A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF.
*
* If FACT = 'N', then AF is an output argument and on exit
* returns the block diagonal matrix D and the multipliers used
* to obtain the factor U or L from the factorization
* A = U*D*U**T or A = L*D*L**T.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input or output) INTEGER array, dimension (N)
* If FACT = 'F', then IPIV is an input argument and on entry
* contains details of the interchanges and the block structure
* of D, as determined by ZSYTRF.
* If IPIV(k) > 0, then rows and columns k and IPIV(k) were
* interchanged and D(k,k) is a 1-by-1 diagonal block.
* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
*
* If FACT = 'N', then IPIV is an output argument and on exit
* contains details of the interchanges and the block structure
* of D, as determined by ZSYTRF.
*
* B (input) COMPLEX*16 array, dimension (LDB,NRHS)
* The N-by-NRHS right hand side matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (output) COMPLEX*16 array, dimension (LDX,NRHS)
* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* RCOND (output) DOUBLE PRECISION
* The estimate of the reciprocal condition number of the matrix
* A. 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.
*
* FERR (output) DOUBLE PRECISION array, dimension (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.
*
* BERR (output) DOUBLE PRECISION array, dimension (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).
*
* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The length of WORK. LWORK >= max(1,2*N), and for best
* performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
* NB is the optimal blocksize for ZSYTRF.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* RWORK (workspace) DOUBLE PRECISION array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, and i is
* <= N: D(i,i) is exactly zero. The factorization
* has been completed but the factor D is exactly
* singular, so the solution and error bounds could
* not be computed. RCOND = 0 is returned.
* = N+1: D 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.
*
* =====================================================================
*
go to the page top
zsysvxx
USAGE:
x, rcond, rpvgrw, berr, err_bnds_norm, err_bnds_comp, info, a, af, ipiv, equed, s, b, params = NumRu::Lapack.zsysvxx( fact, uplo, a, af, ipiv, equed, s, b, params, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO )
* Purpose
* =======
*
* ZSYSVXX uses the diagonal pivoting factorization to compute the
* solution to a complex*16 system of linear equations A * X = B, where
* A is an N-by-N symmetric matrix and X and B are N-by-NRHS
* matrices.
*
* If requested, both normwise and maximum componentwise error bounds
* are returned. ZSYSVXX will return a solution with a tiny
* guaranteed error (O(eps) where eps is the working machine
* precision) unless the matrix is very ill-conditioned, in which
* case a warning is returned. Relevant condition numbers also are
* calculated and returned.
*
* ZSYSVXX accepts user-provided factorizations and equilibration
* factors; see the definitions of the FACT and EQUED options.
* Solving with refinement and using a factorization from a previous
* ZSYSVXX call will also produce a solution with either O(eps)
* errors or warnings, but we cannot make that claim for general
* user-provided factorizations and equilibration factors if they
* differ from what ZSYSVXX would itself produce.
*
* Description
* ===========
*
* The following steps are performed:
*
* 1. If FACT = 'E', double precision scaling factors are computed to equilibrate
* the system:
*
* diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*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(S)*A*diag(S) and B by diag(S)*B.
*
* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
* the matrix A (after equilibration if FACT = 'E') as
*
* A = U * D * U**T, if UPLO = 'U', or
* A = L * D * L**T, if UPLO = 'L',
*
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, and D is symmetric and block diagonal with
* 1-by-1 and 2-by-2 diagonal blocks.
*
* 3. If some D(i,i)=0, so that D 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 (see
* argument RCOND). If the reciprocal of the condition number is
* less than machine precision, the routine still goes on to solve
* for X and compute error bounds as described below.
*
* 4. The system of equations is solved for X using the factored form
* of A.
*
* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
* the routine will use iterative refinement to try to get a small
* error and error bounds. Refinement calculates the residual to at
* least twice the working precision.
*
* 6. If equilibration was used, the matrix X is premultiplied by
* diag(R) so that it solves the original system before
* equilibration.
*
* Arguments
* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* defaults are acceptable. If the defaults are acceptable, users
* can pass NPARAMS = 0 which prevents the source code from accessing
* the PARAMS argument.
*
* FACT (input) CHARACTER*1
* Specifies whether or not the factored form of the matrix A is
* supplied on entry, and if not, whether the matrix A should be
* equilibrated before it is factored.
* = 'F': On entry, AF and IPIV contain the factored form of A.
* If EQUED is not 'N', the matrix A has been
* equilibrated with scaling factors given by S.
* A, AF, and IPIV are not modified.
* = 'N': The matrix A will be copied to AF and factored.
* = 'E': The matrix A will be equilibrated if necessary, then
* copied to AF and factored.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
* The symmetric matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of A contains the upper triangular
* part of the matrix A, and the strictly lower triangular
* part of A is not referenced. If UPLO = 'L', the leading
* N-by-N lower triangular part of A contains the lower
* triangular part of the matrix A, and the strictly upper
* triangular part of A is not referenced.
*
* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
* diag(S)*A*diag(S).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
* If FACT = 'F', then AF is an input argument and on entry
* contains the block diagonal matrix D and the multipliers
* used to obtain the factor U or L from the factorization A =
* U*D*U**T or A = L*D*L**T as computed by DSYTRF.
*
* If FACT = 'N', then AF is an output argument and on exit
* returns the block diagonal matrix D and the multipliers
* used to obtain the factor U or L from the factorization A =
* U*D*U**T or A = L*D*L**T.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input or output) INTEGER array, dimension (N)
* If FACT = 'F', then IPIV is an input argument and on entry
* contains details of the interchanges and the block
* structure of D, as determined by DSYTRF. If IPIV(k) > 0,
* then rows and columns k and IPIV(k) were interchanged and
* D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and
* IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and
* -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2
* diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,
* then rows and columns k+1 and -IPIV(k) were interchanged
* and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
*
* If FACT = 'N', then IPIV is an output argument and on exit
* contains details of the interchanges and the block
* structure of D, as determined by DSYTRF.
*
* EQUED (input or output) CHARACTER*1
* Specifies the form of equilibration that was done.
* = 'N': No equilibration (always true if FACT = 'N').
* = 'Y': Both row and column equilibration, i.e., A has been
* replaced by diag(S) * A * diag(S).
* EQUED is an input argument if FACT = 'F'; otherwise, it is an
* output argument.
*
* S (input or output) DOUBLE PRECISION array, dimension (N)
* The scale factors for A. If EQUED = 'Y', A is multiplied on
* the left and right by diag(S). S is an input argument if FACT =
* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
* = 'Y', each element of S must be positive. If S is output, each
* element of S is a power of the radix. If S is input, each element
* of S should be a power of the radix to ensure a reliable solution
* and error estimates. Scaling by powers of the radix does not cause
* rounding errors unless the result underflows or overflows.
* Rounding errors during scaling lead to refining with a matrix that
* is not equivalent to the input matrix, producing error estimates
* that may not be reliable.
*
* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
* On entry, the N-by-NRHS right hand side matrix B.
* On exit,
* if EQUED = 'N', B is not modified;
* if EQUED = 'Y', B is overwritten by diag(S)*B;
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (output) COMPLEX*16 array, dimension (LDX,NRHS)
* If INFO = 0, the N-by-NRHS solution matrix X to the original
* system of equations. Note that A and B are modified on exit if
* EQUED .ne. 'N', and the solution to the equilibrated system is
* inv(diag(S))*X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* RCOND (output) DOUBLE PRECISION
* Reciprocal scaled condition number. This is an estimate of the
* reciprocal Skeel condition number of the matrix A after
* equilibration (if done). If this is less than the machine
* precision (in particular, if it is zero), the matrix is singular
* to working precision. Note that the error may still be small even
* if this number is very small and the matrix appears ill-
* conditioned.
*
* RPVGRW (output) DOUBLE PRECISION
* Reciprocal pivot growth. On exit, this contains the reciprocal
* pivot growth factor norm(A)/norm(U). The "max absolute element"
* norm is used. If this 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, estimated condition numbers,
* and error bounds could be unreliable. If factorization fails with
* 0 0 and <= N: U(INFO,INFO) 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. RCOND = 0
* is returned.
* = N+J: The solution corresponding to the Jth right-hand side is
* not guaranteed. The solutions corresponding to other right-
* hand sides K with K > J may not be guaranteed as well, but
* only the first such right-hand side is reported. If a small
* componentwise error is not requested (PARAMS(3) = 0.0) then
* the Jth right-hand side is the first with a normwise error
* bound that is not guaranteed (the smallest J such
* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
* the Jth right-hand side is the first with either a normwise or
* componentwise error bound that is not guaranteed (the smallest
* J such that either ERR_BNDS_NORM(J,1) = 0.0 or
* ERR_BNDS_COMP(J,1) = 0.0). See the definition of
* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
* ==================================================================
*
go to the page top
zsyswapr
USAGE:
a = NumRu::Lapack.zsyswapr( uplo, a, i1, i2, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYSWAPR( UPLO, N, A, I1, I2)
* Purpose
* =======
*
* ZSYSWAPR applies an elementary permutation on the rows and the columns of
* a symmetric matrix.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) DOUBLE COMPLEX array, dimension (LDA,N)
* On entry, the NB diagonal matrix D and the multipliers
* used to obtain the factor U or L as computed by ZSYTRF.
*
* On exit, if INFO = 0, the (symmetric) inverse of the original
* matrix. If UPLO = 'U', the upper triangular part of the
* inverse is formed and the part of A below the diagonal is not
* referenced; if UPLO = 'L' the lower triangular part of the
* inverse is formed and the part of A above the diagonal is
* not referenced.
*
* I1 (input) INTEGER
* Index of the first row to swap
*
* I2 (input) INTEGER
* Index of the second row to swap
*
* =====================================================================
*
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I
DOUBLE COMPLEX TMP
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL ZSWAP
* ..
go to the page top
zsytf2
USAGE:
ipiv, info, a = NumRu::Lapack.zsytf2( uplo, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO )
* Purpose
* =======
*
* ZSYTF2 computes the factorization of a complex symmetric matrix A
* using the Bunch-Kaufman diagonal pivoting method:
*
* A = U*D*U' or A = L*D*L'
*
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, U' is the transpose of U, and D is symmetric and
* block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
*
* This is the unblocked version of the algorithm, calling Level 2 BLAS.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the upper or lower triangular part of the
* symmetric matrix A is stored:
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
* n-by-n upper triangular part of A contains the upper
* triangular part of the matrix A, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading n-by-n lower triangular part of A contains the lower
* triangular part of the matrix A, and the strictly upper
* triangular part of A is not referenced.
*
* On exit, the block diagonal matrix D and the multipliers used
* to obtain the factor U or L (see below for further details).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (output) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D.
* If IPIV(k) > 0, then rows and columns k and IPIV(k) were
* interchanged and D(k,k) is a 1-by-1 diagonal block.
* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -k, the k-th argument had an illegal value
* > 0: if INFO = k, D(k,k) is exactly zero. The factorization
* has been completed, but the block diagonal matrix D is
* exactly singular, and division by zero will occur if it
* is used to solve a system of equations.
*
* Further Details
* ===============
*
* 09-29-06 - patch from
* Bobby Cheng, MathWorks
*
* Replace l.209 and l.377
* IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
* by
* IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
*
* 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
* Company
*
* If UPLO = 'U', then A = U*D*U', where
* U = P(n)*U(n)* ... *P(k)U(k)* ...,
* i.e., U is a product of terms P(k)*U(k), where k decreases from n to
* 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
* defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
*
* ( I v 0 ) k-s
* U(k) = ( 0 I 0 ) s
* ( 0 0 I ) n-k
* k-s s n-k
*
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
* If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
* and A(k,k), and v overwrites A(1:k-2,k-1:k).
*
* If UPLO = 'L', then A = L*D*L', where
* L = P(1)*L(1)* ... *P(k)*L(k)* ...,
* i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
* n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
* defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
*
* ( I 0 0 ) k-1
* L(k) = ( 0 I 0 ) s
* ( 0 v I ) n-k-s+1
* k-1 s n-k-s+1
*
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
* If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
* and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
*
* =====================================================================
*
go to the page top
zsytrf
USAGE:
ipiv, work, info, a = NumRu::Lapack.zsytrf( uplo, a, lwork, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
* Purpose
* =======
*
* ZSYTRF computes the factorization of a complex symmetric matrix A
* using the Bunch-Kaufman diagonal pivoting method. The form of the
* factorization is
*
* A = U*D*U**T or A = L*D*L**T
*
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, and D is symmetric and block diagonal with
* with 1-by-1 and 2-by-2 diagonal blocks.
*
* This is the blocked version of the algorithm, calling Level 3 BLAS.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* = 'L': Lower triangle of A is stored.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
* N-by-N upper triangular part of A contains the upper
* triangular part of the matrix A, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading N-by-N lower triangular part of A contains the lower
* triangular part of the matrix A, and the strictly upper
* triangular part of A is not referenced.
*
* On exit, the block diagonal matrix D and the multipliers used
* to obtain the factor U or L (see below for further details).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (output) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D.
* If IPIV(k) > 0, then rows and columns k and IPIV(k) were
* interchanged and D(k,k) is a 1-by-1 diagonal block.
* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
*
* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The length of WORK. LWORK >=1. For best performance
* LWORK >= N*NB, where NB is the block size returned by ILAENV.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, D(i,i) is exactly zero. The factorization
* has been completed, but the block diagonal matrix D is
* exactly singular, and division by zero will occur if it
* is used to solve a system of equations.
*
* Further Details
* ===============
*
* If UPLO = 'U', then A = U*D*U', where
* U = P(n)*U(n)* ... *P(k)U(k)* ...,
* i.e., U is a product of terms P(k)*U(k), where k decreases from n to
* 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
* defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
*
* ( I v 0 ) k-s
* U(k) = ( 0 I 0 ) s
* ( 0 0 I ) n-k
* k-s s n-k
*
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
* If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
* and A(k,k), and v overwrites A(1:k-2,k-1:k).
*
* If UPLO = 'L', then A = L*D*L', where
* L = P(1)*L(1)* ... *P(k)*L(k)* ...,
* i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
* n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
* defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
*
* ( I 0 0 ) k-1
* L(k) = ( 0 I 0 ) s
* ( 0 v I ) n-k-s+1
* k-1 s n-k-s+1
*
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
* If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
* and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL LQUERY, UPPER
INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
EXTERNAL XERBLA, ZLASYF, ZSYTF2
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
go to the page top
zsytri
USAGE:
info, a = NumRu::Lapack.zsytri( uplo, a, ipiv, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
* Purpose
* =======
*
* ZSYTRI computes the inverse of a complex symmetric indefinite matrix
* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
* ZSYTRF.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
* On entry, the block diagonal matrix D and the multipliers
* used to obtain the factor U or L as computed by ZSYTRF.
*
* On exit, if INFO = 0, the (symmetric) inverse of the original
* matrix. If UPLO = 'U', the upper triangular part of the
* inverse is formed and the part of A below the diagonal is not
* referenced; if UPLO = 'L' the lower triangular part of the
* inverse is formed and the part of A above the diagonal is
* not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* WORK (workspace) COMPLEX*16 array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
* inverse could not be computed.
*
* =====================================================================
*
go to the page top
zsytri2
USAGE:
info, a = NumRu::Lapack.zsytri2( uplo, a, ipiv, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRI2( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
* Purpose
* =======
*
* ZSYTRI2 computes the inverse of a complex symmetric indefinite matrix
* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
* ZSYTRF. ZSYTRI2 sets the LEADING DIMENSION of the workspace
* before calling ZSYTRI2X that actually computes the inverse.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) DOUBLE COMPLEX array, dimension (LDA,N)
* On entry, the NB diagonal matrix D and the multipliers
* used to obtain the factor U or L as computed by ZSYTRF.
*
* On exit, if INFO = 0, the (symmetric) inverse of the original
* matrix. If UPLO = 'U', the upper triangular part of the
* inverse is formed and the part of A below the diagonal is not
* referenced; if UPLO = 'L' the lower triangular part of the
* inverse is formed and the part of A above the diagonal is
* not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the NB structure of D
* as determined by ZSYTRF.
*
* WORK (workspace) DOUBLE COMPLEX array, dimension (N+NB+1)*(NB+3)
*
* LWORK (input) INTEGER
* The dimension of the array WORK.
* WORK is size >= (N+NB+1)*(NB+3)
* If LDWORK = -1, then a workspace query is assumed; the routine
* calculates:
* - the optimal size of the WORK array, returns
* this value as the first entry of the WORK array,
* - and no error message related to LDWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
* inverse could not be computed.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL UPPER, LQUERY
INTEGER MINSIZE, NBMAX
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
EXTERNAL ZSYTRI2X
* ..
go to the page top
zsytri2x
USAGE:
info, a = NumRu::Lapack.zsytri2x( uplo, a, ipiv, nb, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
* Purpose
* =======
*
* ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
* ZSYTRF.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) DOUBLE COMPLEX array, dimension (LDA,N)
* On entry, the NNB diagonal matrix D and the multipliers
* used to obtain the factor U or L as computed by ZSYTRF.
*
* On exit, if INFO = 0, the (symmetric) inverse of the original
* matrix. If UPLO = 'U', the upper triangular part of the
* inverse is formed and the part of A below the diagonal is not
* referenced; if UPLO = 'L' the lower triangular part of the
* inverse is formed and the part of A above the diagonal is
* not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the NNB structure of D
* as determined by ZSYTRF.
*
* WORK (workspace) DOUBLE COMPLEX array, dimension (N+NNB+1,NNB+3)
*
* NB (input) INTEGER
* Block size
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
* inverse could not be computed.
*
* =====================================================================
*
go to the page top
zsytrs
USAGE:
info, b = NumRu::Lapack.zsytrs( uplo, a, ipiv, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
* Purpose
* =======
*
* ZSYTRS solves a system of linear equations A*X = B with a complex
* symmetric matrix A using the factorization A = U*D*U**T or
* A = L*D*L**T computed by ZSYTRF.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input) COMPLEX*16 array, dimension (LDA,N)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by ZSYTRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
* On entry, the right hand side matrix B.
* On exit, the solution matrix X.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
go to the page top
zsytrs2
USAGE:
info, b = NumRu::Lapack.zsytrs2( uplo, a, ipiv, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
* Purpose
* =======
*
* ZSYTRS2 solves a system of linear equations A*X = B with a real
* symmetric matrix A using the factorization A = U*D*U**T or
* A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the details of the factorization are stored
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input) DOUBLE COMPLEX array, dimension (LDA,N)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by ZSYTRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by ZSYTRF.
*
* B (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS)
* On entry, the right hand side matrix B.
* On exit, the solution matrix X.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* WORK (workspace) REAL array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
go to the page top
back to matrix types
back to data types