COMPLEX routines for triangular (or in some cases quasi-triangular) matrix
ctrcon
USAGE:
rcond, info = NumRu::Lapack.ctrcon( norm, uplo, diag, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK, RWORK, INFO )
* Purpose
* =======
*
* CTRCON estimates the reciprocal of the condition number of a
* triangular matrix A, in either the 1-norm or the infinity-norm.
*
* The norm of A is computed and an estimate is obtained for
* norm(inv(A)), then the reciprocal of the condition number is
* computed as
* RCOND = 1 / ( norm(A) * norm(inv(A)) ).
*
* Arguments
* =========
*
* NORM (input) CHARACTER*1
* Specifies whether the 1-norm condition number or the
* infinity-norm condition number is required:
* = '1' or 'O': 1-norm;
* = 'I': Infinity-norm.
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) COMPLEX array, dimension (LDA,N)
* The triangular matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of the array A contains the upper
* triangular matrix, and the strictly lower triangular part of
* A is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of the array A contains the lower triangular
* matrix, and the strictly upper triangular part of A is not
* referenced. If DIAG = 'U', the diagonal elements of A are
* also not referenced and are assumed to be 1.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* RCOND (output) REAL
* The reciprocal of the condition number of the matrix A,
* computed as RCOND = 1/(norm(A) * norm(inv(A))).
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* RWORK (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
ctrevc
USAGE:
m, info, t, vl, vr = NumRu::Lapack.ctrevc( side, howmny, select, t, vl, vr, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
* Purpose
* =======
*
* CTREVC computes some or all of the right and/or left eigenvectors of
* a complex upper triangular matrix T.
* Matrices of this type are produced by the Schur factorization of
* a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR.
*
* The right eigenvector x and the left eigenvector y of T corresponding
* to an eigenvalue w are defined by:
*
* T*x = w*x, (y**H)*T = w*(y**H)
*
* where y**H denotes the conjugate transpose of the vector y.
* The eigenvalues are not input to this routine, but are read directly
* from the diagonal of T.
*
* This routine returns the matrices X and/or Y of right and left
* eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
* input matrix. If Q is the unitary factor that reduces a matrix A to
* Schur form T, then Q*X and Q*Y are the matrices of right and left
* eigenvectors of A.
*
* Arguments
* =========
*
* SIDE (input) CHARACTER*1
* = 'R': compute right eigenvectors only;
* = 'L': compute left eigenvectors only;
* = 'B': compute both right and left eigenvectors.
*
* HOWMNY (input) CHARACTER*1
* = 'A': compute all right and/or left eigenvectors;
* = 'B': compute all right and/or left eigenvectors,
* backtransformed using the matrices supplied in
* VR and/or VL;
* = 'S': compute selected right and/or left eigenvectors,
* as indicated by the logical array SELECT.
*
* SELECT (input) LOGICAL array, dimension (N)
* If HOWMNY = 'S', SELECT specifies the eigenvectors to be
* computed.
* The eigenvector corresponding to the j-th eigenvalue is
* computed if SELECT(j) = .TRUE..
* Not referenced if HOWMNY = 'A' or 'B'.
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input/output) COMPLEX array, dimension (LDT,N)
* The upper triangular matrix T. T is modified, but restored
* on exit.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* VL (input/output) COMPLEX array, dimension (LDVL,MM)
* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
* contain an N-by-N matrix Q (usually the unitary matrix Q of
* Schur vectors returned by CHSEQR).
* On exit, if SIDE = 'L' or 'B', VL contains:
* if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
* if HOWMNY = 'B', the matrix Q*Y;
* if HOWMNY = 'S', the left eigenvectors of T specified by
* SELECT, stored consecutively in the columns
* of VL, in the same order as their
* eigenvalues.
* Not referenced if SIDE = 'R'.
*
* LDVL (input) INTEGER
* The leading dimension of the array VL. LDVL >= 1, and if
* SIDE = 'L' or 'B', LDVL >= N.
*
* VR (input/output) COMPLEX array, dimension (LDVR,MM)
* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
* contain an N-by-N matrix Q (usually the unitary matrix Q of
* Schur vectors returned by CHSEQR).
* On exit, if SIDE = 'R' or 'B', VR contains:
* if HOWMNY = 'A', the matrix X of right eigenvectors of T;
* if HOWMNY = 'B', the matrix Q*X;
* if HOWMNY = 'S', the right eigenvectors of T specified by
* SELECT, stored consecutively in the columns
* of VR, in the same order as their
* eigenvalues.
* Not referenced if SIDE = 'L'.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR. LDVR >= 1, and if
* SIDE = 'R' or 'B'; LDVR >= N.
*
* MM (input) INTEGER
* The number of columns in the arrays VL and/or VR. MM >= M.
*
* M (output) INTEGER
* The number of columns in the arrays VL and/or VR actually
* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
* is set to N. Each selected eigenvector occupies one
* column.
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* RWORK (workspace) REAL array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* The algorithm used in this program is basically backward (forward)
* substitution, with scaling to make the the code robust against
* possible overflow.
*
* Each eigenvector is normalized so that the element of largest
* magnitude has magnitude 1; here the magnitude of a complex number
* (x,y) is taken to be |x| + |y|.
*
* =====================================================================
*
go to the page top
ctrexc
USAGE:
info, t, q = NumRu::Lapack.ctrexc( compq, t, q, ifst, ilst, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO )
* Purpose
* =======
*
* CTREXC reorders the Schur factorization of a complex matrix
* A = Q*T*Q**H, so that the diagonal element of T with row index IFST
* is moved to row ILST.
*
* The Schur form T is reordered by a unitary similarity transformation
* Z**H*T*Z, and optionally the matrix Q of Schur vectors is updated by
* postmultplying it with Z.
*
* Arguments
* =========
*
* COMPQ (input) CHARACTER*1
* = 'V': update the matrix Q of Schur vectors;
* = 'N': do not update Q.
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input/output) COMPLEX array, dimension (LDT,N)
* On entry, the upper triangular matrix T.
* On exit, the reordered upper triangular matrix.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* Q (input/output) COMPLEX array, dimension (LDQ,N)
* On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
* On exit, if COMPQ = 'V', Q has been postmultiplied by the
* unitary transformation matrix Z which reorders T.
* If COMPQ = 'N', Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= max(1,N).
*
* IFST (input) INTEGER
* ILST (input) INTEGER
* Specify the reordering of the diagonal elements of T:
* The element with row index IFST is moved to row ILST by a
* sequence of transpositions between adjacent elements.
* 1 <= IFST <= N; 1 <= ILST <= N.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL WANTQ
INTEGER K, M1, M2, M3
REAL CS
COMPLEX SN, T11, T22, TEMP
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL CLARTG, CROT, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC CONJG, MAX
* ..
go to the page top
ctrrfs
USAGE:
ferr, berr, info = NumRu::Lapack.ctrrfs( uplo, trans, diag, a, b, x, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
* Purpose
* =======
*
* CTRRFS provides error bounds and backward error estimates for the
* solution to a system of linear equations with a triangular
* coefficient matrix.
*
* The solution matrix X must be computed by CTRTRS or some other
* means before entering this routine. CTRRFS does not do iterative
* refinement because doing so cannot improve the backward error.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* TRANS (input) CHARACTER*1
* 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)
*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
*
* 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 array, dimension (LDA,N)
* The triangular matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of the array A contains the upper
* triangular matrix, and the strictly lower triangular part of
* A is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of the array A contains the lower triangular
* matrix, and the strictly upper triangular part of A is not
* referenced. If DIAG = 'U', the diagonal elements of A are
* also not referenced and are assumed to be 1.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input) COMPLEX 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) COMPLEX array, dimension (LDX,NRHS)
* The solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* FERR (output) REAL 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) REAL 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 array, dimension (2*N)
*
* RWORK (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
ctrsen
USAGE:
w, m, s, sep, work, info, t, q = NumRu::Lapack.ctrsen( job, compq, select, t, q, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO )
* Purpose
* =======
*
* CTRSEN reorders the Schur factorization of a complex matrix
* A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
* the leading positions on the diagonal of the upper triangular matrix
* T, and the leading columns of Q form an orthonormal basis of the
* corresponding right invariant subspace.
*
* Optionally the routine computes the reciprocal condition numbers of
* the cluster of eigenvalues and/or the invariant subspace.
*
* Arguments
* =========
*
* JOB (input) CHARACTER*1
* Specifies whether condition numbers are required for the
* cluster of eigenvalues (S) or the invariant subspace (SEP):
* = 'N': none;
* = 'E': for eigenvalues only (S);
* = 'V': for invariant subspace only (SEP);
* = 'B': for both eigenvalues and invariant subspace (S and
* SEP).
*
* COMPQ (input) CHARACTER*1
* = 'V': update the matrix Q of Schur vectors;
* = 'N': do not update Q.
*
* SELECT (input) LOGICAL array, dimension (N)
* SELECT specifies the eigenvalues in the selected cluster. To
* select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input/output) COMPLEX array, dimension (LDT,N)
* On entry, the upper triangular matrix T.
* On exit, T is overwritten by the reordered matrix T, with the
* selected eigenvalues as the leading diagonal elements.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* Q (input/output) COMPLEX array, dimension (LDQ,N)
* On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
* On exit, if COMPQ = 'V', Q has been postmultiplied by the
* unitary transformation matrix which reorders T; the leading M
* columns of Q form an orthonormal basis for the specified
* invariant subspace.
* If COMPQ = 'N', Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q.
* LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
*
* W (output) COMPLEX array, dimension (N)
* The reordered eigenvalues of T, in the same order as they
* appear on the diagonal of T.
*
* M (output) INTEGER
* The dimension of the specified invariant subspace.
* 0 <= M <= N.
*
* S (output) REAL
* If JOB = 'E' or 'B', S is a lower bound on the reciprocal
* condition number for the selected cluster of eigenvalues.
* S cannot underestimate the true reciprocal condition number
* by more than a factor of sqrt(N). If M = 0 or N, S = 1.
* If JOB = 'N' or 'V', S is not referenced.
*
* SEP (output) REAL
* If JOB = 'V' or 'B', SEP is the estimated reciprocal
* condition number of the specified invariant subspace. If
* M = 0 or N, SEP = norm(T).
* If JOB = 'N' or 'E', SEP is not referenced.
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK.
* If JOB = 'N', LWORK >= 1;
* if JOB = 'E', LWORK = max(1,M*(N-M));
* if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
*
* 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
*
* Further Details
* ===============
*
* CTRSEN first collects the selected eigenvalues by computing a unitary
* transformation Z to move them to the top left corner of T. In other
* words, the selected eigenvalues are the eigenvalues of T11 in:
*
* Z'*T*Z = ( T11 T12 ) n1
* ( 0 T22 ) n2
* n1 n2
*
* where N = n1+n2 and Z' means the conjugate transpose of Z. The first
* n1 columns of Z span the specified invariant subspace of T.
*
* If T has been obtained from the Schur factorization of a matrix
* A = Q*T*Q', then the reordered Schur factorization of A is given by
* A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
* corresponding invariant subspace of A.
*
* The reciprocal condition number of the average of the eigenvalues of
* T11 may be returned in S. S lies between 0 (very badly conditioned)
* and 1 (very well conditioned). It is computed as follows. First we
* compute R so that
*
* P = ( I R ) n1
* ( 0 0 ) n2
* n1 n2
*
* is the projector on the invariant subspace associated with T11.
* R is the solution of the Sylvester equation:
*
* T11*R - R*T22 = T12.
*
* Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
* the two-norm of M. Then S is computed as the lower bound
*
* (1 + F-norm(R)**2)**(-1/2)
*
* on the reciprocal of 2-norm(P), the true reciprocal condition number.
* S cannot underestimate 1 / 2-norm(P) by more than a factor of
* sqrt(N).
*
* An approximate error bound for the computed average of the
* eigenvalues of T11 is
*
* EPS * norm(T) / S
*
* where EPS is the machine precision.
*
* The reciprocal condition number of the right invariant subspace
* spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
* SEP is defined as the separation of T11 and T22:
*
* sep( T11, T22 ) = sigma-min( C )
*
* where sigma-min(C) is the smallest singular value of the
* n1*n2-by-n1*n2 matrix
*
* C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
*
* I(m) is an m by m identity matrix, and kprod denotes the Kronecker
* product. We estimate sigma-min(C) by the reciprocal of an estimate of
* the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
* cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
*
* When SEP is small, small changes in T can cause large changes in
* the invariant subspace. An approximate bound on the maximum angular
* error in the computed right invariant subspace is
*
* EPS * norm(T) / SEP
*
* =====================================================================
*
go to the page top
ctrsna
USAGE:
s, sep, m, info = NumRu::Lapack.ctrsna( job, howmny, select, t, vl, vr, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO )
* Purpose
* =======
*
* CTRSNA estimates reciprocal condition numbers for specified
* eigenvalues and/or right eigenvectors of a complex upper triangular
* matrix T (or of any matrix Q*T*Q**H with Q unitary).
*
* Arguments
* =========
*
* JOB (input) CHARACTER*1
* Specifies whether condition numbers are required for
* eigenvalues (S) or eigenvectors (SEP):
* = 'E': for eigenvalues only (S);
* = 'V': for eigenvectors only (SEP);
* = 'B': for both eigenvalues and eigenvectors (S and SEP).
*
* HOWMNY (input) CHARACTER*1
* = 'A': compute condition numbers for all eigenpairs;
* = 'S': compute condition numbers for selected eigenpairs
* specified by the array SELECT.
*
* SELECT (input) LOGICAL array, dimension (N)
* If HOWMNY = 'S', SELECT specifies the eigenpairs for which
* condition numbers are required. To select condition numbers
* for the j-th eigenpair, SELECT(j) must be set to .TRUE..
* If HOWMNY = 'A', SELECT is not referenced.
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input) COMPLEX array, dimension (LDT,N)
* The upper triangular matrix T.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* VL (input) COMPLEX array, dimension (LDVL,M)
* If JOB = 'E' or 'B', VL must contain left eigenvectors of T
* (or of any Q*T*Q**H with Q unitary), corresponding to the
* eigenpairs specified by HOWMNY and SELECT. The eigenvectors
* must be stored in consecutive columns of VL, as returned by
* CHSEIN or CTREVC.
* If JOB = 'V', VL is not referenced.
*
* LDVL (input) INTEGER
* The leading dimension of the array VL.
* LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
*
* VR (input) COMPLEX array, dimension (LDVR,M)
* If JOB = 'E' or 'B', VR must contain right eigenvectors of T
* (or of any Q*T*Q**H with Q unitary), corresponding to the
* eigenpairs specified by HOWMNY and SELECT. The eigenvectors
* must be stored in consecutive columns of VR, as returned by
* CHSEIN or CTREVC.
* If JOB = 'V', VR is not referenced.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR.
* LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
*
* S (output) REAL array, dimension (MM)
* If JOB = 'E' or 'B', the reciprocal condition numbers of the
* selected eigenvalues, stored in consecutive elements of the
* array. Thus S(j), SEP(j), and the j-th columns of VL and VR
* all correspond to the same eigenpair (but not in general the
* j-th eigenpair, unless all eigenpairs are selected).
* If JOB = 'V', S is not referenced.
*
* SEP (output) REAL array, dimension (MM)
* If JOB = 'V' or 'B', the estimated reciprocal condition
* numbers of the selected eigenvectors, stored in consecutive
* elements of the array.
* If JOB = 'E', SEP is not referenced.
*
* MM (input) INTEGER
* The number of elements in the arrays S (if JOB = 'E' or 'B')
* and/or SEP (if JOB = 'V' or 'B'). MM >= M.
*
* M (output) INTEGER
* The number of elements of the arrays S and/or SEP actually
* used to store the estimated condition numbers.
* If HOWMNY = 'A', M is set to N.
*
* WORK (workspace) COMPLEX array, dimension (LDWORK,N+6)
* If JOB = 'E', WORK is not referenced.
*
* LDWORK (input) INTEGER
* The leading dimension of the array WORK.
* LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
*
* RWORK (workspace) REAL array, dimension (N)
* If JOB = 'E', RWORK is not referenced.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* The reciprocal of the condition number of an eigenvalue lambda is
* defined as
*
* S(lambda) = |v'*u| / (norm(u)*norm(v))
*
* where u and v are the right and left eigenvectors of T corresponding
* to lambda; v' denotes the conjugate transpose of v, and norm(u)
* denotes the Euclidean norm. These reciprocal condition numbers always
* lie between zero (very badly conditioned) and one (very well
* conditioned). If n = 1, S(lambda) is defined to be 1.
*
* An approximate error bound for a computed eigenvalue W(i) is given by
*
* EPS * norm(T) / S(i)
*
* where EPS is the machine precision.
*
* The reciprocal of the condition number of the right eigenvector u
* corresponding to lambda is defined as follows. Suppose
*
* T = ( lambda c )
* ( 0 T22 )
*
* Then the reciprocal condition number is
*
* SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
*
* where sigma-min denotes the smallest singular value. We approximate
* the smallest singular value by the reciprocal of an estimate of the
* one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
* defined to be abs(T(1,1)).
*
* An approximate error bound for a computed right eigenvector VR(i)
* is given by
*
* EPS * norm(T) / SEP(i)
*
* =====================================================================
*
go to the page top
ctrsyl
USAGE:
scale, info, c = NumRu::Lapack.ctrsyl( trana, tranb, isgn, a, b, c, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO )
* Purpose
* =======
*
* CTRSYL solves the complex Sylvester matrix equation:
*
* op(A)*X + X*op(B) = scale*C or
* op(A)*X - X*op(B) = scale*C,
*
* where op(A) = A or A**H, and A and B are both upper triangular. A is
* M-by-M and B is N-by-N; the right hand side C and the solution X are
* M-by-N; and scale is an output scale factor, set <= 1 to avoid
* overflow in X.
*
* Arguments
* =========
*
* TRANA (input) CHARACTER*1
* Specifies the option op(A):
* = 'N': op(A) = A (No transpose)
* = 'C': op(A) = A**H (Conjugate transpose)
*
* TRANB (input) CHARACTER*1
* Specifies the option op(B):
* = 'N': op(B) = B (No transpose)
* = 'C': op(B) = B**H (Conjugate transpose)
*
* ISGN (input) INTEGER
* Specifies the sign in the equation:
* = +1: solve op(A)*X + X*op(B) = scale*C
* = -1: solve op(A)*X - X*op(B) = scale*C
*
* M (input) INTEGER
* The order of the matrix A, and the number of rows in the
* matrices X and C. M >= 0.
*
* N (input) INTEGER
* The order of the matrix B, and the number of columns in the
* matrices X and C. N >= 0.
*
* A (input) COMPLEX array, dimension (LDA,M)
* The upper triangular matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
*
* B (input) COMPLEX array, dimension (LDB,N)
* The upper triangular matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* C (input/output) COMPLEX array, dimension (LDC,N)
* On entry, the M-by-N right hand side matrix C.
* On exit, C is overwritten by the solution matrix X.
*
* LDC (input) INTEGER
* The leading dimension of the array C. LDC >= max(1,M)
*
* SCALE (output) REAL
* The scale factor, scale, set <= 1 to avoid overflow in X.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* = 1: A and B have common or very close eigenvalues; perturbed
* values were used to solve the equation (but the matrices
* A and B are unchanged).
*
* =====================================================================
*
go to the page top
ctrti2
USAGE:
info, a = NumRu::Lapack.ctrti2( uplo, diag, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRTI2( UPLO, DIAG, N, A, LDA, INFO )
* Purpose
* =======
*
* CTRTI2 computes the inverse of a complex upper or lower triangular
* matrix.
*
* This is the Level 2 BLAS version of the algorithm.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the matrix A is upper or lower triangular.
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* DIAG (input) CHARACTER*1
* Specifies whether or not the matrix A is unit triangular.
* = 'N': Non-unit triangular
* = 'U': Unit triangular
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the triangular matrix A. If UPLO = 'U', the
* leading n by n upper triangular part of the array A contains
* the upper triangular matrix, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading n by n lower triangular part of the array A contains
* the lower triangular matrix, and the strictly upper
* triangular part of A is not referenced. If DIAG = 'U', the
* diagonal elements of A are also not referenced and are
* assumed to be 1.
*
* On exit, the (triangular) inverse of the original matrix, in
* the same storage format.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* =====================================================================
*
go to the page top
ctrtri
USAGE:
info, a = NumRu::Lapack.ctrtri( uplo, diag, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRTRI( UPLO, DIAG, N, A, LDA, INFO )
* Purpose
* =======
*
* CTRTRI computes the inverse of a complex upper or lower triangular
* matrix A.
*
* This is the Level 3 BLAS version of the algorithm.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the triangular matrix A. If UPLO = 'U', the
* leading N-by-N upper triangular part of the array A contains
* the upper triangular matrix, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading N-by-N lower triangular part of the array A contains
* the lower triangular matrix, and the strictly upper
* triangular part of A is not referenced. If DIAG = 'U', the
* diagonal elements of A are also not referenced and are
* assumed to be 1.
* On exit, the (triangular) inverse of the original matrix, in
* the same storage format.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, A(i,i) is exactly zero. The triangular
* matrix is singular and its inverse can not be computed.
*
* =====================================================================
*
go to the page top
ctrtrs
USAGE:
info, b = NumRu::Lapack.ctrtrs( uplo, trans, diag, a, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO )
* Purpose
* =======
*
* CTRTRS solves a triangular system of the form
*
* A * X = B, A**T * X = B, or A**H * X = B,
*
* where A is a triangular matrix of order N, and B is an N-by-NRHS
* matrix. A check is made to verify that A is nonsingular.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* TRANS (input) CHARACTER*1
* 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)
*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
*
* 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 array, dimension (LDA,N)
* The triangular matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of the array A contains the upper
* triangular matrix, and the strictly lower triangular part of
* A is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of the array A contains the lower triangular
* matrix, and the strictly upper triangular part of A is not
* referenced. If DIAG = 'U', the diagonal elements of A are
* also not referenced and are assumed to be 1.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input/output) COMPLEX array, dimension (LDB,NRHS)
* On entry, the right hand side matrix B.
* On exit, if INFO = 0, 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
* > 0: if INFO = i, the i-th diagonal element of A is zero,
* indicating that the matrix is singular and the solutions
* X have not been computed.
*
* =====================================================================
*
go to the page top
ctrttf
USAGE:
arf, info = NumRu::Lapack.ctrttf( transr, uplo, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRTTF( TRANSR, UPLO, N, A, LDA, ARF, INFO )
* Purpose
* =======
*
* CTRTTF copies a triangular matrix A from standard full format (TR)
* to rectangular full packed format (TF) .
*
* Arguments
* =========
*
* TRANSR (input) CHARACTER*1
* = 'N': ARF in Normal mode is wanted;
* = 'C': ARF in Conjugate Transpose mode is wanted;
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) COMPLEX array, dimension ( LDA, N )
* On entry, the triangular matrix A. If UPLO = 'U', the
* leading N-by-N upper triangular part of the array A contains
* the upper triangular matrix, and the strictly lower
* triangular part of A is not referenced. If UPLO = 'L', the
* leading N-by-N lower triangular part of the array A contains
* the lower triangular matrix, and the strictly upper
* triangular part of A is not referenced.
*
* LDA (input) INTEGER
* The leading dimension of the matrix A. LDA >= max(1,N).
*
* ARF (output) COMPLEX*16 array, dimension ( N*(N+1)/2 ),
* On exit, the upper or lower triangular matrix A stored in
* RFP format. For a further discussion see Notes below.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* We first consider Standard Packed Format when N is even.
* We give an example where N = 6.
*
* AP is Upper AP is Lower
*
* 00 01 02 03 04 05 00
* 11 12 13 14 15 10 11
* 22 23 24 25 20 21 22
* 33 34 35 30 31 32 33
* 44 45 40 41 42 43 44
* 55 50 51 52 53 54 55
*
*
* Let TRANSR = `N'. RFP holds AP as follows:
* For UPLO = `U' the upper trapezoid A(0:5,0:2) consists of the last
* three columns of AP upper. The lower triangle A(4:6,0:2) consists of
* conjugate-transpose of the first three columns of AP upper.
* For UPLO = `L' the lower trapezoid A(1:6,0:2) consists of the first
* three columns of AP lower. The upper triangle A(0:2,0:2) consists of
* conjugate-transpose of the last three columns of AP lower.
* To denote conjugate we place -- above the element. This covers the
* case N even and TRANSR = `N'.
*
* RFP A RFP A
*
* -- -- --
* 03 04 05 33 43 53
* -- --
* 13 14 15 00 44 54
* --
* 23 24 25 10 11 55
*
* 33 34 35 20 21 22
* --
* 00 44 45 30 31 32
* -- --
* 01 11 55 40 41 42
* -- -- --
* 02 12 22 50 51 52
*
* Now let TRANSR = `C'. RFP A in both UPLO cases is just the conjugate-
* transpose of RFP A above. One therefore gets:
*
*
* RFP A RFP A
*
* -- -- -- -- -- -- -- -- -- --
* 03 13 23 33 00 01 02 33 00 10 20 30 40 50
* -- -- -- -- -- -- -- -- -- --
* 04 14 24 34 44 11 12 43 44 11 21 31 41 51
* -- -- -- -- -- -- -- -- -- --
* 05 15 25 35 45 55 22 53 54 55 22 32 42 52
*
*
* We next consider Standard Packed Format when N is odd.
* We give an example where N = 5.
*
* AP is Upper AP is Lower
*
* 00 01 02 03 04 00
* 11 12 13 14 10 11
* 22 23 24 20 21 22
* 33 34 30 31 32 33
* 44 40 41 42 43 44
*
*
* Let TRANSR = `N'. RFP holds AP as follows:
* For UPLO = `U' the upper trapezoid A(0:4,0:2) consists of the last
* three columns of AP upper. The lower triangle A(3:4,0:1) consists of
* conjugate-transpose of the first two columns of AP upper.
* For UPLO = `L' the lower trapezoid A(0:4,0:2) consists of the first
* three columns of AP lower. The upper triangle A(0:1,1:2) consists of
* conjugate-transpose of the last two columns of AP lower.
* To denote conjugate we place -- above the element. This covers the
* case N odd and TRANSR = `N'.
*
* RFP A RFP A
*
* -- --
* 02 03 04 00 33 43
* --
* 12 13 14 10 11 44
*
* 22 23 24 20 21 22
* --
* 00 33 34 30 31 32
* -- --
* 01 11 44 40 41 42
*
* Now let TRANSR = `C'. RFP A in both UPLO cases is just the conjugate-
* transpose of RFP A above. One therefore gets:
*
*
* RFP A RFP A
*
* -- -- -- -- -- -- -- -- --
* 02 12 22 00 01 00 10 20 30 40 50
* -- -- -- -- -- -- -- -- --
* 03 13 23 33 11 33 11 21 31 41 51
* -- -- -- -- -- -- -- -- --
* 04 14 24 34 44 43 44 22 32 42 52
*
* =====================================================================
*
go to the page top
ctrttp
USAGE:
ap, info = NumRu::Lapack.ctrttp( uplo, a, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTRTTP( UPLO, N, A, LDA, AP, INFO )
* Purpose
* =======
*
* CTRTTP copies a triangular matrix A from full format (TR) to standard
* packed format (TP).
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* N (input) INTEGER
* The order of the matrices AP and A. N >= 0.
*
* A (input) COMPLEX array, dimension (LDA,N)
* On entry, the triangular 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).
*
* AP (output) COMPLEX array, dimension ( N*(N+1)/2 ),
* On exit, the upper or lower triangular matrix A, packed
* columnwise in a linear array. The j-th column of A is stored
* in the array AP as follows:
* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=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