自学教程：C++ sscal_函数代码示例

/* Subroutine */ int spotf2_(char *uplo, integer *n, real *a, integer *lda, 	integer *info){/*  -- LAPACK routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          February 29, 1992       Purpose       =======       SPOTF2 computes the Cholesky factorization of a real symmetric       positive definite matrix A.       The factorization has the form          A = U' * U ,  if UPLO = 'U', or          A = L  * L',  if UPLO = 'L',       where U is an upper triangular matrix and L is lower triangular.       This is the unblocked version of the algorithm, calling Level 2 BLAS.       Arguments       =========       UPLO    (input) CHARACTER*1               Specifies whether the upper or lower triangular part of the               symmetric matrix A is stored.               = 'U':  Upper triangular               = 'L':  Lower triangular       N       (input) INTEGER               The order of the matrix A.  N >= 0.       A       (input/output) REAL array, dimension (LDA,N)               On entry, the symmetric matrix A.  If UPLO = 'U', the leading               n by n upper triangular part of A contains the upper               triangular part of the matrix A, and the strictly lower               triangular part of A is not referenced.  If UPLO = 'L', the               leading n by n lower triangular part of A contains the lower               triangular part of the matrix A, and the strictly upper               triangular part of A is not referenced.               On exit, if INFO = 0, the factor U or L from the Cholesky               factorization A = U'*U  or A = L*L'.       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               > 0: if INFO = k, the leading minor of order k is not                    positive definite, and the factorization could not be                    completed.       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;    static real c_b10 = -1.f;    static real c_b12 = 1.f;        /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    real r__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);    static integer j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 	    real *, integer *, real *, real *, integer *);    static logical upper;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real ajj;#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    if (! upper && ! lsame_(uplo, "L")) {	*info = -1;    } else if (*n < 0) {	*info = -2;    } else if (*lda < max(1,*n)) {	*info = -4;    }    if (*info != 0) {//.........这里部分代码省略.........

/* Subroutine */ int sstemr_(char *jobz, char *range, integer *n, real *d__, 	real *e, real *vl, real *vu, integer *il, integer *iu, integer *m, 	real *w, real *z__, integer *ldz, integer *nzc, integer *isuppz, 	logical *tryrac, real *work, integer *lwork, integer *iwork, integer *	liwork, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j;    real r1, r2;    integer jj;    real cs;    integer in;    real sn, wl, wu;    integer iil, iiu;    real eps, tmp;    integer indd, iend, jblk, wend;    real rmin, rmax;    integer itmp;    real tnrm;    integer inde2;    extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)	    ;    integer itmp2;    real rtol1, rtol2, scale;    integer indgp;    extern logical lsame_(char *, char *);    integer iinfo;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    integer iindw, ilast, lwmin;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *), sswap_(integer *, real *, integer *, real *, integer *);    logical wantz;    extern /* Subroutine */ int slaev2_(real *, real *, real *, real *, real *, real *, real *);    logical alleig;    integer ibegin;    logical indeig;    integer iindbl;    logical valeig;    extern doublereal slamch_(char *);    integer wbegin;    real safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    real bignum;    integer inderr, iindwk, indgrs, offset;    extern /* Subroutine */ int slarrc_(char *, integer *, real *, real *, 	    real *, real *, real *, integer *, integer *, integer *, integer *), slarre_(char *, integer *, real *, real *, integer *, 	    integer *, real *, real *, real *, real *, real *, real *, 	    integer *, integer *, integer *, real *, real *, real *, integer *, integer *, real *, real *, real *, integer *, integer *)	    ;    real thresh;    integer iinspl, indwrk, ifirst, liwmin, nzcmin;    real pivmin;    extern doublereal slanst_(char *, integer *, real *, real *);    extern /* Subroutine */ int slarrj_(integer *, real *, real *, integer *, 	    integer *, real *, integer *, real *, real *, real *, integer *, 	    real *, real *, integer *), slarrr_(integer *, real *, real *, 	    integer *);    integer nsplit;    extern /* Subroutine */ int slarrv_(integer *, real *, real *, real *, 	    real *, real *, integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, integer *, integer *, 	    real *, real *, integer *, integer *, real *, integer *, integer *);    real smlnum;    extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);    logical lquery, zquery;/*  -- LAPACK computational routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SSTEMR computes selected eigenvalues and, optionally, eigenvectors *//*  of a real symmetric tridiagonal matrix T. Any such unreduced matrix has *//*  a well defined set of pairwise different real eigenvalues, the corresponding *//*  real eigenvectors are pairwise orthogonal. *//*  The spectrum may be computed either completely or partially by specifying *//*  either an interval (VL,VU] or a range of indices IL:IU for the desired *//*  eigenvalues. *///.........这里部分代码省略.........

