/* Subroutine */ int cgbrfs_(char *trans, integer *n, integer *kl, integer *
ku, integer *nrhs, complex *ab, integer *ldab, complex *afb, integer *
ldafb, integer *ipiv, complex *b, integer *ldb, complex *x, integer *
ldx, real *ferr, real *berr, complex *work, real *rwork, integer *
info)
{
/* System generated locals */
integer ab_dim1, ab_offset, afb_dim1, afb_offset, b_dim1, b_offset,
x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;
real r__1, r__2, r__3, r__4;
complex q__1;
/* Builtin functions */
double r_imag(complex *);
/* Local variables */
integer i__, j, k;
real s;
integer kk;
real xk;
integer nz;
real eps;
integer kase;
real safe1, safe2;
extern /* Subroutine */ int cgbmv_(char *, integer *, integer *, integer *
, integer *, complex *, complex *, integer *, complex *, integer *
, complex *, complex *, integer *);
extern logical lsame_(char *, char *);
integer isave[3];
extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
complex *, integer *), caxpy_(integer *, complex *, complex *,
integer *, complex *, integer *);
integer count;
extern /* Subroutine */ int clacn2_(integer *, complex *, complex *, real
*, integer *, integer *);
extern doublereal slamch_(char *);
real safmin;
extern /* Subroutine */ int xerbla_(char *, integer *), cgbtrs_(
char *, integer *, integer *, integer *, integer *, complex *,
integer *, integer *, complex *, integer *, integer *);
logical notran;
char transn[1], transt[1];
real lstres;
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CGBRFS improves the computed solution to a system of linear */
/* equations when the coefficient matrix is banded, and provides */
/* error bounds and backward error estimates for the solution. */
/* Arguments */
/* ========= */
/* 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) */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* KL (input) INTEGER */
/* The number of subdiagonals within the band of A. KL >= 0. */
/* KU (input) INTEGER */
/* The number of superdiagonals within the band of A. KU >= 0. */
/* NRHS (input) INTEGER */
/* The number of right hand sides, i.e., the number of columns */
/* of the matrices B and X. NRHS >= 0. */
/* AB (input) COMPLEX array, dimension (LDAB,N) */
/* The original band matrix A, stored in rows 1 to KL+KU+1. */
/* The j-th column of A is stored in the j-th column of the */
/* array AB as follows: */
/* AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl). */
/* LDAB (input) INTEGER */
/* The leading dimension of the array AB. LDAB >= KL+KU+1. */
/* AFB (input) COMPLEX array, dimension (LDAFB,N) */
/* Details of the LU factorization of the band matrix A, as */
/* computed by CGBTRF. U is stored as an upper triangular band */
/* matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and */
/* the multipliers used during the factorization are stored in */
/* rows KL+KU+2 to 2*KL+KU+1. */
//.........这里部分代码省略.........
/* DECK CSIDI */
/* Subroutine */ int csidi_(complex *a, integer *lda, integer *n, integer *
kpvt, complex *det, complex *work, integer *job)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
real r__1, r__2;
complex q__1, q__2, q__3;
/* Local variables */
static complex d__;
static integer j, k;
static complex t, ak;
static integer jb, ks, km1;
static real ten;
static complex akp1, temp, akkp1;
static logical nodet;
extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
complex *, integer *);
extern /* Complex */ void cdotu_(complex *, integer *, complex *, integer
*, complex *, integer *);
extern /* Subroutine */ int cswap_(integer *, complex *, integer *,
complex *, integer *), caxpy_(integer *, complex *, complex *,
integer *, complex *, integer *);
static integer kstep;
static logical noinv;
/* ***BEGIN PROLOGUE CSIDI */
/* ***PURPOSE Compute the determinant and inverse of a complex symmetric */
/* matrix using the factors from CSIFA. */
/* ***LIBRARY SLATEC (LINPACK) */
/* ***CATEGORY D2C1, D3C1 */
/* ***TYPE COMPLEX (SSIDI-S, DSIDI-D, CHIDI-C, CSIDI-C) */
/* ***KEYWORDS DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX, */
/* SYMMETRIC */
/* ***AUTHOR Bunch, J., (UCSD) */
/* ***DESCRIPTION */
/* CSIDI computes the determinant and inverse */
/* of a complex symmetric matrix using the factors from CSIFA. */
/* On Entry */
/* A COMPLEX(LDA,N) */
/* the output from CSIFA. */
/* LDA INTEGER */
/* the leading dimension of the array A . */
/* N INTEGER */
/* the order of the matrix A . */
/* KVPT INTEGER(N) */
/* the pivot vector from CSIFA. */
/* WORK COMPLEX(N) */
/* work vector. Contents destroyed. */
/* JOB INTEGER */
/* JOB has the decimal expansion AB where */
/* If B .NE. 0, the inverse is computed, */
/* If A .NE. 0, the determinant is computed, */
/* For example, JOB = 11 gives both. */
/* On Return */
/* Variables not requested by JOB are not used. */
/* A contains the upper triangle of the inverse of */
/* the original matrix. The strict lower triangle */
/* is never referenced. */
/* DET COMPLEX(2) */
/* determinant of original matrix. */
/* Determinant = DET(1) * 10.0**DET(2) */
/* with 1.0 .LE. ABS(DET(1)) .LT. 10.0 */
/* or DET(1) = 0.0. */
/* Error Condition */
/* A division by zero may occur if the inverse is requested */
/* and CSICO has set RCOND .EQ. 0.0 */
/* or CSIFA has set INFO .NE. 0 . */
/* ***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. */
/* Stewart, LINPACK Users' Guide, SIAM, 1979. */
/* ***ROUTINES CALLED CAXPY, CCOPY, CDOTU, CSWAP */
/* ***REVISION HISTORY (YYMMDD) */
/* 780814 DATE WRITTEN */
/* 890531 Changed all specific intrinsics to generic. (WRB) */
/* 890831 Modified array declarations. (WRB) */
/* 891107 Corrected category and modified routine equivalence */
/* list. (WRB) */
/* 891107 REVISION DATE from Version 3.2 */
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
/* 900326 Removed duplicate information from DESCRIPTION section. */
/* (WRB) */
/* 920501 Reformatted the REFERENCES section. (WRB) */
/* ***END PROLOGUE CSIDI */
//.........这里部分代码省略.........
/* Subroutine */ int clarz_(char *side, integer *m, integer *n, integer *l,
complex *v, integer *incv, complex *tau, complex *c__, integer *ldc,
complex *work)
{
/* System generated locals */
integer c_dim1, c_offset;
complex q__1;
/* Local variables */
extern /* Subroutine */ int cgerc_(integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, integer *),
cgemv_(char *, integer *, integer *, complex *, complex *,
integer *, complex *, integer *, complex *, complex *, integer *);
extern logical lsame_(char *, char *);
extern /* Subroutine */ int cgeru_(integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, integer *),
ccopy_(integer *, complex *, integer *, complex *, integer *),
caxpy_(integer *, complex *, complex *, integer *, complex *,
integer *), clacgv_(integer *, complex *, integer *);
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLARZ applies a complex elementary reflector H to a complex */
/* M-by-N matrix C, from either the left or the right. H is represented */
/* in the form */
/* H = I - tau * v * v' */
/* where tau is a complex scalar and v is a complex vector. */
/* If tau = 0, then H is taken to be the unit matrix. */
/* To apply H' (the conjugate transpose of H), supply conjg(tau) instead */
/* tau. */
/* H is a product of k elementary reflectors as returned by CTZRZF. */
/* Arguments */
/* ========= */
/* SIDE (input) CHARACTER*1 */
/* = 'L': form H * C */
/* = 'R': form C * H */
/* M (input) INTEGER */
/* The number of rows of the matrix C. */
/* N (input) INTEGER */
/* The number of columns of the matrix C. */
/* L (input) INTEGER */
/* The number of entries of the vector V containing */
/* the meaningful part of the Householder vectors. */
/* If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. */
/* V (input) COMPLEX array, dimension (1+(L-1)*abs(INCV)) */
/* The vector v in the representation of H as returned by */
/* CTZRZF. V is not used if TAU = 0. */
/* INCV (input) INTEGER */
/* The increment between elements of v. INCV <> 0. */
/* TAU (input) COMPLEX */
/* The value tau in the representation of H. */
/* C (input/output) COMPLEX array, dimension (LDC,N) */
/* On entry, the M-by-N matrix C. */
/* On exit, C is overwritten by the matrix H * C if SIDE = 'L', */
/* or C * H if SIDE = 'R'. */
/* LDC (input) INTEGER */
/* The leading dimension of the array C. LDC >= max(1,M). */
/* WORK (workspace) COMPLEX array, dimension */
/* (N) if SIDE = 'L' */
/* or (M) if SIDE = 'R' */
/* Further Details */
/* =============== */
/* Based on contributions by */
/* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
//.........这里部分代码省略.........
/* Subroutine */ int clahrd_(integer *n, integer *k, integer *nb, complex *a,
integer *lda, complex *tau, complex *t, integer *ldt, complex *y,
integer *ldy)
{
/* System generated locals */
integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2,
i__3;
complex q__1;
/* Local variables */
integer i__;
complex ei;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *), cgemv_(char *, integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, complex *,
integer *), ccopy_(integer *, complex *, integer *,
complex *, integer *), caxpy_(integer *, complex *, complex *,
integer *, complex *, integer *), ctrmv_(char *, char *, char *,
integer *, complex *, integer *, complex *, integer *), clarfg_(integer *, complex *, complex *, integer
*, complex *), clacgv_(integer *, complex *, integer *);
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLAHRD reduces the first NB columns of a complex general n-by-(n-k+1) */
/* matrix A so that elements below the k-th subdiagonal are zero. The */
/* reduction is performed by a unitary similarity transformation */
/* Q' * A * Q. The routine returns the matrices V and T which determine */
/* Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. */
/* This is an OBSOLETE auxiliary routine. */
/* This routine will be 'deprecated' in a future release. */
/* Please use the new routine CLAHR2 instead. */
/* Arguments */
/* ========= */
/* N (input) INTEGER */
/* The order of the matrix A. */
/* K (input) INTEGER */
/* The offset for the reduction. Elements below the k-th */
/* subdiagonal in the first NB columns are reduced to zero. */
/* NB (input) INTEGER */
/* The number of columns to be reduced. */
/* A (input/output) COMPLEX array, dimension (LDA,N-K+1) */
/* On entry, the n-by-(n-k+1) general matrix A. */
/* On exit, the elements on and above the k-th subdiagonal in */
/* the first NB columns are overwritten with the corresponding */
/* elements of the reduced matrix; the elements below the k-th */
/* subdiagonal, with the array TAU, represent the matrix Q as a */
/* product of elementary reflectors. The other columns of A are */
/* unchanged. See Further Details. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* TAU (output) COMPLEX array, dimension (NB) */
/* The scalar factors of the elementary reflectors. See Further */
/* Details. */
/* T (output) COMPLEX array, dimension (LDT,NB) */
/* The upper triangular matrix T. */
/* LDT (input) INTEGER */
/* The leading dimension of the array T. LDT >= NB. */
/* Y (output) COMPLEX array, dimension (LDY,NB) */
/* The n-by-nb matrix Y. */
/* LDY (input) INTEGER */
/* The leading dimension of the array Y. LDY >= max(1,N). */
/* Further Details */
/* =============== */
/* The matrix Q is represented as a product of nb elementary reflectors */
/* Q = H(1) H(2) . . . H(nb). */
/* 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+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in */
/* A(i+k+1:n,i), and tau in TAU(i). */
//.........这里部分代码省略.........
/* Subroutine */ int claed8_(integer *k, integer *n, integer *qsiz, complex *
q, integer *ldq, real *d__, real *rho, integer *cutpnt, real *z__,
real *dlamda, complex *q2, integer *ldq2, real *w, integer *indxp,
integer *indx, integer *indxq, integer *perm, integer *givptr,
integer *givcol, real *givnum, integer *info)
{
/* System generated locals */
integer q_dim1, q_offset, q2_dim1, q2_offset, i__1;
real r__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
real c__;
integer i__, j;
real s, t;
integer k2, n1, n2, jp, n1p1;
real eps, tau, tol;
integer jlam, imax, jmax;
extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
ccopy_(integer *, complex *, integer *, complex *, integer *),
csrot_(integer *, complex *, integer *, complex *, integer *,
real *, real *), scopy_(integer *, real *, integer *, real *,
integer *);
extern doublereal slapy2_(real *, real *), slamch_(char *);
extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex
*, integer *, complex *, integer *), xerbla_(char *,
integer *);
extern integer isamax_(integer *, real *, integer *);
extern /* Subroutine */ int slamrg_(integer *, integer *, real *, integer
*, integer *, integer *);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLAED8 merges the two sets of eigenvalues together into a single */
/* sorted set. Then it tries to deflate the size of the problem. */
/* There are two ways in which deflation can occur: when two or more */
/* eigenvalues are close together or if there is a tiny element in the */
/* Z vector. For each such occurrence the order of the related secular */
/* equation problem is reduced by one. */
/* Arguments */
/* ========= */
/* K (output) INTEGER */
/* Contains the number of non-deflated eigenvalues. */
/* This is the order of the related secular equation. */
/* N (input) INTEGER */
/* The dimension of the symmetric tridiagonal matrix. N >= 0. */
/* QSIZ (input) INTEGER */
/* The dimension of the unitary matrix used to reduce */
/* the dense or band matrix to tridiagonal form. */
/* QSIZ >= N if ICOMPQ = 1. */
/* Q (input/output) COMPLEX array, dimension (LDQ,N) */
/* On entry, Q contains the eigenvectors of the partially solved */
/* system which has been previously updated in matrix */
/* multiplies with other partially solved eigensystems. */
/* On exit, Q contains the trailing (N-K) updated eigenvectors */
/* (those which were deflated) in its last N-K columns. */
/* LDQ (input) INTEGER */
/* The leading dimension of the array Q. LDQ >= max( 1, N ). */
/* D (input/output) REAL array, dimension (N) */
/* On entry, D contains the eigenvalues of the two submatrices to */
/* be combined. On exit, D contains the trailing (N-K) updated */
/* eigenvalues (those which were deflated) sorted into increasing */
/* order. */
/* RHO (input/output) REAL */
/* Contains the off diagonal element associated with the rank-1 */
/* cut which originally split the two submatrices which are now */
/* being recombined. RHO is modified during the computation to */
/* the value required by SLAED3. */
/* CUTPNT (input) INTEGER */
/* Contains the location of the last eigenvalue in the leading */
/* sub-matrix. MIN(1,N) <= CUTPNT <= N. */
/* Z (input) REAL array, dimension (N) */
/* On input this vector contains the updating vector (the last */
/* row of the first sub-eigenvector matrix and the first row of */
/* the second sub-eigenvector matrix). The contents of Z are */
/* destroyed during the updating process. */
//.........这里部分代码省略.........
/* Subroutine */ int clahqr_(logical *wantt, logical *wantz, integer *n,
integer *ilo, integer *ihi, complex *h__, integer *ldh, complex *w,
integer *iloz, integer *ihiz, complex *z__, integer *ldz, integer *
info)
{
/* System generated locals */
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4;
real r__1, r__2, r__3, r__4, r__5, r__6;
complex q__1, q__2, q__3, q__4, q__5, q__6, q__7;
/* Builtin functions */
double r_imag(complex *);
void r_cnjg(complex *, complex *);
double c_abs(complex *);
void c_sqrt(complex *, complex *), pow_ci(complex *, complex *, integer *)
;
/* Local variables */
integer i__, j, k, l, m;
real s;
complex t, u, v[2], x, y;
integer i1, i2;
complex t1;
real t2;
complex v2;
real aa, ab, ba, bb, h10;
complex h11;
real h21;
complex h22, sc;
integer nh, nz;
real sx;
integer jhi;
complex h11s;
integer jlo, its;
real ulp;
complex sum;
real tst;
complex temp;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *), ccopy_(integer *, complex *, integer *, complex *,
integer *);
real rtemp;
extern /* Subroutine */ int slabad_(real *, real *), clarfg_(integer *,
complex *, complex *, integer *, complex *);
extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);
extern doublereal slamch_(char *);
real safmin, safmax, smlnum;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLAHQR is an auxiliary routine called by CHSEQR to update the */
/* eigenvalues and Schur decomposition already computed by CHSEQR, by */
/* dealing with the Hessenberg submatrix in rows and columns ILO to */
/* IHI. */
/* Arguments */
/* ========= */
/* WANTT (input) LOGICAL */
/* = .TRUE. : the full Schur form T is required; */
/* = .FALSE.: only eigenvalues are required. */
/* WANTZ (input) LOGICAL */
/* = .TRUE. : the matrix of Schur vectors Z is required; */
/* = .FALSE.: Schur vectors are not required. */
/* N (input) INTEGER */
/* The order of the matrix H. N >= 0. */
/* ILO (input) INTEGER */
/* IHI (input) INTEGER */
/* It is assumed that H is already upper triangular in rows and */
/* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). */
/* CLAHQR works primarily with the Hessenberg submatrix in rows */
/* and columns ILO to IHI, but applies transformations to all of */
/* H if WANTT is .TRUE.. */
/* 1 <= ILO <= max(1,IHI); IHI <= N. */
/* H (input/output) COMPLEX array, dimension (LDH,N) */
/* On entry, the upper Hessenberg matrix H. */
/* On exit, if INFO is zero and if WANTT is .TRUE., then H */
/* is upper triangular in rows and columns ILO:IHI. If INFO */
/* is zero and if WANTT is .FALSE., then the contents of H */
/* are unspecified on exit. The output state of H in case */
/* INF is positive is below under the description of INFO. */
/* LDH (input) INTEGER */
/* The leading dimension of the array H. LDH >= max(1,N). */
//.........这里部分代码省略.........
请发表评论