COMPLEX routines for (complex) Hermitian, packed storage matrix
chpcon
USAGE:
rcond, info = NumRu::Lapack.chpcon( uplo, ap, ipiv, anorm, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, INFO )
* Purpose
* =======
*
* CHPCON estimates the reciprocal of the condition number of a complex
* Hermitian packed matrix A using the factorization A = U*D*U**H or
* A = L*D*L**H computed by CHPTRF.
*
* 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**H;
* = 'L': Lower triangular, form is A = L*D*L**H.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* AP (input) COMPLEX array, dimension (N*(N+1)/2)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by CHPTRF, stored as a
* packed triangular matrix.
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by CHPTRF.
*
* ANORM (input) REAL
* The 1-norm of the original matrix A.
*
* RCOND (output) REAL
* 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 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
chpev
USAGE:
w, z, info, ap = NumRu::Lapack.chpev( jobz, uplo, ap, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, INFO )
* Purpose
* =======
*
* CHPEV computes all the eigenvalues and, optionally, eigenvectors of a
* complex Hermitian matrix in packed storage.
*
* Arguments
* =========
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, AP is overwritten by values generated during the
* reduction to tridiagonal form. If UPLO = 'U', the diagonal
* and first superdiagonal of the tridiagonal matrix T overwrite
* the corresponding elements of A, and if UPLO = 'L', the
* diagonal and first subdiagonal of T overwrite the
* corresponding elements of A.
*
* W (output) REAL array, dimension (N)
* If INFO = 0, the eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, N)
* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
* eigenvectors of the matrix A, with the i-th column of Z
* holding the eigenvector associated with W(i).
* If JOBZ = 'N', then Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace) COMPLEX array, dimension (max(1, 2*N-1))
*
* RWORK (workspace) REAL array, dimension (max(1, 3*N-2))
*
* INFO (output) INTEGER
* = 0: successful exit.
* < 0: if INFO = -i, the i-th argument had an illegal value.
* > 0: if INFO = i, the algorithm failed to converge; i
* off-diagonal elements of an intermediate tridiagonal
* form did not converge to zero.
*
* =====================================================================
*
go to the page top
chpevd
USAGE:
w, z, work, rwork, iwork, info, ap = NumRu::Lapack.chpevd( jobz, uplo, ap, [:lwork => lwork, :lrwork => lrwork, :liwork => liwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
* Purpose
* =======
*
* CHPEVD computes all the eigenvalues and, optionally, eigenvectors of
* a complex Hermitian matrix A in packed storage. If eigenvectors are
* desired, it uses a divide and conquer algorithm.
*
* The divide and conquer algorithm makes very mild assumptions about
* floating point arithmetic. It will work on machines with a guard
* digit in add/subtract, or on those binary machines without guard
* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
* Cray-2. It could conceivably fail on hexadecimal or decimal machines
* without guard digits, but we know of none.
*
* Arguments
* =========
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, AP is overwritten by values generated during the
* reduction to tridiagonal form. If UPLO = 'U', the diagonal
* and first superdiagonal of the tridiagonal matrix T overwrite
* the corresponding elements of A, and if UPLO = 'L', the
* diagonal and first subdiagonal of T overwrite the
* corresponding elements of A.
*
* W (output) REAL array, dimension (N)
* If INFO = 0, the eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, N)
* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
* eigenvectors of the matrix A, with the i-th column of Z
* holding the eigenvector associated with W(i).
* If JOBZ = 'N', then Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the required LWORK.
*
* LWORK (input) INTEGER
* The dimension of array WORK.
* If N <= 1, LWORK must be at least 1.
* If JOBZ = 'N' and N > 1, LWORK must be at least N.
* If JOBZ = 'V' and N > 1, LWORK must be at least 2*N.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the required sizes of the WORK, RWORK and
* IWORK arrays, returns these values as the first entries of
* the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK is issued by XERBLA.
*
* RWORK (workspace/output) REAL array, dimension (MAX(1,LRWORK))
* On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
*
* LRWORK (input) INTEGER
* The dimension of array RWORK.
* If N <= 1, LRWORK must be at least 1.
* If JOBZ = 'N' and N > 1, LRWORK must be at least N.
* If JOBZ = 'V' and N > 1, LRWORK must be at least
* 1 + 5*N + 2*N**2.
*
* If LRWORK = -1, then a workspace query is assumed; the
* routine only calculates the required sizes of the WORK, RWORK
* and IWORK arrays, returns these values as the first entries
* of the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK is issued by XERBLA.
*
* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
* On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
*
* LIWORK (input) INTEGER
* The dimension of array IWORK.
* If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
* If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
*
* If LIWORK = -1, then a workspace query is assumed; the
* routine only calculates the required sizes of the WORK, RWORK
* and IWORK arrays, returns these values as the first entries
* of the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK 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, the algorithm failed to converge; i
* off-diagonal elements of an intermediate tridiagonal
* form did not converge to zero.
*
* =====================================================================
*
go to the page top
chpevx
USAGE:
m, w, z, ifail, info, ap = NumRu::Lapack.chpevx( jobz, range, uplo, ap, vl, vu, il, iu, abstol, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO )
* Purpose
* =======
*
* CHPEVX computes selected eigenvalues and, optionally, eigenvectors
* of a complex Hermitian matrix A in packed storage.
* Eigenvalues/vectors can be selected by specifying either a range of
* values or a range of indices for the desired eigenvalues.
*
* Arguments
* =========
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* RANGE (input) CHARACTER*1
* = 'A': all eigenvalues will be found;
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found;
* = 'I': the IL-th through IU-th eigenvalues will be found.
*
* 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, AP is overwritten by values generated during the
* reduction to tridiagonal form. If UPLO = 'U', the diagonal
* and first superdiagonal of the tridiagonal matrix T overwrite
* the corresponding elements of A, and if UPLO = 'L', the
* diagonal and first subdiagonal of T overwrite the
* corresponding elements of A.
*
* VL (input) REAL
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
* Not referenced if RANGE = 'A' or 'V'.
*
* ABSTOL (input) REAL
* The absolute error tolerance for the eigenvalues.
* An approximate eigenvalue is accepted as converged
* when it is determined to lie in an interval [a,b]
* of width less than or equal to
*
* ABSTOL + EPS * max( |a|,|b| ) ,
*
* where EPS is the machine precision. If ABSTOL is less than
* or equal to zero, then EPS*|T| will be used in its place,
* where |T| is the 1-norm of the tridiagonal matrix obtained
* by reducing AP to tridiagonal form.
*
* Eigenvalues will be computed most accurately when ABSTOL is
* set to twice the underflow threshold 2*SLAMCH('S'), not zero.
* If this routine returns with INFO>0, indicating that some
* eigenvectors did not converge, try setting ABSTOL to
* 2*SLAMCH('S').
*
* See "Computing Small Singular Values of Bidiagonal Matrices
* with Guaranteed High Relative Accuracy," by Demmel and
* Kahan, LAPACK Working Note #3.
*
* M (output) INTEGER
* The total number of eigenvalues found. 0 <= M <= N.
* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*
* W (output) REAL array, dimension (N)
* If INFO = 0, the selected eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, max(1,M))
* If JOBZ = 'V', then if INFO = 0, the first M columns of Z
* contain the orthonormal eigenvectors of the matrix A
* corresponding to the selected eigenvalues, with the i-th
* column of Z holding the eigenvector associated with W(i).
* If an eigenvector fails to converge, then that column of Z
* contains the latest approximation to the eigenvector, and
* the index of the eigenvector is returned in IFAIL.
* If JOBZ = 'N', then Z is not referenced.
* Note: the user must ensure that at least max(1,M) columns are
* supplied in the array Z; if RANGE = 'V', the exact value of M
* is not known in advance and an upper bound must be used.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* RWORK (workspace) REAL array, dimension (7*N)
*
* IWORK (workspace) INTEGER array, dimension (5*N)
*
* IFAIL (output) INTEGER array, dimension (N)
* If JOBZ = 'V', then if INFO = 0, the first M elements of
* IFAIL are zero. If INFO > 0, then IFAIL contains the
* indices of the eigenvectors that failed to converge.
* If JOBZ = 'N', then IFAIL is not referenced.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, then i eigenvectors failed to converge.
* Their indices are stored in array IFAIL.
*
* =====================================================================
*
go to the page top
chpgst
USAGE:
info, ap = NumRu::Lapack.chpgst( itype, uplo, n, ap, bp, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPGST( ITYPE, UPLO, N, AP, BP, INFO )
* Purpose
* =======
*
* CHPGST reduces a complex Hermitian-definite generalized
* eigenproblem to standard form, using packed storage.
*
* If ITYPE = 1, the problem is A*x = lambda*B*x,
* and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
*
* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
* B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
*
* B must have been previously factorized as U**H*U or L*L**H by CPPTRF.
*
* Arguments
* =========
*
* ITYPE (input) INTEGER
* = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
* = 2 or 3: compute U*A*U**H or L**H*A*L.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored and B is factored as
* U**H*U;
* = 'L': Lower triangle of A is stored and B is factored as
* L*L**H.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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.
*
* On exit, if INFO = 0, the transformed matrix, stored in the
* same format as A.
*
* BP (input) COMPLEX array, dimension (N*(N+1)/2)
* The triangular factor from the Cholesky factorization of B,
* stored in the same format as A, as returned by CPPTRF.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
go to the page top
chpgv
USAGE:
w, z, info, ap, bp = NumRu::Lapack.chpgv( itype, jobz, uplo, ap, bp, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPGV( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, RWORK, INFO )
* Purpose
* =======
*
* CHPGV computes all the eigenvalues and, optionally, the eigenvectors
* of a complex generalized Hermitian-definite eigenproblem, of the form
* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
* Here A and B are assumed to be Hermitian, stored in packed format,
* and B is also positive definite.
*
* Arguments
* =========
*
* ITYPE (input) INTEGER
* Specifies the problem type to be solved:
* = 1: A*x = (lambda)*B*x
* = 2: A*B*x = (lambda)*x
* = 3: B*A*x = (lambda)*x
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangles of A and B are stored;
* = 'L': Lower triangles of A and B are stored.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, the contents of AP are destroyed.
*
* BP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian matrix
* B, packed columnwise in a linear array. The j-th column of B
* is stored in the array BP as follows:
* if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
* if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
*
* On exit, the triangular factor U or L from the Cholesky
* factorization B = U**H*U or B = L*L**H, in the same storage
* format as B.
*
* W (output) REAL array, dimension (N)
* If INFO = 0, the eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, N)
* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
* eigenvectors. The eigenvectors are normalized as follows:
* if ITYPE = 1 or 2, Z**H*B*Z = I;
* if ITYPE = 3, Z**H*inv(B)*Z = I.
* If JOBZ = 'N', then Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace) COMPLEX array, dimension (max(1, 2*N-1))
*
* RWORK (workspace) REAL array, dimension (max(1, 3*N-2))
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: CPPTRF or CHPEV returned an error code:
* <= N: if INFO = i, CHPEV failed to converge;
* i off-diagonal elements of an intermediate
* tridiagonal form did not convergeto zero;
* > N: if INFO = N + i, for 1 <= i <= n, then the leading
* minor of order i of B is not positive definite.
* The factorization of B could not be completed and
* no eigenvalues or eigenvectors were computed.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL UPPER, WANTZ
CHARACTER TRANS
INTEGER J, NEIG
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL CHPEV, CHPGST, CPPTRF, CTPMV, CTPSV, XERBLA
* ..
go to the page top
chpgvd
USAGE:
w, z, iwork, info, ap, bp = NumRu::Lapack.chpgvd( itype, jobz, uplo, ap, bp, [:lwork => lwork, :lrwork => lrwork, :liwork => liwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPGVD( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
* Purpose
* =======
*
* CHPGVD computes all the eigenvalues and, optionally, the eigenvectors
* of a complex generalized Hermitian-definite eigenproblem, of the form
* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
* B are assumed to be Hermitian, stored in packed format, and B is also
* positive definite.
* If eigenvectors are desired, it uses a divide and conquer algorithm.
*
* The divide and conquer algorithm makes very mild assumptions about
* floating point arithmetic. It will work on machines with a guard
* digit in add/subtract, or on those binary machines without guard
* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
* Cray-2. It could conceivably fail on hexadecimal or decimal machines
* without guard digits, but we know of none.
*
* Arguments
* =========
*
* ITYPE (input) INTEGER
* Specifies the problem type to be solved:
* = 1: A*x = (lambda)*B*x
* = 2: A*B*x = (lambda)*x
* = 3: B*A*x = (lambda)*x
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangles of A and B are stored;
* = 'L': Lower triangles of A and B are stored.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, the contents of AP are destroyed.
*
* BP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian matrix
* B, packed columnwise in a linear array. The j-th column of B
* is stored in the array BP as follows:
* if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
* if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
*
* On exit, the triangular factor U or L from the Cholesky
* factorization B = U**H*U or B = L*L**H, in the same storage
* format as B.
*
* W (output) REAL array, dimension (N)
* If INFO = 0, the eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, N)
* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
* eigenvectors. The eigenvectors are normalized as follows:
* if ITYPE = 1 or 2, Z**H*B*Z = I;
* if ITYPE = 3, Z**H*inv(B)*Z = I.
* If JOBZ = 'N', then Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the required LWORK.
*
* LWORK (input) INTEGER
* The dimension of array WORK.
* If N <= 1, LWORK >= 1.
* If JOBZ = 'N' and N > 1, LWORK >= N.
* If JOBZ = 'V' and N > 1, LWORK >= 2*N.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the required sizes of the WORK, RWORK and
* IWORK arrays, returns these values as the first entries of
* the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK is issued by XERBLA.
*
* RWORK (workspace) REAL array, dimension (MAX(1,LRWORK))
* On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
*
* LRWORK (input) INTEGER
* The dimension of array RWORK.
* If N <= 1, LRWORK >= 1.
* If JOBZ = 'N' and N > 1, LRWORK >= N.
* If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
*
* If LRWORK = -1, then a workspace query is assumed; the
* routine only calculates the required sizes of the WORK, RWORK
* and IWORK arrays, returns these values as the first entries
* of the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK is issued by XERBLA.
*
* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
* On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
*
* LIWORK (input) INTEGER
* The dimension of array IWORK.
* If JOBZ = 'N' or N <= 1, LIWORK >= 1.
* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
*
* If LIWORK = -1, then a workspace query is assumed; the
* routine only calculates the required sizes of the WORK, RWORK
* and IWORK arrays, returns these values as the first entries
* of the WORK, RWORK and IWORK arrays, and no error message
* related to LWORK or LRWORK or LIWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: CPPTRF or CHPEVD returned an error code:
* <= N: if INFO = i, CHPEVD failed to converge;
* i off-diagonal elements of an intermediate
* tridiagonal form did not convergeto zero;
* > N: if INFO = N + i, for 1 <= i <= n, then the leading
* minor of order i of B is not positive definite.
* The factorization of B could not be completed and
* no eigenvalues or eigenvectors were computed.
*
* Further Details
* ===============
*
* Based on contributions by
* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL LQUERY, UPPER, WANTZ
CHARACTER TRANS
INTEGER J, LIWMIN, LRWMIN, LWMIN, NEIG
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL CHPEVD, CHPGST, CPPTRF, CTPMV, CTPSV, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, REAL
* ..
go to the page top
chpgvx
USAGE:
m, w, z, ifail, info, ap, bp = NumRu::Lapack.chpgvx( itype, jobz, range, uplo, ap, bp, vl, vu, il, iu, abstol, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO )
* Purpose
* =======
*
* CHPGVX computes selected eigenvalues and, optionally, eigenvectors
* of a complex generalized Hermitian-definite eigenproblem, of the form
* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
* B are assumed to be Hermitian, stored in packed format, and B is also
* positive definite. Eigenvalues and eigenvectors can be selected by
* specifying either a range of values or a range of indices for the
* desired eigenvalues.
*
* Arguments
* =========
*
* ITYPE (input) INTEGER
* Specifies the problem type to be solved:
* = 1: A*x = (lambda)*B*x
* = 2: A*B*x = (lambda)*x
* = 3: B*A*x = (lambda)*x
*
* JOBZ (input) CHARACTER*1
* = 'N': Compute eigenvalues only;
* = 'V': Compute eigenvalues and eigenvectors.
*
* RANGE (input) CHARACTER*1
* = 'A': all eigenvalues will be found;
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found;
* = 'I': the IL-th through IU-th eigenvalues will be found.
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangles of A and B are stored;
* = 'L': Lower triangles of A and B are stored.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, the contents of AP are destroyed.
*
* BP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian matrix
* B, packed columnwise in a linear array. The j-th column of B
* is stored in the array BP as follows:
* if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
* if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
*
* On exit, the triangular factor U or L from the Cholesky
* factorization B = U**H*U or B = L*L**H, in the same storage
* format as B.
*
* VL (input) REAL
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
* Not referenced if RANGE = 'A' or 'V'.
*
* ABSTOL (input) REAL
* The absolute error tolerance for the eigenvalues.
* An approximate eigenvalue is accepted as converged
* when it is determined to lie in an interval [a,b]
* of width less than or equal to
*
* ABSTOL + EPS * max( |a|,|b| ) ,
*
* where EPS is the machine precision. If ABSTOL is less than
* or equal to zero, then EPS*|T| will be used in its place,
* where |T| is the 1-norm of the tridiagonal matrix obtained
* by reducing AP to tridiagonal form.
*
* Eigenvalues will be computed most accurately when ABSTOL is
* set to twice the underflow threshold 2*SLAMCH('S'), not zero.
* If this routine returns with INFO>0, indicating that some
* eigenvectors did not converge, try setting ABSTOL to
* 2*SLAMCH('S').
*
* M (output) INTEGER
* The total number of eigenvalues found. 0 <= M <= N.
* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*
* W (output) REAL array, dimension (N)
* On normal exit, the first M elements contain the selected
* eigenvalues in ascending order.
*
* Z (output) COMPLEX array, dimension (LDZ, N)
* If JOBZ = 'N', then Z is not referenced.
* If JOBZ = 'V', then if INFO = 0, the first M columns of Z
* contain the orthonormal eigenvectors of the matrix A
* corresponding to the selected eigenvalues, with the i-th
* column of Z holding the eigenvector associated with W(i).
* The eigenvectors are normalized as follows:
* if ITYPE = 1 or 2, Z**H*B*Z = I;
* if ITYPE = 3, Z**H*inv(B)*Z = I.
*
* If an eigenvector fails to converge, then that column of Z
* contains the latest approximation to the eigenvector, and the
* index of the eigenvector is returned in IFAIL.
* Note: the user must ensure that at least max(1,M) columns are
* supplied in the array Z; if RANGE = 'V', the exact value of M
* is not known in advance and an upper bound must be used.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1, and if
* JOBZ = 'V', LDZ >= max(1,N).
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* RWORK (workspace) REAL array, dimension (7*N)
*
* IWORK (workspace) INTEGER array, dimension (5*N)
*
* IFAIL (output) INTEGER array, dimension (N)
* If JOBZ = 'V', then if INFO = 0, the first M elements of
* IFAIL are zero. If INFO > 0, then IFAIL contains the
* indices of the eigenvectors that failed to converge.
* If JOBZ = 'N', then IFAIL is not referenced.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: CPPTRF or CHPEVX returned an error code:
* <= N: if INFO = i, CHPEVX failed to converge;
* i eigenvectors failed to converge. Their indices
* are stored in array IFAIL.
* > N: if INFO = N + i, for 1 <= i <= n, then the leading
* minor of order i of B is not positive definite.
* The factorization of B could not be completed and
* no eigenvalues or eigenvectors were computed.
*
* Further Details
* ===============
*
* Based on contributions by
* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
CHARACTER TRANS
INTEGER J
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL CHPEVX, CHPGST, CPPTRF, CTPMV, CTPSV, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MIN
* ..
go to the page top
chprfs
USAGE:
ferr, berr, info, x = NumRu::Lapack.chprfs( uplo, ap, afp, ipiv, b, x, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
* Purpose
* =======
*
* CHPRFS improves the computed solution to a system of linear
* equations when the coefficient matrix is Hermitian indefinite
* and packed, 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.
*
* AP (input) COMPLEX array, dimension (N*(N+1)/2)
* The upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
*
* AFP (input) COMPLEX array, dimension (N*(N+1)/2)
* The factored form of the matrix A. AFP contains the block
* diagonal matrix D and the multipliers used to obtain the
* factor U or L from the factorization A = U*D*U**H or
* A = L*D*L**H as computed by CHPTRF, stored as a packed
* triangular matrix.
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by CHPTRF.
*
* 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/output) COMPLEX array, dimension (LDX,NRHS)
* On entry, the solution matrix X, as computed by CHPTRS.
* On exit, the improved 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
*
* Internal Parameters
* ===================
*
* ITMAX is the maximum number of steps of iterative refinement.
*
* =====================================================================
*
go to the page top
chpsv
USAGE:
ipiv, info, ap, b = NumRu::Lapack.chpsv( uplo, ap, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPSV( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
* Purpose
* =======
*
* CHPSV computes the solution to a complex system of linear equations
* A * X = B,
* where A is an N-by-N Hermitian matrix stored in packed format and X
* and B are N-by-NRHS matrices.
*
* The diagonal pivoting method is used to factor A as
* A = U * D * U**H, if UPLO = 'U', or
* A = L * D * L**H, if UPLO = 'L',
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, D is Hermitian 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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.
* See below for further details.
*
* On exit, the block diagonal matrix D and the multipliers used
* to obtain the factor U or L from the factorization
* A = U*D*U**H or A = L*D*L**H as computed by CHPTRF, stored as
* a packed triangular matrix in the same storage format as A.
*
* IPIV (output) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D, as
* determined by CHPTRF. 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 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).
*
* 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.
*
* Further Details
* ===============
*
* The packed storage scheme is illustrated by the following example
* when N = 4, UPLO = 'U':
*
* Two-dimensional storage of the Hermitian matrix A:
*
* a11 a12 a13 a14
* a22 a23 a24
* a33 a34 (aij = conjg(aji))
* a44
*
* Packed storage of the upper triangle of A:
*
* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL CHPTRF, CHPTRS, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
go to the page top
chpsvx
USAGE:
x, rcond, ferr, berr, info, afp, ipiv = NumRu::Lapack.chpsvx( fact, uplo, ap, afp, ipiv, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
* Purpose
* =======
*
* CHPSVX uses the diagonal pivoting factorization A = U*D*U**H or
* A = L*D*L**H to compute the solution to a complex system of linear
* equations A * X = B, where A is an N-by-N Hermitian matrix stored
* in packed format 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 as
* A = U * D * U**H, if UPLO = 'U', or
* A = L * D * L**H, if UPLO = 'L',
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices and D is Hermitian 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, AFP and IPIV contain the factored form of
* A. AFP and IPIV will not be modified.
* = 'N': The matrix A will be copied to AFP 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.
*
* AP (input) COMPLEX array, dimension (N*(N+1)/2)
* The upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
* See below for further details.
*
* AFP (input or output) COMPLEX array, dimension (N*(N+1)/2)
* If FACT = 'F', then AFP 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**H or A = L*D*L**H as computed by CHPTRF, stored as
* a packed triangular matrix in the same storage format as A.
*
* If FACT = 'N', then AFP is an output argument and on exit
* contains the block diagonal matrix D and the multipliers used
* to obtain the factor U or L from the factorization
* A = U*D*U**H or A = L*D*L**H as computed by CHPTRF, stored as
* a packed triangular matrix in the same storage format as A.
*
* 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 CHPTRF.
* 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 CHPTRF.
*
* B (input) COMPLEX 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 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) REAL
* 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) 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
* > 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.
*
* Further Details
* ===============
*
* The packed storage scheme is illustrated by the following example
* when N = 4, UPLO = 'U':
*
* Two-dimensional storage of the Hermitian matrix A:
*
* a11 a12 a13 a14
* a22 a23 a24
* a33 a34 (aij = conjg(aji))
* a44
*
* Packed storage of the upper triangle of A:
*
* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
*
* =====================================================================
*
go to the page top
chptrd
USAGE:
d, e, tau, info, ap = NumRu::Lapack.chptrd( uplo, ap, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPTRD( UPLO, N, AP, D, E, TAU, INFO )
* Purpose
* =======
*
* CHPTRD reduces a complex Hermitian matrix A stored in packed form to
* real symmetric tridiagonal form T by a unitary similarity
* transformation: Q**H * A * Q = T.
*
* 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
* On exit, if UPLO = 'U', the diagonal and first superdiagonal
* of A are overwritten by the corresponding elements of the
* tridiagonal matrix T, and the elements above the first
* superdiagonal, with the array TAU, represent the unitary
* matrix Q as a product of elementary reflectors; if UPLO
* = 'L', the diagonal and first subdiagonal of A are over-
* written by the corresponding elements of the tridiagonal
* matrix T, and the elements below the first subdiagonal, with
* the array TAU, represent the unitary matrix Q as a product
* of elementary reflectors. See Further Details.
*
* D (output) REAL array, dimension (N)
* The diagonal elements of the tridiagonal matrix T:
* D(i) = A(i,i).
*
* E (output) REAL array, dimension (N-1)
* The off-diagonal elements of the tridiagonal matrix T:
* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
*
* TAU (output) COMPLEX array, dimension (N-1)
* The scalar factors of the elementary reflectors (see Further
* Details).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* If UPLO = 'U', the matrix Q is represented as a product of elementary
* reflectors
*
* Q = H(n-1) . . . H(2) H(1).
*
* Each H(i) has the form
*
* H(i) = I - tau * v * v'
*
* where tau is a complex scalar, and v is a complex vector with
* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
* overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
*
* If UPLO = 'L', the matrix Q is represented as a product of elementary
* reflectors
*
* Q = H(1) H(2) . . . H(n-1).
*
* Each H(i) has the form
*
* H(i) = I - tau * v * v'
*
* where tau is a complex scalar, and v is a complex vector with
* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
* overwriting A(i+2:n,i), and tau is stored in TAU(i).
*
* =====================================================================
*
go to the page top
chptrf
USAGE:
ipiv, info, ap = NumRu::Lapack.chptrf( uplo, ap, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO )
* Purpose
* =======
*
* CHPTRF computes the factorization of a complex Hermitian packed
* matrix A using the Bunch-Kaufman diagonal pivoting method:
*
* A = U*D*U**H or A = L*D*L**H
*
* where U (or L) is a product of permutation and unit upper (lower)
* triangular matrices, and D is Hermitian and block diagonal with
* 1-by-1 and 2-by-2 diagonal blocks.
*
* 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.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the Hermitian 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.
*
* On exit, the block diagonal matrix D and the multipliers used
* to obtain the factor U or L, stored as a packed triangular
* matrix overwriting A (see below for further details).
*
* 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 = -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
* ===============
*
* 5-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
chptri
USAGE:
info, ap = NumRu::Lapack.chptri( uplo, ap, ipiv, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
* Purpose
* =======
*
* CHPTRI computes the inverse of a complex Hermitian indefinite matrix
* A in packed storage using the factorization A = U*D*U**H or
* A = L*D*L**H computed by CHPTRF.
*
* 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**H;
* = 'L': Lower triangular, form is A = L*D*L**H.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)
* On entry, the block diagonal matrix D and the multipliers
* used to obtain the factor U or L as computed by CHPTRF,
* stored as a packed triangular matrix.
*
* On exit, if INFO = 0, the (Hermitian) inverse of the original
* matrix, stored as a packed triangular matrix. The j-th column
* of inv(A) is stored in the array AP as follows:
* if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
* if UPLO = 'L',
* AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by CHPTRF.
*
* WORK (workspace) COMPLEX 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, D(i,i) = 0; the matrix is singular and its
* inverse could not be computed.
*
* =====================================================================
*
go to the page top
chptrs
USAGE:
info, b = NumRu::Lapack.chptrs( uplo, ap, ipiv, b, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CHPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
* Purpose
* =======
*
* CHPTRS solves a system of linear equations A*X = B with a complex
* Hermitian matrix A stored in packed format using the factorization
* A = U*D*U**H or A = L*D*L**H computed by CHPTRF.
*
* 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**H;
* = 'L': Lower triangular, form is A = L*D*L**H.
*
* 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.
*
* AP (input) COMPLEX array, dimension (N*(N+1)/2)
* The block diagonal matrix D and the multipliers used to
* obtain the factor U or L as computed by CHPTRF, stored as a
* packed triangular matrix.
*
* IPIV (input) INTEGER array, dimension (N)
* Details of the interchanges and the block structure of D
* as determined by CHPTRF.
*
* B (input/output) 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).
*
* 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