/* Subroutine */ int slaror_(char *side, char *init, integer *m, integer *n,        real *a, integer *lda, integer *iseed, real *x, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real r__1;    /* Builtin functions */    double r_sign(real *, real *);    /* Local variables */    static integer kbeg, jcol;    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *,            integer *, real *, integer *, real *, integer *);    static integer irow;    extern real snrm2_(integer *, real *, integer *);    static integer j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),            sgemv_(char *, integer *, integer *, real *, real *, integer *,            real *, integer *, real *, real *, integer *);    static integer ixfrm, itype, nxfrm;    static real xnorm;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real factor;    extern doublereal slarnd_(integer *, integer *);    extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *,            real *, real *, integer *);    static real xnorms;/*  -- LAPACK auxiliary test routine (version 2.0) --       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,       Courant Institute, Argonne National Lab, and Rice University       September 30, 1994    Purpose    =======    SLAROR pre- or post-multiplies an M by N matrix A by a random    orthogonal matrix U, overwriting A.  A may optionally be initialized    to the identity matrix before multiplying by U.  U is generated using    the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).    Arguments    =========    SIDE    (input) CHARACTER*1            Specifies whether A is multiplied on the left or right by U.            = 'L':         Multiply A on the left (premultiply) by U            = 'R':         Multiply A on the right (postmultiply) by U'            = 'C' or 'T':  Multiply A on the left by U and the right                            by U' (Here, U' means U-transpose.)    INIT    (input) CHARACTER*1            Specifies whether or not A should be initialized to the            identity matrix.            = 'I':  Initialize A to (a section of) the identity matrix                     before applying U.            = 'N':  No initialization.  Apply U to the input matrix A.            INIT = 'I' may be used to generate square or rectangular            orthogonal matrices:            For M = N and SIDE = 'L' or 'R', the rows will be orthogonal            to each other, as will the columns.            If M < N, SIDE = 'R' produces a dense matrix whose rows are            orthogonal and whose columns are not, while SIDE = 'L'            produces a matrix whose rows are orthogonal, and whose first            M columns are orthogonal, and whose remaining columns are            zero.            If M > N, SIDE = 'L' produces a dense matrix whose columns            are orthogonal and whose rows are not, while SIDE = 'R'            produces a matrix whose columns are orthogonal, and whose            first M rows are orthogonal, and whose remaining rows are            zero.    M       (input) INTEGER            The number of rows of A.    N       (input) INTEGER            The number of columns of A.    A       (input/output) REAL array, dimension (LDA, N)            On entry, the array A.            On exit, overwritten by U A ( if SIDE = 'L' ),             or by A U ( if SIDE = 'R' ),             or by U A U' ( if SIDE = 'C' or 'T').    LDA     (input) INTEGER            The leading dimension of the array A.  LDA >= max(1,M).//.........这里部分代码省略.........

//.........这里部分代码省略.........			1], &c__1, &tau[i__ - 1]);		e[i__ - 1] = a[i__ - 1 + i__ * a_dim1];		a[i__ - 1 + i__ * a_dim1] = 1.f;/*              Compute W(1:i-1,i) */		i__2 = i__ - 1;		ssymv_("Upper", &i__2, &c_b6, &a[a_offset], lda, &a[i__ * 			a_dim1 + 1], &c__1, &c_b16, &w[iw * w_dim1 + 1], &			c__1);		if (i__ < *n) {		    i__2 = i__ - 1;		    i__3 = *n - i__;		    sgemv_("Transpose", &i__2, &i__3, &c_b6, &w[(iw + 1) * 			    w_dim1 + 1], ldw, &a[i__ * a_dim1 + 1], &c__1, &			    c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);		    i__2 = i__ - 1;		    i__3 = *n - i__;		    sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[(i__ + 1) *			     a_dim1 + 1], lda, &w[i__ + 1 + iw * w_dim1], &			    c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);		    i__2 = i__ - 1;		    i__3 = *n - i__;		    sgemv_("Transpose", &i__2, &i__3, &c_b6, &a[(i__ + 1) * 			    a_dim1 + 1], lda, &a[i__ * a_dim1 + 1], &c__1, &			    c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);		    i__2 = i__ - 1;		    i__3 = *n - i__;		    sgemv_("No transpose", &i__2, &i__3, &c_b5, &w[(iw + 1) * 			    w_dim1 + 1], ldw, &w[i__ + 1 + iw * w_dim1], &			    c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);		}		i__2 = i__ - 1;		sscal_(&i__2, &tau[i__ - 1], &w[iw * w_dim1 + 1], &c__1);		i__2 = i__ - 1;		alpha = tau[i__ - 1] * -.5f * sdot_(&i__2, &w[iw * w_dim1 + 1], &c__1, &a[i__ * a_dim1 + 1], &c__1);		i__2 = i__ - 1;		saxpy_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &w[iw * 			w_dim1 + 1], &c__1);	    }	}    } else {/*        Reduce first NB columns of lower triangle */	i__1 = *nb;	for (i__ = 1; i__ <= i__1; ++i__) {/*           Update A(i:n,i) */	    i__2 = *n - i__ + 1;	    i__3 = i__ - 1;	    sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + a_dim1], lda, 		     &w[i__ + w_dim1], ldw, &c_b6, &a[i__ + i__ * a_dim1], &		    c__1);	    i__2 = *n - i__ + 1;	    i__3 = i__ - 1;	    sgemv_("No transpose", &i__2, &i__3, &c_b5, &w[i__ + w_dim1], ldw, 		     &a[i__ + a_dim1], lda, &c_b6, &a[i__ + i__ * a_dim1], &		    c__1);	    if (i__ < *n) {/*              Generate elementary reflector H(i) to annihilate *//*              A(i+2:n,i) */

 int sgebal_(char *job, int *n, float *a, int *lda, 	int *ilo, int *ihi, float *scale, int *info){    /* System generated locals */    int a_dim1, a_offset, i__1, i__2;    float r__1, r__2;    /* Local variables */    float c__, f, g;    int i__, j, k, l, m;    float r__, s, ca, ra;    int ica, ira, iexc;    extern int lsame_(char *, char *);    extern  int sscal_(int *, float *, float *, int *), 	    sswap_(int *, float *, int *, float *, int *);    float sfmin1, sfmin2, sfmax1, sfmax2;    extern double slamch_(char *);    extern  int xerbla_(char *, int *);    extern int isamax_(int *, float *, int *);    int noconv;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SGEBAL balances a general float matrix A.  This involves, first, *//*  permuting A by a similarity transformation to isolate eigenvalues *//*  in the first 1 to ILO-1 and last IHI+1 to N elements on the *//*  diagonal; and second, applying a diagonal similarity transformation *//*  to rows and columns ILO to IHI to make the rows and columns as *//*  close in norm as possible.  Both steps are optional. *//*  Balancing may reduce the 1-norm of the matrix, and improve the *//*  accuracy of the computed eigenvalues and/or eigenvectors. *//*  Arguments *//*  ========= *//*  JOB     (input) CHARACTER*1 *//*          Specifies the operations to be performed on A: *//*          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0 *//*                  for i = 1,...,N; *//*          = 'P':  permute only; *//*          = 'S':  scale only; *//*          = 'B':  both permute and scale. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  A       (input/output) REAL array, dimension (LDA,N) *//*          On entry, the input matrix A. *//*          On exit,  A is overwritten by the balanced matrix. *//*          If JOB = 'N', A is not referenced. *//*          See Further Details. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= MAX(1,N). *//*  ILO     (output) INTEGER *//*  IHI     (output) INTEGER *//*          ILO and IHI are set to ints such that on exit *//*          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N. *//*          If JOB = 'N' or 'S', ILO = 1 and IHI = N. *//*  SCALE   (output) REAL array, dimension (N) *//*          Details of the permutations and scaling factors applied to *//*          A.  If P(j) is the index of the row and column interchanged *//*          with row and column j and D(j) is the scaling factor *//*          applied to row and column j, then *//*          SCALE(j) = P(j)    for j = 1,...,ILO-1 *//*                   = D(j)    for j = ILO,...,IHI *//*                   = P(j)    for j = IHI+1,...,N. *//*          The order in which the interchanges are made is N to IHI+1, *//*          then 1 to ILO-1. *//*  INFO    (output) INTEGER *//*          = 0:  successful exit. *//*          < 0:  if INFO = -i, the i-th argument had an illegal value. *//*  Further Details *//*  =============== *//*  The permutations consist of row and column interchanges which put *//*  the matrix in the form *//*             ( T1   X   Y  ) *//*     P A P = (  0   B   Z  ) *//*             (  0   0   T2 ) *//*  where T1 and T2 are upper triangular matrices whose eigenvalues lie *//*  along the diagonal.  The column indices ILO and IHI mark the starting *///.........这里部分代码省略.........

/* Subroutine */ int sgbtrf_(integer *m, integer *n, integer *kl, integer *ku, 	 real *ab, integer *ldab, integer *ipiv, integer *info){    /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;    real r__1;    /* Local variables */    integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju, 	    kv, nw;    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 	    integer *, real *, integer *, real *, integer *);    real temp;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemm_(char *, char *, integer *, integer *, integer *, real *, 	    real *, integer *, real *, integer *, real *, real *, integer *);    real work13[4160]	/* was [65][64] */, work31[4160]	/* was [65][	    64] */;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *), sswap_(integer *, real *, integer *, real *, integer *), strsm_(char *, char *, char *, char *, integer *, integer *, 	    real *, real *, integer *, real *, integer *), sgbtf2_(integer *, integer *, integer *, integer 	    *, real *, integer *, integer *, integer *), xerbla_(char *, 	    integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *), isamax_(integer *, real *, 	    integer *);    extern /* Subroutine */ int slaswp_(integer *, real *, integer *, integer 	    *, integer *, integer *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SGBTRF computes an LU factorization of a real m-by-n band matrix A *//*  using partial pivoting with row interchanges. *//*  This is the blocked version of the algorithm, calling Level 3 BLAS. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A.  M >= 0. *//*  N       (input) INTEGER *//*          The number of columns 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. *//*  AB      (input/output) REAL array, dimension (LDAB,N) *//*          On entry, the matrix A in band storage, in rows KL+1 to *//*          2*KL+KU+1; rows 1 to KL of the array need not be set. *//*          The j-th column of A is stored in the j-th column of the *//*          array AB as follows: *//*          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) *//*          On exit, details of the factorization: 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. *//*          See below for further details. *//*  LDAB    (input) INTEGER *//*          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1. *//*  IPIV    (output) INTEGER array, dimension (min(M,N)) *//*          The pivot indices; for 1 <= i <= min(M,N), row i of the *//*          matrix was interchanged with row IPIV(i). *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -i, the i-th argument had an illegal value *//*          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization *//*               has been completed, but the factor U is exactly *//*               singular, and division by zero will occur if it is used *//*               to solve a system of equations. *//*  Further Details *//*  =============== *//*  The band storage scheme is illustrated by the following example, when *//*  M = N = 6, KL = 2, KU = 1: *//*  On entry:                       On exit: *//*      *    *    *    +    +    +       *    *    *   u14  u25  u36 *///.........这里部分代码省略.........

/* Subroutine */ int sgetf2_(integer *m, integer *n, real *a, integer *lda, 	integer *ipiv, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    real r__1;    /* Local variables */    integer i__, j, jp;    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 	    integer *, real *, integer *, real *, integer *), sscal_(integer *, real *, real *, integer *);    real sfmin;    extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *, 	    integer *);    extern doublereal slamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    extern integer isamax_(integer *, real *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SGETF2 computes an LU factorization of a general m-by-n matrix A *//*  using partial pivoting with row interchanges. *//*  The factorization has the form *//*     A = P * L * U *//*  where P is a permutation matrix, L is lower triangular with unit *//*  diagonal elements (lower trapezoidal if m > n), and U is upper *//*  triangular (upper trapezoidal if m < n). *//*  This is the right-looking Level 2 BLAS version of the algorithm. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A.  M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A.  N >= 0. *//*  A       (input/output) REAL array, dimension (LDA,N) *//*          On entry, the m by n matrix to be factored. *//*          On exit, the factors L and U from the factorization *//*          A = P*L*U; the unit diagonal elements of L are not stored. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,M). *//*  IPIV    (output) INTEGER array, dimension (min(M,N)) *//*          The pivot indices; for 1 <= i <= min(M,N), row i of the *//*          matrix was interchanged with row IPIV(i). *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -k, the k-th argument had an illegal value *//*          > 0: if INFO = k, U(k,k) is exactly zero. The factorization *//*               has been completed, but the factor U is exactly *//*               singular, and division by zero will occur if it is used *//*               to solve a system of equations. *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --ipiv;    /* Function Body */    *info = 0;    if (*m < 0) {	*info = -1;    } else if (*n < 0) {	*info = -2;//.........这里部分代码省略.........

/* Subroutine */ int sgeevx_(char *balanc, char *jobvl, char *jobvr, char *	sense, integer *n, real *a, integer *lda, real *wr, real *wi, real *	vl, integer *ldvl, real *vr, integer *ldvr, integer *ilo, integer *	ihi, real *scale, real *abnrm, real *rconde, real *rcondv, real *work, 	 integer *lwork, integer *iwork, integer *info){    /* System generated locals */    integer a_dim1, a_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, 	    i__2, i__3;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, k;    real r__, cs, sn;    char job[1];    real scl, dum[1], eps;    char side[1];    real anrm;    integer ierr, itau, iwrk, nout;    extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 	    integer *, real *, real *);    extern doublereal snrm2_(integer *, real *, integer *);    integer icond;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    extern doublereal slapy2_(real *, real *);    extern /* Subroutine */ int slabad_(real *, real *);    logical scalea;    real cscale;    extern /* Subroutine */ int sgebak_(char *, char *, integer *, integer *, 	    integer *, real *, integer *, real *, integer *, integer *), sgebal_(char *, integer *, real *, integer *, 	    integer *, integer *, real *, integer *);    extern doublereal slamch_(char *), slange_(char *, integer *, 	    integer *, real *, integer *, real *);    extern /* Subroutine */ int sgehrd_(integer *, integer *, integer *, real 	    *, integer *, real *, real *, integer *, integer *), xerbla_(char 	    *, integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *);    logical select[1];    real bignum;    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *);    extern integer isamax_(integer *, real *, integer *);    extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 	    integer *, real *, integer *), slartg_(real *, real *, 	    real *, real *, real *), sorghr_(integer *, integer *, integer *, 	    real *, integer *, real *, real *, integer *, integer *), shseqr_(	    char *, char *, integer *, integer *, integer *, real *, integer *, real *, real *, real *, integer *, real *, integer *, integer *), strevc_(char *, char *, logical *, integer *, 	    real *, integer *, real *, integer *, real *, integer *, integer *, integer *, real *, integer *);    integer minwrk, maxwrk;    extern /* Subroutine */ int strsna_(char *, char *, logical *, integer *, 	    real *, integer *, real *, integer *, real *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *, 	    integer *);    logical wantvl, wntsnb;    integer hswork;    logical wntsne;    real smlnum;    logical lquery, wantvr, wntsnn, wntsnv;/*  -- LAPACK driver routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SGEEVX computes for an N-by-N real nonsymmetric matrix A, the *//*  eigenvalues and, optionally, the left and/or right eigenvectors. *//*  Optionally also, it computes a balancing transformation to improve *//*  the conditioning of the eigenvalues and eigenvectors (ILO, IHI, *//*  SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues *//*  (RCONDE), and reciprocal condition numbers for the right *//*  eigenvectors (RCONDV). *//*  The right eigenvector v(j) of A satisfies *//*                   A * v(j) = lambda(j) * v(j) *//*  where lambda(j) is its eigenvalue. *//*  The left eigenvector u(j) of A satisfies *//*                u(j)**H * A = lambda(j) * u(j)**H *//*  where u(j)**H denotes the conjugate transpose of u(j). *//*  The computed eigenvectors are normalized to have Euclidean norm *//*  equal to 1 and largest component real. *//*  Balancing a matrix means permuting the rows and columns to make it *//*  more nearly upper triangular, and applying a diagonal similarity *///.........这里部分代码省略.........

/* Subroutine */int clatps_(char *uplo, char *trans, char *diag, char * normin, integer *n, complex *ap, complex *x, real *scale, real *cnorm, integer *info){    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    real r__1, r__2, r__3, r__4;    complex q__1, q__2, q__3, q__4;    /* Builtin functions */    double r_imag(complex *);    void r_cnjg(complex *, complex *);    /* Local variables */    integer i__, j, ip;    real xj, rec, tjj;    integer jinc, jlen;    real xbnd;    integer imax;    real tmax;    complex tjjs;    real xmax, grow;    extern /* Complex */    VOID cdotc_f2c_(complex *, integer *, complex *, integer *, complex *, integer *);    extern logical lsame_(char *, char *);    extern /* Subroutine */    int sscal_(integer *, real *, real *, integer *);    real tscal;    complex uscal;    integer jlast;    extern /* Complex */    VOID cdotu_f2c_(complex *, integer *, complex *, integer *, complex *, integer *);    complex csumj;    extern /* Subroutine */    int caxpy_(integer *, complex *, complex *, integer *, complex *, integer *);    logical upper;    extern /* Subroutine */    int ctpsv_(char *, char *, char *, integer *, complex *, complex *, integer *), slabad_( real *, real *);    extern integer icamax_(integer *, complex *, integer *);    extern /* Complex */    VOID cladiv_(complex *, complex *, complex *);    extern real slamch_(char *);    extern /* Subroutine */    int csscal_(integer *, real *, complex *, integer *), xerbla_(char *, integer *);    real bignum;    extern integer isamax_(integer *, real *, integer *);    extern real scasum_(integer *, complex *, integer *);    logical notran;    integer jfirst;    real smlnum;    logical nounit;    /* -- LAPACK auxiliary routine (version 3.4.2) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* September 2012 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Statement Functions .. */    /* .. */    /* .. Statement Function definitions .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    --cnorm;    --x;    --ap;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    notran = lsame_(trans, "N");    nounit = lsame_(diag, "N");    /* Test the input parameters. */    if (! upper && ! lsame_(uplo, "L"))    {        *info = -1;    }    else if (! notran && ! lsame_(trans, "T") && ! lsame_(trans, "C"))    {        *info = -2;    }    else if (! nounit && ! lsame_(diag, "U"))    {        *info = -3;    }    else if (! lsame_(normin, "Y") && ! lsame_(normin, "N"))    {        *info = -4;    }    else if (*n < 0)    {//.........这里部分代码省略.........

int sspevx_(char *jobz, char *range, char *uplo, int *n,            float *ap, float *vl, float *vu, int *il, int *iu, float *abstol,            int *m, float *w, float *z__, int *ldz, float *work, int *            iwork, int *ifail, int *info){    /* System generated locals */    int z_dim1, z_offset, i__1, i__2;    float r__1, r__2;    /* Builtin functions */    double sqrt(double);    /* Local variables */    int i__, j, jj;    float eps, vll, vuu, tmp1;    int indd, inde;    float anrm;    int imax;    float rmin, rmax;    int test;    int itmp1, indee;    float sigma;    extern int lsame_(char *, char *);    int iinfo;    extern  int sscal_(int *, float *, float *, int *);    char order[1];    extern  int scopy_(int *, float *, int *, float *,                       int *), sswap_(int *, float *, int *, float *, int *                                     );    int wantz, alleig, indeig;    int iscale, indibl;    int valeig;    extern double slamch_(char *);    float safmin;    extern  int xerbla_(char *, int *);    float abstll, bignum;    int indtau, indisp, indiwo, indwrk;    extern double slansp_(char *, char *, int *, float *, float *);    extern  int sstein_(int *, float *, float *, int *,                        float *, int *, int *, float *, int *, float *, int *                        , int *, int *), ssterf_(int *, float *, float *,                                int *);    int nsplit;    extern  int sstebz_(char *, char *, int *, float *,                        float *, int *, int *, float *, float *, float *, int *,                        int *, float *, int *, int *, float *, int *,                        int *);    float smlnum;    extern  int sopgtr_(char *, int *, float *, float *,                        float *, int *, float *, int *), ssptrd_(char *,                                int *, float *, float *, float *, float *, int *),                                    ssteqr_(char *, int *, float *, float *, float *, int *,                                            float *, int *), sopmtr_(char *, char *, char *,                                                    int *, int *, float *, float *, float *, int *, float *,                                                    int *);    /*  -- LAPACK driver routine (version 3.2) -- */    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */    /*     November 2006 */    /*     .. Scalar Arguments .. */    /*     .. */    /*     .. Array Arguments .. */    /*     .. */    /*  Purpose */    /*  ======= */    /*  SSPEVX computes selected eigenvalues and, optionally, eigenvectors */    /*  of a float symmetric 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) REAL array, dimension (N*(N+1)/2) */    /*          On entry, the upper or lower triangle of the symmetric 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. *///.........这里部分代码省略.........

/* Subroutine */ int sspgst_(integer *itype, char *uplo, integer *n, real *ap,	 real *bp, integer *info){/*  -- LAPACK routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          March 31, 1993       Purpose       =======       SSPGST reduces a real symmetric-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**T)*A*inv(U) or inv(L)*A*inv(L**T)       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**T or L**T*A*L.       B must have been previously factorized as U**T*U or L*L**T by SPPTRF.       Arguments       =========       ITYPE   (input) INTEGER               = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);               = 2 or 3: compute U*A*U**T or L**T*A*L.       UPLO    (input) CHARACTER               = 'U':  Upper triangle of A is stored and B is factored as                       U**T*U;               = 'L':  Lower triangle of A is stored and B is factored as                       L*L**T.       N       (input) INTEGER               The order of the matrices A and B.  N >= 0.       AP      (input/output) REAL array, dimension (N*(N+1)/2)               On entry, the upper or lower triangle of the symmetric 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) REAL 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 SPPTRF.       INFO    (output) INTEGER               = 0:  successful exit               < 0:  if INFO = -i, the i-th argument had an illegal value       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;    static real c_b9 = -1.f;    static real c_b11 = 1.f;        /* System generated locals */    integer i__1, i__2;    real r__1;    /* Local variables */    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);    extern /* Subroutine */ int sspr2_(char *, integer *, real *, real *, 	    integer *, real *, integer *, real *);    static integer j, k;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static logical upper;    static integer j1, k1;    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 	    real *, integer *), sspmv_(char *, integer *, real *, real *, 	    real *, integer *, real *, real *, integer *), stpmv_(	    char *, char *, char *, integer *, real *, real *, integer *), stpsv_(char *, char *, char *, integer *,	     real *, real *, integer *);    static integer jj, kk;    static real ct;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real ajj;    static integer j1j1;    static real akk;    static integer k1k1;    static real bjj, bkk;    --bp;    --ap;    /* Function Body */    *info = 0;//.........这里部分代码省略.........

/* Subroutine */ int sorgr2_(integer *m, integer *n, integer *k, real *a, 	integer *lda, real *tau, real *work, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    real r__1;    /* Local variables */    integer i__, j, l, ii;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    slarf_(char *, integer *, integer *, real *, integer *, real *, 	    real *, integer *, real *), xerbla_(char *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SORGR2 generates an m by n real matrix Q with orthonormal rows, *//*  which is defined as the last m rows of a product of k elementary *//*  reflectors of order n *//*        Q  =  H(1) H(2) . . . H(k) *//*  as returned by SGERQF. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix Q. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix Q. N >= M. *//*  K       (input) INTEGER *//*          The number of elementary reflectors whose product defines the *//*          matrix Q. M >= K >= 0. *//*  A       (input/output) REAL array, dimension (LDA,N) *//*          On entry, the (m-k+i)-th row must contain the vector which *//*          defines the elementary reflector H(i), for i = 1,2,...,k, as *//*          returned by SGERQF in the last k rows of its array argument *//*          A. *//*          On exit, the m by n matrix Q. *//*  LDA     (input) INTEGER *//*          The first dimension of the array A. LDA >= max(1,M). *//*  TAU     (input) REAL array, dimension (K) *//*          TAU(i) must contain the scalar factor of the elementary *//*          reflector H(i), as returned by SGERQF. *//*  WORK    (workspace) REAL array, dimension (M) *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -i, the i-th argument has an illegal value *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input arguments */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --tau;    --work;    /* Function Body */    *info = 0;    if (*m < 0) {	*info = -1;    } else if (*n < *m) {	*info = -2;    } else if (*k < 0 || *k > *m) {	*info = -3;    } else if (*lda < max(1,*m)) {	*info = -5;    }    if (*info != 0) {//.........这里部分代码省略.........

//.........这里部分代码省略.........    The contents of A on exit are illustrated by the following examples       with nb = 2:       m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):         (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )         (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )         (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )         (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )         (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )         (  v1  v2  a   a   a  )       where a denotes an element of the original matrix which is unchanged,       vi denotes an element of the vector defining H(i), and ui an element       of the vector defining G(i).       =====================================================================          Quick return if possible          Parameter adjustments */    /* Table of constant values */    static real c_b4 = -1.f;    static real c_b5 = 1.f;    static integer c__1 = 1;    static real c_b16 = 0.f;        /* System generated locals */    integer a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset, i__1, i__2, 	    i__3;    /* Local variables */    static integer i__;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 	    real *, integer *, real *, real *, integer *), slarfg_(	    integer *, real *, real *, integer *, real *);#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]#define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1]#define y_ref(a_1,a_2) y[(a_2)*y_dim1 + a_1]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --d__;    --e;    --tauq;    --taup;    x_dim1 = *ldx;    x_offset = 1 + x_dim1 * 1;    x -= x_offset;    y_dim1 = *ldy;    y_offset = 1 + y_dim1 * 1;    y -= y_offset;    /* Function Body */    if (*m <= 0 || *n <= 0) {	return 0;    }    if (*m >= *n) {/*        Reduce to upper bidiagonal form */	i__1 = *nb;

/* Subroutine */ int ssyevr_(char *jobz, char *range, char *uplo, integer *n, 	real *a, integer *lda, real *vl, real *vu, integer *il, integer *iu, 	real *abstol, integer *m, real *w, real *z__, integer *ldz, integer *	isuppz, real *work, integer *lwork, integer *iwork, integer *liwork, 	integer *info){    /* System generated locals */    integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j, nb, jj;    real eps, vll, vuu, tmp1;    integer indd, inde;    real anrm;    integer imax;    real rmin, rmax;    logical test;    integer inddd, indee;    real sigma;    extern logical lsame_(char *, char *);    integer iinfo;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    char order[1];    integer indwk, lwmin;    logical lower;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *), sswap_(integer *, real *, integer *, real *, integer *);    logical wantz, alleig, indeig;    integer iscale, ieeeok, indibl, indifl;    logical valeig;    extern doublereal slamch_(char *);    real safmin;    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *);    extern /* Subroutine */ int xerbla_(char *, integer *);    real abstll, bignum;    integer indtau, indisp, indiwo, indwkn, liwmin;    logical tryrac;    extern /* Subroutine */ int sstein_(integer *, real *, real *, integer *, 	    real *, integer *, integer *, real *, integer *, real *, integer *, integer *, integer *), ssterf_(integer *, real *, real *, 	    integer *);    integer llwrkn, llwork, nsplit;    real smlnum;    extern doublereal slansy_(char *, char *, integer *, real *, integer *, 	    real *);    extern /* Subroutine */ int sstebz_(char *, char *, integer *, real *, 	    real *, integer *, integer *, real *, real *, real *, integer *, 	    integer *, real *, integer *, integer *, real *, integer *, 	    integer *), sstemr_(char *, char *, integer *, 	    real *, real *, real *, real *, integer *, integer *, integer *, 	    real *, real *, integer *, integer *, integer *, logical *, real *, integer *, integer *, integer *, integer *);    integer lwkopt;    logical lquery;    extern /* Subroutine */ int sormtr_(char *, char *, char *, integer *, 	    integer *, real *, integer *, real *, real *, integer *, real *, 	    integer *, integer *), ssytrd_(char *, 	    integer *, real *, integer *, real *, real *, real *, real *, 	    integer *, integer *);/*  -- LAPACK driver routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SSYEVR computes selected eigenvalues and, optionally, eigenvectors *//*  of a real symmetric matrix A.  Eigenvalues and eigenvectors can be *//*  selected by specifying either a range of values or a range of *//*  indices for the desired eigenvalues. *//*  SSYEVR first reduces the matrix A to tridiagonal form T with a call *//*  to SSYTRD.  Then, whenever possible, SSYEVR calls SSTEMR to compute *//*  the eigenspectrum using Relatively Robust Representations.  SSTEMR *//*  computes eigenvalues by the dqds algorithm, while orthogonal *//*  eigenvectors are computed from various "good" L D L^T representations *//*  (also known as Relatively Robust Representations). Gram-Schmidt *//*  orthogonalization is avoided as far as possible. More specifically, *//*  the various steps of the algorithm are as follows. *//*  For each unreduced block (submatrix) of T, *//*     (a) Compute T - sigma I  = L D L^T, so that L and D *//*         define all the wanted eigenvalues to high relative accuracy. *//*         This means that small relative changes in the entries of D and L *//*         cause only small relative changes in the eigenvalues and *//*         eigenvectors. The standard (unfactored) representation of the *//*         tridiagonal matrix T does not have this property in general. *///.........这里部分代码省略.........

/* Subroutine */int chpevd_(char *jobz, char *uplo, integer *n, complex *ap, real *w, complex *z__, integer *ldz, complex *work, integer *lwork, real *rwork, integer *lrwork, integer *iwork, integer *liwork, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1;    real r__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    real eps;    integer inde;    real anrm;    integer imax;    real rmin, rmax, sigma;    extern logical lsame_(char *, char *);    integer iinfo;    extern /* Subroutine */    int sscal_(integer *, real *, real *, integer *);    integer lwmin, llrwk, llwrk;    logical wantz;    integer iscale;    extern real clanhp_(char *, char *, integer *, complex *, real *);    extern /* Subroutine */    int cstedc_(char *, integer *, real *, real *, complex *, integer *, complex *, integer *, real *, integer *, integer *, integer *, integer *);    extern real slamch_(char *);    extern /* Subroutine */    int csscal_(integer *, real *, complex *, integer *);    real safmin;    extern /* Subroutine */    int xerbla_(char *, integer *);    real bignum;    integer indtau;    extern /* Subroutine */    int chptrd_(char *, integer *, complex *, real *, real *, complex *, integer *);    integer indrwk, indwrk, liwmin;    extern /* Subroutine */    int ssterf_(integer *, real *, real *, integer *);    integer lrwmin;    extern /* Subroutine */    int cupmtr_(char *, char *, char *, integer *, integer *, complex *, complex *, complex *, integer *, complex *, integer *);    real smlnum;    logical lquery;    /* -- LAPACK driver routine (version 3.4.0) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* November 2011 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Executable Statements .. */    /* Test the input parameters. */    /* Parameter adjustments */    --ap;    --w;    z_dim1 = *ldz;    z_offset = 1 + z_dim1;    z__ -= z_offset;    --work;    --rwork;    --iwork;    /* Function Body */    wantz = lsame_(jobz, "V");    lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1;    *info = 0;    if (! (wantz || lsame_(jobz, "N")))    {        *info = -1;    }    else if (! (lsame_(uplo, "L") || lsame_(uplo, "U")))    {        *info = -2;    }    else if (*n < 0)    {        *info = -3;    }    else if (*ldz < 1 || wantz && *ldz < *n)    {        *info = -7;    }    if (*info == 0)    {        if (*n <= 1)        {            lwmin = 1;            liwmin = 1;            lrwmin = 1;        }//.........这里部分代码省略.........

//.........这里部分代码省略.........            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.       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer z_dim1, z_offset, i__1, i__2;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    static integer indd, inde;    static real anrm;    static integer imax;    static real rmin, rmax;    static integer itmp1, i__, j, indee;    static real sigma;    extern logical lsame_(char *, char *);    static integer iinfo;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static char order[1];    extern /* Subroutine */ int cswap_(integer *, complex *, integer *, 	    complex *, integer *), scopy_(integer *, real *, integer *, real *	    , integer *);    static logical wantz;    static integer jj;    static logical alleig, indeig;    static integer iscale, indibl;    extern doublereal clanhp_(char *, char *, integer *, complex *, real *);    static logical valeig;    extern doublereal slamch_(char *);    extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 	    *);    static real safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real abstll, bignum;    static integer indiwk, indisp, indtau;    extern /* Subroutine */ int chptrd_(char *, integer *, complex *, real *, 	    real *, complex *, integer *), cstein_(integer *, real *, 	    real *, integer *, real *, integer *, integer *, complex *, 	    integer *, real *, integer *, integer *, integer *);    static integer indrwk, indwrk;    extern /* Subroutine */ int csteqr_(char *, integer *, real *, real *, 	    complex *, integer *, real *, integer *), cupgtr_(char *, 	    integer *, complex *, complex *, complex *, integer *, complex *, 	    integer *), ssterf_(integer *, real *, real *, integer *);    static integer nsplit;    extern /* Subroutine */ int cupmtr_(char *, char *, char *, integer *, 	    integer *, complex *, complex *, complex *, integer *, complex *, 	    integer *);    static real smlnum;    extern /* Subroutine */ int sstebz_(char *, char *, integer *, real *, 

/* Subroutine */ int sggbak_(char *job, char *side, integer *n, integer *ilo, 	integer *ihi, real *lscale, real *rscale, integer *m, real *v, 	integer *ldv, integer *info){/*  -- LAPACK routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          September 30, 1994       Purpose       =======       SGGBAK forms the right or left eigenvectors of a real generalized       eigenvalue problem A*x = lambda*B*x, by backward transformation on       the computed eigenvectors of the balanced pair of matrices output by       SGGBAL.       Arguments       =========       JOB     (input) CHARACTER*1               Specifies the type of backward transformation required:               = 'N':  do nothing, return immediately;               = 'P':  do backward transformation for permutation only;               = 'S':  do backward transformation for scaling only;               = 'B':  do backward transformations for both permutation and                       scaling.               JOB must be the same as the argument JOB supplied to SGGBAL.       SIDE    (input) CHARACTER*1               = 'R':  V contains right eigenvectors;               = 'L':  V contains left eigenvectors.       N       (input) INTEGER               The number of rows of the matrix V.  N >= 0.       ILO     (input) INTEGER       IHI     (input) INTEGER               The integers ILO and IHI determined by SGGBAL.               1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.       LSCALE  (input) REAL array, dimension (N)               Details of the permutations and/or scaling factors applied               to the left side of A and B, as returned by SGGBAL.       RSCALE  (input) REAL array, dimension (N)               Details of the permutations and/or scaling factors applied               to the right side of A and B, as returned by SGGBAL.       M       (input) INTEGER               The number of columns of the matrix V.  M >= 0.       V       (input/output) REAL array, dimension (LDV,M)               On entry, the matrix of right or left eigenvectors to be               transformed, as returned by STGEVC.               On exit, V is overwritten by the transformed eigenvectors.       LDV     (input) INTEGER               The leading dimension of the matrix V. LDV >= max(1,N).       INFO    (output) INTEGER               = 0:  successful exit.               < 0:  if INFO = -i, the i-th argument had an illegal value.       Further Details       ===============       See R.C. Ward, Balancing the generalized eigenvalue problem,                      SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.       =====================================================================          Test the input parameters          Parameter adjustments */    /* System generated locals */    integer v_dim1, v_offset, i__1;    /* Local variables */    static integer i__, k;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static logical leftv;    extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *, 	    integer *), xerbla_(char *, integer *);    static logical rightv;#define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1]    --lscale;    --rscale;    v_dim1 = *ldv;    v_offset = 1 + v_dim1 * 1;    v -= v_offset;    /* Function Body */    rightv = lsame_(side, "R");    leftv = lsame_(side, "L");    *info = 0;//.........这里部分代码省略.........

//.........这里部分代码省略.........            If JOBZ  = 'N' or N <= 1 then LIWORK must be at least 1.               If JOBZ  = 'V' and N > 1 then LIWORK must be at least 2+5*N.       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 E did not converge to zero.       =====================================================================          Test the input parameters.          Parameter adjustments          Function Body */    /* Table of constant values */    static integer c__2 = 2;    static integer c__1 = 1;        /* System generated locals */    integer z_dim1, z_offset, i__1;    real r__1;    /* Builtin functions */    double log(doublereal);    integer pow_ii(integer *, integer *);    double sqrt(doublereal);    /* Local variables */    static real rmin, rmax, tnrm, sigma;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static integer lwmin;    static logical wantz;    static integer iscale;    extern doublereal slamch_(char *);    static real safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real bignum;    extern /* Subroutine */ int sstedc_(char *, integer *, real *, real *, 	    real *, integer *, real *, integer *, integer *, integer *, 	    integer *);    static integer liwmin;    extern doublereal slanst_(char *, integer *, real *, real *);    extern /* Subroutine */ int ssterf_(integer *, real *, real *, integer *);    static real smlnum;    static integer lgn;    static real eps;#define D(I) d[(I)-1]#define E(I) e[(I)-1]#define WORK(I) work[(I)-1]#define IWORK(I) iwork[(I)-1]#define Z(I,J) z[(I)-1 + ((J)-1)* ( *ldz)]    wantz = lsame_(jobz, "V");    *info = 0;    liwmin = 1;    lwmin = 1;    if (! (wantz || lsame_(jobz, "N"))) {

/* Subroutine */ int chpev_(char *jobz, char *uplo, integer *n, complex *ap, 	real *w, complex *z, integer *ldz, complex *work, real *rwork, 	integer *info){/*  -- LAPACK driver routine (version 2.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          March 31, 1993       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.       =====================================================================          Test the input parameters.          Parameter adjustments          Function Body */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer z_dim1, z_offset, i__1;    real r__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    static integer inde;    static real anrm;    static integer imax;    static real rmin, rmax, sigma;    extern logical lsame_(char *, char *);    static integer iinfo;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static logical wantz;    static integer iscale;//.........这里部分代码省略.........

/* Subroutine */int cheev_(char *jobz, char *uplo, integer *n, complex *a, integer *lda, real *w, complex *work, integer *lwork, real *rwork, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real r__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer nb;    real eps;    integer inde;    real anrm;    integer imax;    real rmin, rmax, sigma;    extern logical lsame_(char *, char *);    integer iinfo;    extern /* Subroutine */    int sscal_(integer *, real *, real *, integer *);    logical lower, wantz;    extern real clanhe_(char *, char *, integer *, complex *, integer *, real *);    integer iscale;    extern /* Subroutine */    int clascl_(char *, integer *, integer *, real *, real *, integer *, integer *, complex *, integer *, integer *);    extern real slamch_(char *);    extern /* Subroutine */    int chetrd_(char *, integer *, complex *, integer *, real *, real *, complex *, complex *, integer *, integer *);    real safmin;    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *);    extern /* Subroutine */    int xerbla_(char *, integer *);    real bignum;    integer indtau, indwrk;    extern /* Subroutine */    int csteqr_(char *, integer *, real *, real *, complex *, integer *, real *, integer *), cungtr_(char *, integer *, complex *, integer *, complex *, complex *, integer *, integer *), ssterf_(integer *, real *, real *, integer *);    integer llwork;    real smlnum;    integer lwkopt;    logical lquery;    /* -- LAPACK driver routine (version 3.4.0) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* November 2011 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Executable Statements .. */    /* Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --w;    --work;    --rwork;    /* Function Body */    wantz = lsame_(jobz, "V");    lower = lsame_(uplo, "L");    lquery = *lwork == -1;    *info = 0;    if (! (wantz || lsame_(jobz, "N")))    {        *info = -1;    }    else if (! (lower || lsame_(uplo, "U")))    {        *info = -2;    }    else if (*n < 0)    {        *info = -3;    }    else if (*lda < max(1,*n))    {        *info = -5;    }    if (*info == 0)    {        nb = ilaenv_(&c__1, "CHETRD", uplo, n, &c_n1, &c_n1, &c_n1);        /* Computing MAX */        i__1 = 1;        i__2 = (nb + 1) * *n; // , expr subst        lwkopt = max(i__1,i__2);        work[1].r = (real) lwkopt;        work[1].i = 0.f; // , expr subst        /* Computing MAX */        i__1 = 1;        i__2 = (*n << 1) - 1; // , expr subst//.........这里部分代码省略.........

//.........这里部分代码省略.........    a -= a_offset;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    if (! upper && ! lsame_(uplo, "L")) {	*info = -1;    } else if (*n < 0) {	*info = -2;    } else if (*lda < max(1,*n)) {	*info = -4;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("SPOTF2", &i__1);	return 0;    }/*     Quick return if possible */    if (*n == 0) {	return 0;    }    if (upper) {/*        Compute the Cholesky factorization A = U'*U. */	i__1 = *n;	for (j = 1; j <= i__1; ++j) {/*           Compute U(J,J) and test for non-positive-definiteness. */	    i__2 = j - 1;	    ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j * a_dim1 + 1], &c__1, 		    &a[j * a_dim1 + 1], &c__1);	    if (ajj <= 0.f || sisnan_(&ajj)) {		a[j + j * a_dim1] = ajj;		goto L30;	    }	    ajj = sqrt(ajj);	    a[j + j * a_dim1] = ajj;/*           Compute elements J+1:N of row J. */	    if (j < *n) {		i__2 = j - 1;		i__3 = *n - j;		sgemv_("Transpose", &i__2, &i__3, &c_b10, &a[(j + 1) * a_dim1 			+ 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b12, &a[j + (			j + 1) * a_dim1], lda);		i__2 = *n - j;		r__1 = 1.f / ajj;		sscal_(&i__2, &r__1, &a[j + (j + 1) * a_dim1], lda);	    }	}    } else {/*        Compute the Cholesky factorization A = L*L'. */	i__1 = *n;	for (j = 1; j <= i__1; ++j) {/*           Compute L(J,J) and test for non-positive-definiteness. */	    i__2 = j - 1;	    ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j + a_dim1], lda, &a[j 		    + a_dim1], lda);	    if (ajj <= 0.f || sisnan_(&ajj)) {		a[j + j * a_dim1] = ajj;		goto L30;	    }	    ajj = sqrt(ajj);	    a[j + j * a_dim1] = ajj;/*           Compute elements J+1:N of column J. */	    if (j < *n) {		i__2 = *n - j;		i__3 = j - 1;		sgemv_("No transpose", &i__2, &i__3, &c_b10, &a[j + 1 + 			a_dim1], lda, &a[j + a_dim1], lda, &c_b12, &a[j + 1 + 			j * a_dim1], &c__1);		i__2 = *n - j;		r__1 = 1.f / ajj;		sscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1);	    }	}    }    goto L40;L30:    *info = j;L40:    return 0;/*     End of SPOTF2 */} /* spotf2_ */

 int slaein_(int *rightv, int *noinit, int *n, 	float *h__, int *ldh, float *wr, float *wi, float *vr, float *vi, float 	*b, int *ldb, float *work, float *eps3, float *smlnum, float *bignum, 	int *info){    /* System generated locals */    int b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4;    float r__1, r__2, r__3, r__4;    /* Builtin functions */    double sqrt(double);    /* Local variables */    int i__, j;    float w, x, y;    int i1, i2, i3;    float w1, ei, ej, xi, xr, rec;    int its, ierr;    float temp, norm, vmax;    extern double snrm2_(int *, float *, int *);    float scale;    extern  int sscal_(int *, float *, float *, int *);    char trans[1];    float vcrit;    extern double sasum_(int *, float *, int *);    float rootn, vnorm;    extern double slapy2_(float *, float *);    float absbii, absbjj;    extern int isamax_(int *, float *, int *);    extern  int sladiv_(float *, float *, float *, float *, float *, float *);    char normin[1];    float nrmsml;    extern  int slatrs_(char *, char *, char *, char *, 	    int *, float *, int *, float *, float *, float *, int *);    float growto;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLAEIN uses inverse iteration to find a right or left eigenvector *//*  corresponding to the eigenvalue (WR,WI) of a float upper Hessenberg *//*  matrix H. *//*  Arguments *//*  ========= *//*  RIGHTV   (input) LOGICAL *//*          = .TRUE. : compute right eigenvector; *//*          = .FALSE.: compute left eigenvector. *//*  NOINIT   (input) LOGICAL *//*          = .TRUE. : no initial vector supplied in (VR,VI). *//*          = .FALSE.: initial vector supplied in (VR,VI). *//*  N       (input) INTEGER *//*          The order of the matrix H.  N >= 0. *//*  H       (input) REAL array, dimension (LDH,N) *//*          The upper Hessenberg matrix H. *//*  LDH     (input) INTEGER *//*          The leading dimension of the array H.  LDH >= MAX(1,N). *//*  WR      (input) REAL *//*  WI      (input) REAL *//*          The float and imaginary parts of the eigenvalue of H whose *//*          corresponding right or left eigenvector is to be computed. *//*  VR      (input/output) REAL array, dimension (N) *//*  VI      (input/output) REAL array, dimension (N) *//*          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain *//*          a float starting vector for inverse iteration using the float *//*          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI *//*          must contain the float and imaginary parts of a complex *//*          starting vector for inverse iteration using the complex *//*          eigenvalue (WR,WI); otherwise VR and VI need not be set. *//*          On exit, if WI = 0.0 (float eigenvalue), VR contains the *//*          computed float eigenvector; if WI.ne.0.0 (complex eigenvalue), *//*          VR and VI contain the float and imaginary parts of the *//*          computed complex eigenvector. The eigenvector is normalized *//*          so that the component of largest magnitude has magnitude 1; *//*          here the magnitude of a complex number (x,y) is taken to be *//*          |x| + |y|. *//*          VI is not referenced if WI = 0.0. *//*  B       (workspace) REAL array, dimension (LDB,N) *//*  LDB     (input) INTEGER *//*          The leading dimension of the array B.  LDB >= N+1. *///.........这里部分代码省略.........

/* Subroutine */ int shseqr_(char *job, char *compz, integer *n, integer *ilo,	 integer *ihi, real *h__, integer *ldh, real *wr, real *wi, real *z__,	 integer *ldz, real *work, integer *lwork, integer *info, ftnlen 	job_len, ftnlen compz_len){    /* System generated locals */    address a__1[2];    integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3[2], i__4, 	    i__5;    real r__1, r__2;    char ch__1[2];    /* Builtin functions */    /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);    /* Local variables */    static integer i__, j, k, l;    static real s[225]	/* was [15][15] */, v[16];    static integer i1, i2, ii, nh, nr, ns, nv;    static real vv[16];    static integer itn;    static real tau;    static integer its;    static real ulp, tst1;    static integer maxb;    static real absw;    static integer ierr;    static real unfl, temp, ovfl;    extern logical lsame_(char *, char *, ftnlen, ftnlen);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static integer itemp;    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 	    real *, integer *, real *, integer *, real *, real *, integer *, 	    ftnlen);    static logical initz, wantt;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *);    static logical wantz;    extern doublereal slapy2_(real *, real *);    extern /* Subroutine */ int slabad_(real *, real *);    extern doublereal slamch_(char *, ftnlen);    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *, ftnlen, ftnlen);    extern /* Subroutine */ int slarfg_(integer *, real *, real *, integer *, 	    real *);    extern integer isamax_(integer *, real *, integer *);    extern doublereal slanhs_(char *, integer *, real *, integer *, real *, 	    ftnlen);    extern /* Subroutine */ int slahqr_(logical *, logical *, integer *, 	    integer *, integer *, real *, integer *, real *, real *, integer *	    , integer *, real *, integer *, integer *), slacpy_(char *, 	    integer *, integer *, real *, integer *, real *, integer *, 	    ftnlen), slaset_(char *, integer *, integer *, real *, real *, 	    real *, integer *, ftnlen), slarfx_(char *, integer *, integer *, 	    real *, real *, real *, integer *, real *, ftnlen);    static real smlnum;    static logical lquery;/*  -- LAPACK routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     June 30, 1999 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SHSEQR computes the eigenvalues of a real upper Hessenberg matrix H *//*  and, optionally, the matrices T and Z from the Schur decomposition *//*  H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur *//*  form), and Z is the orthogonal matrix of Schur vectors. *//*  Optionally Z may be postmultiplied into an input orthogonal matrix Q, *//*  so that this routine can give the Schur factorization of a matrix A *//*  which has been reduced to the Hessenberg form H by the orthogonal *//*  matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T. *//*  Arguments *//*  ========= *//*  JOB     (input) CHARACTER*1 *//*          = 'E':  compute eigenvalues only; *//*          = 'S':  compute eigenvalues and the Schur form T. *//*  COMPZ   (input) CHARACTER*1 *//*          = 'N':  no Schur vectors are computed; *//*          = 'I':  Z is initialized to the unit matrix and the matrix Z *//*                  of Schur vectors of H is returned; *//*          = 'V':  Z must contain an orthogonal matrix Q on entry, and *//*                  the product Q*Z is returned. *//*  N       (input) INTEGER *//*          The order of the matrix H.  N >= 0. *///.........这里部分代码省略.........

//.........这里部分代码省略.........    Further Details       ===============       Based on contributions by          Ming Gu and Ren-Cang Li, Computer Science Division, University of            California at Berkeley, USA          Osni Marques, LBNL/NERSC, USA       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static real c_b5 = -1.f;    static integer c__1 = 1;    static real c_b11 = 1.f;    static real c_b13 = 0.f;    static integer c__0 = 0;        /* System generated locals */    integer givcol_dim1, givcol_offset, b_dim1, b_offset, bx_dim1, bx_offset, 	    difr_dim1, difr_offset, givnum_dim1, givnum_offset, poles_dim1, 	    poles_offset, i__1, i__2;    real r__1;    /* Local variables */    static real temp;    extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 	    integer *, real *, real *);    extern doublereal snrm2_(integer *, real *, integer *);    static integer i__, j, m, n;    static real diflj, difrj, dsigj;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 	    real *, integer *, real *, real *, integer *), scopy_(	    integer *, real *, integer *, real *, integer *);    extern doublereal slamc3_(real *, real *);    static real dj;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real dsigjp;    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *), slacpy_(char *, integer *, integer *, real *, integer *, 	    real *, integer *);    static integer nlp1;#define difr_ref(a_1,a_2) difr[(a_2)*difr_dim1 + a_1]#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]#define poles_ref(a_1,a_2) poles[(a_2)*poles_dim1 + a_1]#define bx_ref(a_1,a_2) bx[(a_2)*bx_dim1 + a_1]#define givcol_ref(a_1,a_2) givcol[(a_2)*givcol_dim1 + a_1]#define givnum_ref(a_1,a_2) givnum[(a_2)*givnum_dim1 + a_1]    b_dim1 = *ldb;    b_offset = 1 + b_dim1 * 1;    b -= b_offset;    bx_dim1 = *ldbx;    bx_offset = 1 + bx_dim1 * 1;    bx -= bx_offset;    --perm;    givcol_dim1 = *ldgcol;    givcol_offset = 1 + givcol_dim1 * 1;    givcol -= givcol_offset;    difr_dim1 = *ldgnum;    difr_offset = 1 + difr_dim1 * 1;    difr -= difr_offset;

/* Subroutine */ int sget01_(integer *m, integer *n, real *a, integer *lda, 	real *afac, integer *ldafac, integer *ipiv, real *rwork, real *resid){    /* System generated locals */    integer a_dim1, a_offset, afac_dim1, afac_offset, i__1, i__2;    /* Local variables */    integer i__, j, k;    real t, eps;    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    real anorm;    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 	    real *, integer *, real *, integer *, real *, real *, integer *), strmv_(char *, char *, char *, integer *, real *, 	    integer *, real *, integer *);    extern doublereal slamch_(char *), slange_(char *, integer *, 	    integer *, real *, integer *, real *);    extern /* Subroutine */ int slaswp_(integer *, real *, integer *, integer 	    *, integer *, integer *, integer *);/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SGET01 reconstructs a matrix A from its L*U factorization and *//*  computes the residual *//*     norm(L*U - A) / ( N * norm(A) * EPS ), *//*  where EPS is the machine epsilon. *//*  Arguments *//*  ========== *//*  M       (input) INTEGER *//*          The number of rows of the matrix A.  M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A.  N >= 0. *//*  A       (input) REAL array, dimension (LDA,N) *//*          The original M x N matrix A. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,M). *//*  AFAC    (input/output) REAL array, dimension (LDAFAC,N) *//*          The factored form of the matrix A.  AFAC contains the factors *//*          L and U from the L*U factorization as computed by SGETRF. *//*          Overwritten with the reconstructed matrix, and then with the *//*          difference L*U - A. *//*  LDAFAC  (input) INTEGER *//*          The leading dimension of the array AFAC.  LDAFAC >= max(1,M). *//*  IPIV    (input) INTEGER array, dimension (N) *//*          The pivot indices from SGETRF. *//*  RWORK   (workspace) REAL array, dimension (M) *//*  RESID   (output) REAL *//*          norm(L*U - A) / ( N * norm(A) * EPS ) *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Quick exit if M = 0 or N = 0. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    afac_dim1 = *ldafac;    afac_offset = 1 + afac_dim1;    afac -= afac_offset;    --ipiv;    --rwork;    /* Function Body */    if (*m <= 0 || *n <= 0) {	*resid = 0.f;//.........这里部分代码省略.........

/* Subroutine */ int sorgl2_(integer *m, integer *n, integer *k, real *a, 	integer *lda, real *tau, real *work, integer *info){/*  -- LAPACK routine (version 2.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          February 29, 1992       Purpose       =======       SORGL2 generates an m by n real matrix Q with orthonormal rows,       which is defined as the first m rows of a product of k elementary       reflectors of order n             Q  =  H(k) . . . H(2) H(1)       as returned by SGELQF.       Arguments       =========       M       (input) INTEGER               The number of rows of the matrix Q. M >= 0.       N       (input) INTEGER               The number of columns of the matrix Q. N >= M.       K       (input) INTEGER               The number of elementary reflectors whose product defines the               matrix Q. M >= K >= 0.       A       (input/output) REAL array, dimension (LDA,N)               On entry, the i-th row must contain the vector which defines               the elementary reflector H(i), for i = 1,2,...,k, as returned               by SGELQF in the first k rows of its array argument A.               On exit, the m-by-n matrix Q.       LDA     (input) INTEGER               The first dimension of the array A. LDA >= max(1,M).       TAU     (input) REAL array, dimension (K)               TAU(i) must contain the scalar factor of the elementary               reflector H(i), as returned by SGELQF.       WORK    (workspace) REAL array, dimension (M)       INFO    (output) INTEGER               = 0: successful exit               < 0: if INFO = -i, the i-th argument has an illegal value       =====================================================================          Test the input arguments          Parameter adjustments          Function Body */    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real r__1;    /* Local variables */    static integer i, j, l;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    slarf_(char *, integer *, integer *, real *, integer *, real *, 	    real *, integer *, real *), xerbla_(char *, integer *);#define TAU(I) tau[(I)-1]#define WORK(I) work[(I)-1]#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]    *info = 0;    if (*m < 0) {	*info = -1;    } else if (*n < *m) {	*info = -2;    } else if (*k < 0 || *k > *m) {	*info = -3;    } else if (*lda < max(1,*m)) {	*info = -5;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("SORGL2", &i__1);	return 0;    }/*     Quick return if possible */    if (*m <= 0) {	return 0;    }//.........这里部分代码省略.........

/* Subroutine */ int stgsja_(char *jobu, char *jobv, char *jobq, integer *m, 	integer *p, integer *n, integer *k, integer *l, real *a, integer *lda, 	 real *b, integer *ldb, real *tola, real *tolb, real *alpha, real *	beta, real *u, integer *ldu, real *v, integer *ldv, real *q, integer *	ldq, real *work, integer *ncycle, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, u_dim1, 	    u_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;    real r__1;    /* Local variables */    integer i__, j;    real a1, a2, a3, b1, b2, b3, csq, csu, csv, snq, rwk, snu, snv;    extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 	    integer *, real *, real *);    real gamma;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    logical initq, initu, initv, wantq, upper;    real error, ssmin;    logical wantu, wantv;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *), slags2_(logical *, real *, real *, real *, real *, 	    real *, real *, real *, real *, real *, real *, real *, real *);    integer kcycle;    extern /* Subroutine */ int xerbla_(char *, integer *), slapll_(	    integer *, real *, integer *, real *, integer *, real *), slartg_(	    real *, real *, real *, real *, real *), slaset_(char *, integer *, integer *, real *, real *, real *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  STGSJA computes the generalized singular value decomposition (GSVD) *//*  of two real upper triangular (or trapezoidal) matrices A and B. *//*  On entry, it is assumed that matrices A and B have the following *//*  forms, which may be obtained by the preprocessing subroutine SGGSVP *//*  from a general M-by-N matrix A and P-by-N matrix B: *//*               N-K-L  K    L *//*     A =    K ( 0    A12  A13 ) if M-K-L >= 0; *//*            L ( 0     0   A23 ) *//*        M-K-L ( 0     0    0  ) *//*             N-K-L  K    L *//*     A =  K ( 0    A12  A13 ) if M-K-L < 0; *//*        M-K ( 0     0   A23 ) *//*             N-K-L  K    L *//*     B =  L ( 0     0   B13 ) *//*        P-L ( 0     0    0  ) *//*  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular *//*  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, *//*  otherwise A23 is (M-K)-by-L upper trapezoidal. *//*  On exit, *//*              U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R ), *//*  where U, V and Q are orthogonal matrices, Z' denotes the transpose *//*  of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are *//*  ``diagonal'' matrices, which are of the following structures: *//*  If M-K-L >= 0, *//*                      K  L *//*         D1 =     K ( I  0 ) *//*                  L ( 0  C ) *//*              M-K-L ( 0  0 ) *//*                    K  L *//*         D2 = L   ( 0  S ) *//*              P-L ( 0  0 ) *//*                 N-K-L  K    L *//*    ( 0 R ) = K (  0   R11  R12 ) K *//*              L (  0    0   R22 ) L *//*  where *//*    C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), *//*    S = diag( BETA(K+1),  ... , BETA(K+L) ), *//*    C**2 + S**2 = I. *//*    R is stored in A(1:K+L,N-K-L+1:N) on exit. *//*  If M-K-L < 0, *///.........这里部分代码省略.........

/* Subroutine */ int clatbs_(char *uplo, char *trans, char *diag, char *	normin, integer *n, integer *kd, complex *ab, integer *ldab, complex *	x, real *scale, real *cnorm, integer *info){    /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5;    real r__1, r__2, r__3, r__4;    complex q__1, q__2, q__3, q__4;    /* Builtin functions */    double r_imag(complex *);    void r_cnjg(complex *, complex *);    /* Local variables */    integer i__, j;    real xj, rec, tjj;    integer jinc, jlen;    real xbnd;    integer imax;    real tmax;    complex tjjs;    real xmax, grow;    integer maind;    extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 	    *, complex *, integer *);    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    real tscal;    complex uscal;    integer jlast;    extern /* Complex */ VOID cdotu_(complex *, integer *, complex *, integer 	    *, complex *, integer *);    complex csumj;    extern /* Subroutine */ int ctbsv_(char *, char *, char *, integer *, 	    integer *, complex *, integer *, complex *, integer *), caxpy_(integer *, complex *, complex *, integer *, complex *, integer *);    logical upper;    extern /* Subroutine */ int slabad_(real *, real *);    extern integer icamax_(integer *, complex *, integer *);    extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);    extern doublereal slamch_(char *);    extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 	    *), xerbla_(char *, integer *);    real bignum;    extern integer isamax_(integer *, real *, integer *);    extern doublereal scasum_(integer *, complex *, integer *);    logical notran;    integer jfirst;    real smlnum;    logical nounit;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  CLATBS solves one of the triangular systems *//*     A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b, *//*  with scaling to prevent overflow, where A is an upper or lower *//*  triangular band matrix.  Here A' denotes the transpose of A, x and b *//*  are n-element vectors, and s is a scaling factor, usually less than *//*  or equal to 1, chosen so that the components of x will be less than *//*  the overflow threshold.  If the unscaled problem will not cause *//*  overflow, the Level 2 BLAS routine CTBSV is called.  If the matrix A *//*  is singular (A(j,j) = 0 for some j), then s is set to 0 and a *//*  non-trivial solution to A*x = 0 is returned. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the matrix A is upper or lower triangular. *//*          = 'U':  Upper triangular *//*          = 'L':  Lower triangular *//*  TRANS   (input) CHARACTER*1 *//*          Specifies the operation applied to A. *//*          = 'N':  Solve A * x = s*b     (No transpose) *//*          = 'T':  Solve A**T * x = s*b  (Transpose) *//*          = 'C':  Solve A**H * x = s*b  (Conjugate transpose) *//*  DIAG    (input) CHARACTER*1 *//*          Specifies whether or not the matrix A is unit triangular. *//*          = 'N':  Non-unit triangular *//*          = 'U':  Unit triangular *//*  NORMIN  (input) CHARACTER*1 *//*          Specifies whether CNORM has been set or not. *//*          = 'Y':  CNORM contains the column norms on entry *//*          = 'N':  CNORM is not set on entry.  On exit, the norms will *///.........这里部分代码省略.........

//.........这里部分代码省略.........    } else if (*ldv < max(1,*n)) {	*info = -9;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("SGEBAK", &i__1);	return 0;    }/*     Quick return if possible */    if (*n == 0) {	return 0;    }    if (*m == 0) {	return 0;    }    if (lsame_(job, "N")) {	return 0;    }    if (*ilo == *ihi) {	goto L30;    }/*     Backward balance */    if (lsame_(job, "S") || lsame_(job, "B")) {	if (rightv) {	    i__1 = *ihi;	    for (i__ = *ilo; i__ <= i__1; ++i__) {		s = scale[i__];		sscal_(m, &s, &v[i__ + v_dim1], ldv);	    }	}	if (leftv) {	    i__1 = *ihi;	    for (i__ = *ilo; i__ <= i__1; ++i__) {		s = 1.f / scale[i__];		sscal_(m, &s, &v[i__ + v_dim1], ldv);	    }	}    }/*     Backward permutation *//*     For  I = ILO-1 step -1 until 1, *//*              IHI+1 step 1 until N do -- */L30:    if (lsame_(job, "P") || lsame_(job, "B")) {	if (rightv) {	    i__1 = *n;	    for (ii = 1; ii <= i__1; ++ii) {		i__ = ii;		if (i__ >= *ilo && i__ <= *ihi) {		    goto L40;		}		if (i__ < *ilo) {		    i__ = *ilo - ii;		}		k = scale[i__];		if (k == i__) {