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

intf2c_zdscal(integer* N,           doublereal* alpha,           doublecomplex* X, integer* incX){    zdscal_(N, alpha, X, incX);    return 0;}

/*! double*_zgematrix operator */inline _zgematrix operator*(const double& d, const _zgematrix& mat){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] operator*(const double&, const _zgematrix&)"            << std::endl;#endif//CPPL_VERBOSE    zdscal_(mat.M*mat.N, d, mat.Array, 1);  return mat;}

/*! zcovector*=double operator */inline zcovector& zcovector::operator*=(const double& d){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] zcovector::operator*=(const double&)"            << std::endl;#endif//CPPL_VERBOSE    zdscal_(L, d, Array, 1);  return *this;}

 int zpbtf2_(char *uplo, int *n, int *kd, 	doublecomplex *ab, int *ldab, int *info){    /* System generated locals */    int ab_dim1, ab_offset, i__1, i__2, i__3;    double d__1;    /* Builtin functions */    double sqrt(double);    /* Local variables */    int j, kn;    double ajj;    int kld;    extern  int zher_(char *, int *, double *, 	    doublecomplex *, int *, doublecomplex *, int *);    extern int lsame_(char *, char *);    int upper;    extern  int xerbla_(char *, int *), zdscal_(	    int *, double *, doublecomplex *, int *), zlacgv_(	    int *, doublecomplex *, int *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZPBTF2 computes the Cholesky factorization of a complex Hermitian *//*  positive definite band 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, U' is the conjugate transpose *//*  of U, 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 *//*          Hermitian matrix A is stored: *//*          = 'U':  Upper triangular *//*          = 'L':  Lower triangular *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  KD      (input) INTEGER *//*          The number of super-diagonals of the matrix A if UPLO = 'U', *//*          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0. *//*  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N) *//*          On entry, the upper or lower triangle of the Hermitian band *//*          matrix A, stored in the first KD+1 rows of the array.  The *//*          j-th column of A is stored in the j-th column of the array AB *//*          as follows: *//*          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for MAX(1,j-kd)<=i<=j; *//*          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=MIN(n,j+kd). *//*          On exit, if INFO = 0, the triangular factor U or L from the *//*          Cholesky factorization A = U'*U or A = L*L' of the band *//*          matrix A, in the same storage format as A. *//*  LDAB    (input) INTEGER *//*          The leading dimension of the array AB.  LDAB >= KD+1. *//*  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. *//*  Further Details *//*  =============== *//*  The band storage scheme is illustrated by the following example, when *//*  N = 6, KD = 2, and UPLO = 'U': *//*  On entry:                       On exit: *//*      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46 *//*      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56 *//*     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66 *//*  Similarly, if UPLO = 'L' the format of A is as follows: *//*  On entry:                       On exit: *///.........这里部分代码省略.........

/* Subroutine */ int zlarfg_(integer *n, doublecomplex *alpha, doublecomplex *	x, integer *incx, doublecomplex *tau){/*  -- LAPACK auxiliary 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       =======       ZLARFG generates a complex elementary reflector H of order n, such       that             H' * ( alpha ) = ( beta ),   H' * H = I.                  (   x   )   (   0  )       where alpha and beta are scalars, with beta real, and x is an       (n-1)-element complex vector. H is represented in the form             H = I - tau * ( 1 ) * ( 1 v' ) ,                           ( v )       where tau is a complex scalar and v is a complex (n-1)-element       vector. Note that H is not hermitian.       If the elements of x are all zero and alpha is real, then tau = 0       and H is taken to be the unit matrix.       Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .       Arguments       =========       N       (input) INTEGER               The order of the elementary reflector.       ALPHA   (input/output) COMPLEX*16               On entry, the value alpha.               On exit, it is overwritten with the value beta.       X       (input/output) COMPLEX*16 array, dimension                              (1+(N-2)*abs(INCX))               On entry, the vector x.               On exit, it is overwritten with the vector v.       INCX    (input) INTEGER               The increment between elements of X. INCX > 0.       TAU     (output) COMPLEX*16               The value tau.       =====================================================================          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b5 = {1.,0.};        /* System generated locals */    integer i__1;    doublereal d__1, d__2;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *), d_sign(doublereal *, doublereal *);    /* Local variables */    static doublereal beta;    static integer j;    static doublereal alphi, alphr;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal xnorm;    extern doublereal dlapy3_(doublereal *, doublereal *, doublereal *), 	    dznrm2_(integer *, doublecomplex *, integer *), dlamch_(char *);    static doublereal safmin;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    static doublereal rsafmn;    extern /* Double Complex */ VOID zladiv_(doublecomplex *, doublecomplex *,	     doublecomplex *);    static integer knt;    --x;    /* Function Body */    if (*n <= 0) {	tau->r = 0., tau->i = 0.;	return 0;    }    i__1 = *n - 1;    xnorm = dznrm2_(&i__1, &x[1], incx);    alphr = alpha->r;    alphi = d_imag(alpha);    if (xnorm == 0. && alphi == 0.) {/*        H  =  I *///.........这里部分代码省略.........

/* Subroutine */ int znaitr_(integer *ido, char *bmat, integer *n, integer *k,	 integer *np, integer *nb, doublecomplex *resid, doublereal *rnorm, 	doublecomplex *v, integer *ldv, doublecomplex *h__, integer *ldh, 	integer *ipntr, doublecomplex *workd, integer *info, ftnlen bmat_len){    /* Initialized data */    static logical first = TRUE_;    /* System generated locals */    integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2, i__3;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double d_imag(doublecomplex *), sqrt(doublereal);    /* Local variables */    static integer i__, j;    static real t0, t1, t2, t3, t4, t5;    static integer jj, ipj, irj, ivj;    static doublereal ulp, tst1;    static integer ierr, iter;    static doublereal unfl, ovfl;    static integer itry;    static doublereal temp1;    static logical orth1, orth2, step3, step4;    static doublereal betaj;    static integer infol;    static doublecomplex cnorm;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal rtemp[2];    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    static doublereal wnorm;    extern /* Subroutine */ int dvout_(integer *, integer *, doublereal *, 	    integer *, char *, ftnlen), zcopy_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), ivout_(integer *, integer 	    *, integer *, integer *, char *, ftnlen), zaxpy_(integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zmout_(integer *, integer *, integer *, doublecomplex 	    *, integer *, integer *, char *, ftnlen), zvout_(integer *, 	    integer *, doublecomplex *, integer *, char *, ftnlen);    extern doublereal dlapy2_(doublereal *, doublereal *);    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static doublereal rnorm1;    extern /* Subroutine */ int zgetv0_(integer *, char *, integer *, logical 	    *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublereal *, integer *, doublecomplex *, 	    integer *, ftnlen);    extern doublereal dlamch_(char *, ftnlen);    extern /* Subroutine */ int second_(real *), zdscal_(integer *, 	    doublereal *, doublecomplex *, integer *);    static logical rstart;    static integer msglvl;    static doublereal smlnum;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublecomplex *, ftnlen);    extern /* Subroutine */ int zlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublecomplex *,	     integer *, integer *, ftnlen);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *///.........这里部分代码省略.........

//.........这里部分代码省略........./*              Exit the loop if the growth factor is too small. */		if (grow <= smlnum) {		    goto L90;		}/*              G(j) = ( 1 + CNORM(j) )*G(j-1) */		xj = cnorm[j] + 1.;		grow /= xj;	    }	}L90:	;    }    if (grow * tscal > smlnum) {/*        Use the Level 2 BLAS solve if the reciprocal of the bound on *//*        elements of X is not too small. */	ztrsv_(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1);    } else {/*        Use a Level 1 BLAS solve, scaling intermediate results. */	if (xmax > bignum * .5) {/*           Scale X so that its components are less than or equal to *//*           BIGNUM in absolute value. */	    *scale = bignum * .5 / xmax;	    zdscal_(n, scale, &x[1], &c__1);	    xmax = bignum;	} else {	    xmax *= 2.;	}	if (notran) {/*           Solve A * x = b */	    i__1 = jlast;	    i__2 = jinc;	    for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {/*              Compute x(j) = b(j) / A(j,j), scaling x if necessary. */		i__3 = j;		xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j]), 			abs(d__2));		if (nounit) {		    i__3 = j + j * a_dim1;		    z__1.r = tscal * a[i__3].r, z__1.i = tscal * a[i__3].i;		    tjjs.r = z__1.r, tjjs.i = z__1.i;		} else {		    tjjs.r = tscal, tjjs.i = 0.;		    if (tscal == 1.) {			goto L110;		    }		}		tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(			d__2));		if (tjj > smlnum) {

/* Subroutine */ int zlauu2_(char *uplo, integer *n, doublecomplex *a, 	integer *lda, integer *info){/*  -- LAPACK auxiliary 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       =======       ZLAUU2 computes the product U * U' or L' * L, where the triangular       factor U or L is stored in the upper or lower triangular part of       the array A.       If UPLO = 'U' or 'u' then the upper triangle of the result is stored,       overwriting the factor U in A.       If UPLO = 'L' or 'l' then the lower triangle of the result is stored,       overwriting the factor L in A.       This is the unblocked form of the algorithm, calling Level 2 BLAS.       Arguments       =========       UPLO    (input) CHARACTER*1               Specifies whether the triangular factor stored in the array A               is upper or lower triangular:               = 'U':  Upper triangular               = 'L':  Lower triangular       N       (input) INTEGER               The order of the triangular factor U or L.  N >= 0.       A       (input/output) COMPLEX*16 array, dimension (LDA,N)               On entry, the triangular factor U or L.               On exit, if UPLO = 'U', the upper triangle of A is               overwritten with the upper triangle of the product U * U';               if UPLO = 'L', the lower triangle of A is overwritten with               the lower triangle of the product 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       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer i__;    extern logical lsame_(char *, char *);    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical upper;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *), zlacgv_(	    integer *, doublecomplex *, integer *);    static doublereal aii;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]    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) {	i__1 = -(*info);	xerbla_("ZLAUU2", &i__1);	return 0;    }//.........这里部分代码省略.........

void zdscal( int n, double alpha, doublecomplex *x, int incx){    zdscal_(&n, &alpha, x, &incx);}

/* Subroutine */ int zhpevx_(char *jobz, char *range, char *uplo, integer *n, 	doublecomplex *ap, doublereal *vl, doublereal *vu, integer *il, 	integer *iu, doublereal *abstol, integer *m, doublereal *w, 	doublecomplex *z__, integer *ldz, doublecomplex *work, doublereal *	rwork, integer *iwork, integer *ifail, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2;    doublereal d__1, d__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j, jj;    doublereal eps, vll, vuu, tmp1;    integer indd, inde;    doublereal anrm;    integer imax;    doublereal rmin, rmax;    logical test;    integer itmp1, indee;    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 	    integer *);    doublereal sigma;    extern logical lsame_(char *, char *);    integer iinfo;    char order[1];    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 	    doublereal *, integer *);    logical wantz;    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern doublereal dlamch_(char *);    logical alleig, indeig;    integer iscale, indibl;    logical valeig;    doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    doublereal abstll, bignum;    integer indiwk, indisp, indtau;    extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *, 	     integer *), dstebz_(char *, char *, integer *, doublereal *, 	    doublereal *, integer *, integer *, doublereal *, doublereal *, 	    doublereal *, integer *, integer *, doublereal *, integer *, 	    integer *, doublereal *, integer *, integer *);    extern doublereal zlanhp_(char *, char *, integer *, doublecomplex *, 	    doublereal *);    integer indrwk, indwrk, nsplit;    doublereal smlnum;    extern /* Subroutine */ int zhptrd_(char *, integer *, doublecomplex *, 	    doublereal *, doublereal *, doublecomplex *, integer *), 	    zstein_(integer *, doublereal *, doublereal *, integer *, 	    doublereal *, integer *, integer *, doublecomplex *, integer *, 	    doublereal *, integer *, integer *, integer *), zsteqr_(char *, 	    integer *, doublereal *, doublereal *, doublecomplex *, integer *, 	     doublereal *, integer *), zupgtr_(char *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zupmtr_(char *, char *, char 	    *, integer *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);/*  -- LAPACK driver routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZHPEVX computes selected eigenvalues and, optionally, eigenvectors *//*  of a complex Hermitian matrix A in packed storage. *//*  Eigenvalues/vectors can be selected by specifying either a range of *//*  values or a range of indices for the desired eigenvalues. *//*  Arguments *//*  ========= *//*  JOBZ    (input) CHARACTER*1 *//*          = 'N':  Compute eigenvalues only; *//*          = 'V':  Compute eigenvalues and eigenvectors. *//*  RANGE   (input) CHARACTER*1 *//*          = 'A': all eigenvalues will be found; *//*          = 'V': all eigenvalues in the half-open interval (VL,VU] *//*                 will be found; *//*          = 'I': the IL-th through IU-th eigenvalues will be found. *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  Upper triangle of A is stored; *//*          = 'L':  Lower triangle of A is stored. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *///.........这里部分代码省略.........

//.........这里部分代码省略.........    Information about the permutations P and the diagonal matrix D is       returned in the vector SCALE.       This subroutine is based on the EISPACK routine CBAL.       Modified by Tzu-Yi Chen, Computer Science Division, University of         California at Berkeley, USA       =====================================================================          Test the input parameters          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1, d__2;    /* Builtin functions */    double d_imag(doublecomplex *), z_abs(doublecomplex *);    /* Local variables */    static integer iexc;    static doublereal c__, f, g;    static integer i__, j, k, l, m;    static doublereal r__, s;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    static doublereal sfmin1, sfmin2, sfmax1, sfmax2, ca, ra;    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    static logical noconv;    static integer ica, ira;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --scale;    /* Function Body */    *info = 0;    if (! lsame_(job, "N") && ! lsame_(job, "P") && ! lsame_(job, "S") 	    && ! lsame_(job, "B")) {	*info = -1;    } else if (*n < 0) {	*info = -2;    } else if (*lda < max(1,*n)) {	*info = -4;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZGEBAL", &i__1);	return 0;    }    k = 1;    l = *n;

//.........这里部分代码省略.........        },{6.,7.},{6.,7.},{6.,7.},{6.,7.},{6.,7.},{.1,-.3},{8.,9.},{            .5,            -.1        },{2.,5.},{2.,5.},{2.,5.},{2.,5.},{2.,5.},{.1,.1},{3.,6.},{            -.6,            .1        },{4.,7.},{.1,-.3},{7.,2.},{7.,2.},{7.,2.},{.3,.1},{5.,8.},{            .1,            .4        },{6.,9.},{.4,.1},{8.,3.},{.1,.2},{9.,4.}    };    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),            e_wsle(void);    /* Subroutine */ int s_stop(char *, ftnlen);    /* Local variables */    static integer i__;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *,                                       doublecomplex *, integer *), ctest_(integer *, doublecomplex *,                                               doublecomplex *, doublecomplex *, doublereal *);    static doublecomplex mwpcs[5], mwpct[5];    extern /* Subroutine */ int itest1_(integer *, integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int stest1_(doublereal *, doublereal *,                                        doublereal *, doublereal *);    static doublecomplex cx[8];    extern /* Subroutine */ int zdscal_(integer *, doublereal *,                                        doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    static integer np1, len;    /* Fortran I/O blocks */    static cilist io___19 = { 0, 6, 0, 0, 0 };#define ctrue5_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define ctrue5_ref(a_1,a_2,a_3) ctrue5[ctrue5_subscr(a_1,a_2,a_3)]#define ctrue6_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define ctrue6_ref(a_1,a_2,a_3) ctrue6[ctrue6_subscr(a_1,a_2,a_3)]#define cv_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define cv_ref(a_1,a_2,a_3) cv[cv_subscr(a_1,a_2,a_3)]    for (combla_1.incx = 1; combla_1.incx <= 2; ++combla_1.incx) {        for (np1 = 1; np1 <= 5; ++np1) {            combla_1.n = np1 - 1;            len = max(combla_1.n,1) << 1;            i__1 = len;            for (i__ = 1; i__ <= i__1; ++i__) {                i__2 = i__ - 1;                i__3 = cv_subscr(i__, np1, combla_1.incx);                cx[i__2].r = cv[i__3].r, cx[i__2].i = cv[i__3].i;                /* L20: */            }            if (combla_1.icase == 6) {                d__1 = dznrm2_(&combla_1.n, cx, &combla_1.incx);                stest1_(&d__1, &strue2[np1 - 1], &strue2[np1 - 1], sfac);            } else if (combla_1.icase == 7) {

int zhpev_(char *jobz, char *uplo, int *n, doublecomplex           *ap, double *w, doublecomplex *z__, int *ldz, doublecomplex *           work, double *rwork, int *info){    /* System generated locals */    int z_dim1, z_offset, i__1;    double d__1;    /* Builtin functions */    double sqrt(double);    /* Local variables */    double eps;    int inde;    double anrm;    int imax;    double rmin, rmax;    extern  int dscal_(int *, double *, double *,                       int *);    double sigma;    extern int lsame_(char *, char *);    int iinfo;    int wantz;    extern double dlamch_(char *);    int iscale;    double safmin;    extern  int xerbla_(char *, int *), zdscal_(        int *, double *, doublecomplex *, int *);    double bignum;    int indtau;    extern  int dsterf_(int *, double *, double *,                        int *);    extern double zlanhp_(char *, char *, int *, doublecomplex *,                          double *);    int indrwk, indwrk;    double smlnum;    extern  int zhptrd_(char *, int *, doublecomplex *,                        double *, double *, doublecomplex *, int *),                               zsteqr_(char *, int *, double *, double *,                                       doublecomplex *, int *, double *, int *),                               zupgtr_(char *, int *, doublecomplex *, doublecomplex *,                                       doublecomplex *, int *, doublecomplex *, int *);    /*  -- LAPACK driver routine (version 3.2) -- */    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */    /*     November 2006 */    /*     .. Scalar Arguments .. */    /*     .. */    /*     .. Array Arguments .. */    /*     .. */    /*  Purpose */    /*  ======= */    /*  ZHPEV 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*16 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) DOUBLE PRECISION array, dimension (N) */    /*          If INFO = 0, the eigenvalues in ascending order. */    /*  Z       (output) COMPLEX*16 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). *///.........这里部分代码省略.........

/* Subroutine */ int zlarrv_(integer *n, doublereal *vl, doublereal *vu, 	doublereal *d__, doublereal *l, doublereal *pivmin, integer *isplit, 	integer *m, integer *dol, integer *dou, doublereal *minrgp, 	doublereal *rtol1, doublereal *rtol2, doublereal *w, doublereal *werr, 	 doublereal *wgap, integer *iblock, integer *indexw, doublereal *gers, 	 doublecomplex *z__, integer *ldz, integer *isuppz, doublereal *work, 	integer *iwork, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2;    doublecomplex z__1;    logical L__1;    /* Builtin functions */    double log(doublereal);    /* Local variables */    integer minwsize, i__, j, k, p, q, miniwsize, ii;    doublereal gl;    integer im, in;    doublereal gu, gap, eps, tau, tol, tmp;    integer zto;    doublereal ztz;    integer iend, jblk;    doublereal lgap;    integer done;    doublereal rgap, left;    integer wend, iter;    doublereal bstw;    integer itmp1, indld;    doublereal fudge;    integer idone;    doublereal sigma;    integer iinfo, iindr;    doublereal resid;    logical eskip;    doublereal right;    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 	    doublereal *, integer *);    integer nclus, zfrom;    doublereal rqtol;    integer iindc1, iindc2, indin1, indin2;    logical stp2ii;    extern /* Subroutine */ int zlar1v_(integer *, integer *, integer *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublecomplex *, 	    logical *, integer *, doublereal *, doublereal *, integer *, 	    integer *, doublereal *, doublereal *, doublereal *, doublereal *)	    ;    doublereal lambda;    extern doublereal dlamch_(char *);    integer ibegin, indeig;    logical needbs;    integer indlld;    doublereal sgndef, mingma;    extern /* Subroutine */ int dlarrb_(integer *, doublereal *, doublereal *, 	     integer *, integer *, doublereal *, doublereal *, integer *, 	    doublereal *, doublereal *, doublereal *, doublereal *, integer *, 	     doublereal *, doublereal *, integer *, integer *);    integer oldien, oldncl, wbegin;    doublereal spdiam;    integer negcnt;    extern /* Subroutine */ int dlarrf_(integer *, doublereal *, doublereal *, 	     doublereal *, integer *, integer *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, integer *);    integer oldcls;    doublereal savgap;    integer ndepth;    doublereal ssigma;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    logical usedbs;    integer iindwk, offset;    doublereal gaptol;    integer newcls, oldfst, indwrk, windex, oldlst;    logical usedrq;    integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl;    doublereal bstres;    integer newsiz, zusedu, zusedw;    doublereal nrminv;    logical tryrqc;    integer isupmx;    doublereal rqcorr;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *///.........这里部分代码省略.........

/* Subroutine */int zlatps_(char *uplo, char *trans, char *diag, char * normin, integer *n, doublecomplex *ap, doublecomplex *x, doublereal * scale, doublereal *cnorm, integer *info){    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, ip;    doublereal xj, rec, tjj;    integer jinc, jlen;    doublereal xbnd;    integer imax;    doublereal tmax;    doublecomplex tjjs;    doublereal xmax, grow;    extern /* Subroutine */    int dscal_(integer *, doublereal *, doublereal *, integer *);    extern logical lsame_(char *, char *);    doublereal tscal;    doublecomplex uscal;    integer jlast;    doublecomplex csumj;    extern /* Double Complex */    VOID zdotc_f2c_(doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *);    logical upper;    extern /* Double Complex */    VOID zdotu_f2c_(doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *);    extern /* Subroutine */    int zaxpy_(integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *), ztpsv_( char *, char *, char *, integer *, doublecomplex *, doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    extern doublereal dlamch_(char *);    extern integer idamax_(integer *, doublereal *, integer *);    extern /* Subroutine */    int xerbla_(char *, integer *), zdscal_( integer *, doublereal *, doublecomplex *, integer *);    doublereal bignum;    extern integer izamax_(integer *, doublecomplex *, integer *);    extern /* Double Complex */    VOID zladiv_(doublecomplex *, doublecomplex *, doublecomplex *);    logical notran;    integer jfirst;    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    doublereal 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)    {        *info = -5;    }//.........这里部分代码省略.........

/* Subroutine */ int ztgsja_(char *jobu, char *jobv, char *jobq, integer *m, 	integer *p, integer *n, integer *k, integer *l, doublecomplex *a, 	integer *lda, doublecomplex *b, integer *ldb, doublereal *tola, 	doublereal *tolb, doublereal *alpha, doublereal *beta, doublecomplex *	u, integer *ldu, doublecomplex *v, integer *ldv, doublecomplex *q, 	integer *ldq, doublecomplex *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;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j;    doublereal a1, b1, a3, b3;    doublecomplex a2, b2;    doublereal csq, csu, csv;    doublecomplex snq;    doublereal rwk;    doublecomplex snu, snv;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    doublereal gamma;    extern logical lsame_(char *, char *);    logical initq, initu, initv, wantq, upper;    doublereal error, ssmin;    logical wantu, wantv;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zlags2_(logical *, doublereal *, 	    doublecomplex *, doublereal *, doublereal *, doublecomplex *, 	    doublereal *, doublereal *, doublecomplex *, doublereal *, 	    doublecomplex *, doublereal *, doublecomplex *);    integer kcycle;    extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *), xerbla_(char *, 	    integer *), zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), zlapll_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublereal *), zlaset_(	    char *, integer *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTGSJA computes the generalized singular value decomposition (GSVD) *//*  of two complex 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 ZGGSVP *//*  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 unitary matrices, Z' denotes the conjugate *//*  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 ) *///.........这里部分代码省略.........

/* Subroutine */int zheevx_(char *jobz, char *range, char *uplo, integer *n, doublecomplex *a, integer *lda, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal * w, doublecomplex *z__, integer *ldz, doublecomplex *work, integer * lwork, doublereal *rwork, integer *iwork, integer *ifail, integer * info){    /* System generated locals */    integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;    doublereal d__1, d__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j, nb, jj;    doublereal eps, vll, vuu, tmp1;    integer indd, inde;    doublereal anrm;    integer imax;    doublereal rmin, rmax;    logical test;    integer itmp1, indee;    extern /* Subroutine */    int dscal_(integer *, doublereal *, doublereal *, integer *);    doublereal sigma;    extern logical lsame_(char *, char *);    integer iinfo;    char order[1];    extern /* Subroutine */    int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);    logical lower, wantz;    extern /* Subroutine */    int zswap_(integer *, doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    logical alleig, indeig;    integer iscale, indibl;    logical valeig;    doublereal safmin;    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *);    extern /* Subroutine */    int xerbla_(char *, integer *), zdscal_( integer *, doublereal *, doublecomplex *, integer *);    doublereal abstll, bignum;    extern doublereal zlanhe_(char *, char *, integer *, doublecomplex *, integer *, doublereal *);    integer indiwk, indisp, indtau;    extern /* Subroutine */    int dsterf_(integer *, doublereal *, doublereal *, integer *), dstebz_(char *, char *, integer *, doublereal *, doublereal *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, integer *, doublereal *, integer *, integer *, doublereal *, integer *, integer *);    integer indrwk, indwrk;    extern /* Subroutine */    int zhetrd_(char *, integer *, doublecomplex *, integer *, doublereal *, doublereal *, doublecomplex *, doublecomplex *, integer *, integer *);    integer lwkmin;    extern /* Subroutine */    int zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *);    integer llwork, nsplit;    doublereal smlnum;    extern /* Subroutine */    int zstein_(integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *, doublecomplex *, integer *, doublereal *, integer *, integer *, integer *);    integer lwkopt;    logical lquery;    extern /* Subroutine */    int zsteqr_(char *, integer *, doublereal *, doublereal *, doublecomplex *, integer *, doublereal *, integer *), zungtr_(char *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, integer *), zunmtr_(char *, char *, char *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, integer *);    /* -- 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;    z_dim1 = *ldz;    z_offset = 1 + z_dim1;    z__ -= z_offset;    --work;    --rwork;    --iwork;    --ifail;    /* Function Body */    lower = lsame_(uplo, "L");    wantz = lsame_(jobz, "V");    alleig = lsame_(range, "A");    valeig = lsame_(range, "V");    indeig = lsame_(range, "I");    lquery = *lwork == -1;    *info = 0;    if (! (wantz || lsame_(jobz, "N")))    {        *info = -1;    }//.........这里部分代码省略.........

/*<    >*//* Subroutine */ int ztrevc_(char *side, char *howmny, logical *select,         integer *n, doublecomplex *t, integer *ldt, doublecomplex *vl,         integer *ldvl, doublecomplex *vr, integer *ldvr, integer *mm, integer         *m, doublecomplex *work, doublereal *rwork, integer *info, ftnlen         side_len, ftnlen howmny_len){    /* System generated locals */    integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1,             i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k, ii, ki, is;    doublereal ulp;    logical allv;    doublereal unfl, ovfl, smin;    logical over;    doublereal scale;    extern logical lsame_(const char *, const char *, ftnlen, ftnlen);    doublereal remax;    logical leftv, bothv;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *,             doublecomplex *, doublecomplex *, integer *, doublecomplex *,             integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    logical somev;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *,             doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    extern doublereal dlamch_(char *, ftnlen);    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zdscal_(            integer *, doublereal *, doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    logical rightv;    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    doublereal smlnum;    extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *,             integer *, doublecomplex *, integer *, doublecomplex *,             doublereal *, doublereal *, integer *, ftnlen, ftnlen, ftnlen,             ftnlen);    (void)side_len;    (void)howmny_len;/*  -- 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 .. *//*<       CHARACTER          HOWMNY, SIDE >*//*<       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N >*//*     .. *//*     .. Array Arguments .. *//*<       LOGICAL            SELECT( * ) >*//*<       DOUBLE PRECISION   RWORK( * ) >*//*<    >*//*     .. *//*  Purpose *//*  ======= *//*  ZTREVC computes some or all of the right and/or left eigenvectors of *//*  a complex upper triangular matrix T. *//*  The right eigenvector x and the left eigenvector y of T corresponding *//*  to an eigenvalue w are defined by: *//*               T*x = w*x,     y'*T = w*y' *//*  where y' denotes the conjugate transpose of the vector y. *//*  If all eigenvectors are requested, the routine may either return the *//*  matrices X and/or Y of right or left eigenvectors of T, or the *//*  products Q*X and/or Q*Y, where Q is an input unitary *//*  matrix. If T was obtained from the Schur factorization of an *//*  original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of *//*  right or left eigenvectors of A. *//*  Arguments *//*  ========= *//*  SIDE    (input) CHARACTER*1 *//*          = 'R':  compute right eigenvectors only; *//*          = 'L':  compute left eigenvectors only; *//*          = 'B':  compute both right and left eigenvectors. *//*  HOWMNY  (input) CHARACTER*1 *//*          = 'A':  compute all right and/or left eigenvectors; *//*          = 'B':  compute all right and/or left eigenvectors, *//*                  and backtransform them using the input matrices *//*                  supplied in VR and/or VL; *//*          = 'S':  compute selected right and/or left eigenvectors, *//*                  specified by the logical array SELECT. *//*  SELECT  (input) LOGICAL array, dimension (N) *//*          If HOWMNY = 'S', SELECT specifies the eigenvectors to be *///.........这里部分代码省略.........

/* Subroutine */ int zchkgt_(logical *dotype, integer *nn, integer *nval, 	integer *nns, integer *nsval, doublereal *thresh, logical *tsterr, 	doublecomplex *a, doublecomplex *af, doublecomplex *b, doublecomplex *	x, doublecomplex *xact, doublecomplex *work, doublereal *rwork, 	integer *iwork, integer *nout){    /* Initialized data */    static integer iseedy[4] = { 0,0,0,1 };    static char transs[1*3] = "N" "T" "C";    /* Format strings */    static char fmt_9999[] = "(12x,/002N =/002,i5,/002,/002,10x,/002 type"	    " /002,i2,/002, test(/002,i2,/002) = /002,g12.5)";    static char fmt_9997[] = "(/002 NORM ='/002,a1,/002', N =/002,i5,/002"	    ",/002,10x,/002 type /002,i2,/002, test(/002,i2,/002) = /002,g12."	    "5)";    static char fmt_9998[] = "(/002 TRANS='/002,a1,/002', N =/002,i5,/002, N"	    "RHS=/002,i3,/002, type /002,i2,/002, test(/002,i2,/002) = /002,g"	    "12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, j, k, m, n;    doublecomplex z__[3];    integer in, kl, ku, ix, lda;    doublereal cond;    integer mode, koff, imat, info;    char path[3], dist[1];    integer irhs, nrhs;    char norm[1], type__[1];    integer nrun;    extern /* Subroutine */ int alahd_(integer *, char *);    integer nfail, iseed[4];    extern doublereal dget06_(doublereal *, doublereal *);    doublereal rcond;    integer nimat;    doublereal anorm;    integer itran;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublereal *, doublereal *);    char trans[1];    integer izero, nerrs;    extern /* Subroutine */ int zgtt01_(integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *), zgtt02_(char *, integer *, 	     integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *), zgtt05_(char *, integer *, 	    integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *);    logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zlatb4_(char *, integer *, integer *, 	     integer *, char *, integer *, integer *, doublereal *, integer *, 	     doublereal *, char *), alaerh_(char *, 	    char *, integer *, integer *, char *, integer *, integer *, 	    integer *, integer *, integer *, integer *, integer *, integer *, 	    integer *);    doublereal rcondc, rcondi;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), alasum_(char *, integer *, integer *, 	     integer *, integer *);    doublereal rcondo, ainvnm;    logical trfcon;    extern /* Subroutine */ int zerrge_(char *, integer *);    extern doublereal zlangt_(char *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *);    extern /* Subroutine */ int zlagtm_(char *, integer *, integer *, 	    doublereal *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *, 	    integer *), zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zgtcon_(char *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, integer *), 	    zlatms_(integer *, integer *, char *, integer *, char *, 	    doublereal *, integer *, doublereal *, doublereal *, integer *, 	    integer *, char *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zlarnv_(integer *, integer *, 	    integer *, doublecomplex *);    doublereal result[7];    extern /* Subroutine */ int zgtrfs_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, doublereal *, 	    integer *), zgttrf_(integer *, doublecomplex *, //.........这里部分代码省略.........

/* Subroutine */ int zhetf2_(char *uplo, integer *n, doublecomplex *a, 	integer *lda, integer *ipiv, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal d__;    integer i__, j, k;    doublecomplex t;    doublereal r1, d11;    doublecomplex d12;    doublereal d22;    doublecomplex d21;    integer kk, kp;    doublecomplex wk;    doublereal tt;    doublecomplex wkm1, wkp1;    integer imax, jmax;    extern /* Subroutine */ int zher_(char *, integer *, doublereal *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    doublereal alpha;    extern logical lsame_(char *, char *);    integer kstep;    logical upper;    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern doublereal dlapy2_(doublereal *, doublereal *);    doublereal absakk;    extern logical disnan_(doublereal *);    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    doublereal colmax;    extern integer izamax_(integer *, doublecomplex *, integer *);    doublereal rowmax;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZHETF2 computes the factorization of a complex Hermitian matrix A *//*  using the Bunch-Kaufman diagonal pivoting method: *//*     A = U*D*U'  or  A = L*D*L' *//*  where U (or L) is a product of permutation and unit upper (lower) *//*  triangular matrices, U' is the conjugate transpose of U, and D is *//*  Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. *//*  This is the unblocked version of the algorithm, calling Level 2 BLAS. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the upper or lower triangular part of the *//*          Hermitian matrix A is stored: *//*          = 'U':  Upper triangular *//*          = 'L':  Lower triangular *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading *//*          n-by-n upper triangular part of A contains the upper *//*          triangular part of the matrix A, and the strictly lower *//*          triangular part of A is not referenced.  If UPLO = 'L', the *//*          leading n-by-n lower triangular part of A contains the lower *//*          triangular part of the matrix A, and the strictly upper *//*          triangular part of A is not referenced. *//*          On exit, the block diagonal matrix D and the multipliers used *//*          to obtain the factor U or L (see below for further details). *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,N). *//*  IPIV    (output) INTEGER array, dimension (N) *//*          Details of the interchanges and the block structure of D. *//*          If IPIV(k) > 0, then rows and columns k and IPIV(k) were *//*          interchanged and D(k,k) is a 1-by-1 diagonal block. *//*          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and *//*          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *//*          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) = *///.........这里部分代码省略.........

//.........这里部分代码省略.........    Similarly, if UPLO = 'L' the format of A is as follows:       On entry:                       On exit:          a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66          a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *          a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *       Array elements marked * are not used by the routine.       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static doublereal c_b8 = -1.;    static integer c__1 = 1;        /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3;    doublereal d__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    extern /* Subroutine */ int zher_(char *, integer *, doublereal *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static integer j;    extern logical lsame_(char *, char *);    static logical upper;    static integer kn;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *), zlacgv_(	    integer *, doublecomplex *, integer *);    static doublereal ajj;    static integer kld;#define ab_subscr(a_1,a_2) (a_2)*ab_dim1 + a_1#define ab_ref(a_1,a_2) ab[ab_subscr(a_1,a_2)]    ab_dim1 = *ldab;    ab_offset = 1 + ab_dim1 * 1;    ab -= ab_offset;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    if (! upper && ! lsame_(uplo, "L")) {	*info = -1;    } else if (*n < 0) {	*info = -2;    } else if (*kd < 0) {	*info = -3;    } else if (*ldab < *kd + 1) {	*info = -5;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZPBTF2", &i__1);	return 0;    }/*     Quick return if possible */

/* Subroutine */ int zneigh_(doublereal *rnorm, integer *n, doublecomplex *	h__, integer *ldh, doublecomplex *ritz, doublecomplex *bounds, 	doublecomplex *q, integer *ldq, doublecomplex *workl, doublereal *	rwork, integer *ierr){    /* System generated locals */    integer h_dim1, h_offset, q_dim1, q_offset, i__1;    doublereal d__1;    /* Local variables */    static integer j;    static real t0, t1;    static doublecomplex vl[1];    static doublereal temp;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zmout_(integer *, integer *, integer 	    *, doublecomplex *, integer *, integer *, char *, ftnlen), zvout_(	    integer *, integer *, doublecomplex *, integer *, char *, ftnlen);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int second_(real *);    static logical select[1];    static integer msglvl;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, ftnlen), 	    zlahqr_(logical *, logical *, integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *,	     doublecomplex *, integer *, integer *), ztrevc_(char *, char *, 	    logical *, integer *, doublecomplex *, integer *, doublecomplex *,	     integer *, doublecomplex *, integer *, integer *, integer *, 	    doublecomplex *, doublereal *, integer *, ftnlen, ftnlen), 	    zdscal_(integer *, doublereal *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *, ftnlen);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *//*     %------------------------% *//*     | Local Scalars & Arrays | *//*     %------------------------% *//*     %----------------------% *//*     | External Subroutines | *//*     %----------------------% *//*     %--------------------% *//*     | External Functions | *//*     %--------------------% *//*     %-----------------------% *//*     | Executable Statements | *//*     %-----------------------% *//*     %-------------------------------% *//*     | Initialize timing statistics  | *//*     | & message level for debugging | *//*     %-------------------------------% */    /* Parameter adjustments */    --rwork;    --workl;    --bounds;    --ritz;    h_dim1 = *ldh;//.........这里部分代码省略.........

/* Subroutine */ int zdrvpt_(logical *dotype, integer *nn, integer *nval, 	integer *nrhs, doublereal *thresh, logical *tsterr, doublecomplex *a, 	doublereal *d__, doublecomplex *e, doublecomplex *b, doublecomplex *x, 	 doublecomplex *xact, doublecomplex *work, doublereal *rwork, integer 	*nout){    /* Initialized data */    static integer iseedy[4] = { 0,0,0,1 };    /* Format strings */    static char fmt_9999[] = "(1x,a6,/002, N =/002,i5,/002, type /002,i2,"	    "/002, test /002,i2,/002, ratio = /002,g12.5)";    static char fmt_9998[] = "(1x,a6,/002, FACT='/002,a1,/002', N =/002,i5"	    ",/002, type /002,i2,/002, test /002,i2,/002, ratio = /002,g12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    double z_abs(doublecomplex *);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, j, k, n;    doublereal z__[3];    integer k1, ia, in, kl, ku, ix, nt, lda;    char fact[1];    doublereal cond;    integer mode;    doublereal dmax__;    integer imat, info;    char path[3], dist[1], type__[1];    integer nrun, ifact;    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 	    integer *);    integer nfail, iseed[4];    extern doublereal dget06_(doublereal *, doublereal *);    doublereal rcond;    integer nimat;    doublereal anorm;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublereal *, doublereal *), dcopy_(integer *, doublereal *, integer *, doublereal *, 	    integer *);    integer izero, nerrs;    extern /* Subroutine */ int zptt01_(integer *, doublereal *, 	    doublecomplex *, doublereal *, doublecomplex *, doublecomplex *, 	    doublereal *);    logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zptt02_(char *, integer *, integer *, 	     doublereal *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *), zptt05_(	    integer *, integer *, doublereal *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *), zptsv_(integer *, integer *, doublereal *, 	    doublecomplex *, doublecomplex *, integer *, integer *), zlatb4_(	    char *, integer *, integer *, integer *, char *, integer *, 	    integer *, doublereal *, integer *, doublereal *, char *), aladhd_(integer *, char *), alaerh_(char 	    *, char *, integer *, integer *, char *, integer *, integer *, 	    integer *, integer *, integer *, integer *, integer *, integer *, 	    integer *);    extern integer idamax_(integer *, doublereal *, integer *);    doublereal rcondc;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), alasvm_(char *, integer *, integer *, 	     integer *, integer *), dlarnv_(integer *, integer *, 	    integer *, doublereal *);    doublereal ainvnm;    extern doublereal zlanht_(char *, integer *, doublereal *, doublecomplex *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlaptm_(char *, integer *, integer *, doublereal *, 	    doublereal *, doublecomplex *, doublecomplex *, integer *, 	    doublereal *, doublecomplex *, integer *), zlatms_(	    integer *, integer *, char *, integer *, char *, doublereal *, 	    integer *, doublereal *, doublereal *, integer *, integer *, char 	    *, doublecomplex *, integer *, doublecomplex *, integer *), zlarnv_(integer *, integer *, integer *, 	    doublecomplex *);    doublereal result[6];    extern /* Subroutine */ int zpttrf_(integer *, doublereal *, 	    doublecomplex *, integer *), zerrvx_(char *, integer *), 	    zpttrs_(char *, integer *, integer *, doublereal *, doublecomplex 	    *, doublecomplex *, integer *, integer *), zptsvx_(char *, 	     integer *, integer *, doublereal *, doublecomplex *, doublereal *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *);    /* Fortran I/O blocks */    static cilist io___35 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___38 = { 0, 0, 0, fmt_9998, 0 };//.........这里部分代码省略.........

/* ----------------------------------------------------------------------| *//* Subroutine */ int zgexpv(integer *n, integer *m, doublereal *t, 	doublecomplex *v, doublecomplex *w, doublereal *tol, doublereal *	anorm, doublecomplex *wsp, integer *lwsp, integer *iwsp, integer *	liwsp, S_fp matvec, void *matvecdata, integer *itrace, integer *iflag){    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1;    complex q__1;    doublecomplex z__1;    /* Builtin functions */    /* Subroutine */ int s_stop(char *, ftnlen);    double sqrt(doublereal), d_sign(doublereal *, doublereal *), pow_di(	    doublereal *, integer *), pow_dd(doublereal *, doublereal *), 	    d_lg10(doublereal *);    integer i_dnnt(doublereal *);    double d_int(doublereal *);    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_wsle();    double z_abs(doublecomplex *);    /* Local variables */    static integer ibrkflag;    static doublereal step_min__, step_max__;    static integer i__, j;    static doublereal break_tol__;    static integer k1;    static doublereal p1, p2, p3;    static integer ih, mh, iv, ns, mx;    static doublereal xm;    static integer j1v;    static doublecomplex hij;    static doublereal sgn, eps, hj1j, sqr1, beta, hump;    static integer ifree, lfree;    static doublereal t_old__;    static integer iexph;    static doublereal t_new__;    static integer nexph;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal t_now__;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    static integer nstep;    static doublereal t_out__;    static integer nmult;    static doublereal vnorm;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static integer nscale;    static doublereal rndoff;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), zgpadm_(integer *, integer *, 	    doublereal *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *, integer *, integer *, integer *), znchbv_(	    integer *, doublereal *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *);    static doublereal t_step__, avnorm;    static integer ireject;    static doublereal err_loc__;    static integer nreject, mbrkdwn;    static doublereal tbrkdwn, s_error__, x_error__;    /* Fortran I/O blocks */    static cilist io___40 = { 0, 6, 0, 0, 0 };    static cilist io___48 = { 0, 6, 0, 0, 0 };    static cilist io___49 = { 0, 6, 0, 0, 0 };    static cilist io___50 = { 0, 6, 0, 0, 0 };    static cilist io___51 = { 0, 6, 0, 0, 0 };    static cilist io___52 = { 0, 6, 0, 0, 0 };    static cilist io___53 = { 0, 6, 0, 0, 0 };    static cilist io___54 = { 0, 6, 0, 0, 0 };    static cilist io___55 = { 0, 6, 0, 0, 0 };    static cilist io___56 = { 0, 6, 0, 0, 0 };    static cilist io___57 = { 0, 6, 0, 0, 0 };    static cilist io___58 = { 0, 6, 0, 0, 0 };    static cilist io___59 = { 0, 6, 0, 0, 0 };/* -----Purpose----------------------------------------------------------| *//* ---  ZGEXPV computes w = exp(t*A)*v *//*     for a Zomplex (i.e., complex double precision) matrix A *//*     It does not compute the matrix exponential in isolation but *//*     instead, it computes directly the action of the exponential *//*     operator on the operand vector. This way of doing so allows *//*     for addressing large sparse problems. *//*     The method used is based on Krylov subspace projection *//*     techniques and the matrix under consideration interacts only *//*     via the external routine `matvec' performing the matrix-vector *//*     product (matrix-free method). *//* -----Arguments--------------------------------------------------------| *///.........这里部分代码省略.........

//.........这里部分代码省略.........       Decode and test the input parameters          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b2 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, 	    i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static logical allv;    static doublereal unfl, ovfl, smin;    static logical over;    static integer i__, j, k;    static doublereal scale;    extern logical lsame_(char *, char *);    static doublereal remax;    static logical leftv, bothv;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical somev;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    static integer ii, ki;    extern doublereal dlamch_(char *);    static integer is;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    static logical rightv;    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    static doublereal smlnum;    extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublereal *, doublereal *, integer *);    static doublereal ulp;#define t_subscr(a_1,a_2) (a_2)*t_dim1 + a_1#define t_ref(a_1,a_2) t[t_subscr(a_1,a_2)]#define vl_subscr(a_1,a_2) (a_2)*vl_dim1 + a_1#define vl_ref(a_1,a_2) vl[vl_subscr(a_1,a_2)]#define vr_subscr(a_1,a_2) (a_2)*vr_dim1 + a_1#define vr_ref(a_1,a_2) vr[vr_subscr(a_1,a_2)]    --select;    t_dim1 = *ldt;    t_offset = 1 + t_dim1 * 1;    t -= t_offset;    vl_dim1 = *ldvl;    vl_offset = 1 + vl_dim1 * 1;    vl -= vl_offset;    vr_dim1 = *ldvr;    vr_offset = 1 + vr_dim1 * 1;    vr -= vr_offset;    --work;    --rwork;    /* Function Body */    bothv = lsame_(side, "B");

//.........这里部分代码省略.........       Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *), z_div(doublecomplex *, 	    doublecomplex *, doublecomplex *);    /* Local variables */    static integer imax, jmax, j, k;    static doublereal t, alpha;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    static integer kstep;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static doublereal r1;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    static doublecomplex d11, d21, d22;    static integer jb, jj, kk, jp, kp;    static doublereal absakk;    static integer kw;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    static doublereal colmax;    extern /* Subroutine */ int zlacgv_(integer *, doublecomplex *, integer *)	    ;    extern integer izamax_(integer *, doublecomplex *, integer *);    static doublereal rowmax;    static integer kkw;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define w_subscr(a_1,a_2) (a_2)*w_dim1 + a_1#define w_ref(a_1,a_2) w[w_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --ipiv;    w_dim1 = *ldw;    w_offset = 1 + w_dim1 * 1;    w -= w_offset;    /* Function Body */    *info = 0;/*     Initialize ALPHA for use in choosing pivot block size. */    alpha = (sqrt(17.) + 1.) / 8.;    if (lsame_(uplo, "U")) {/*        Factorize the trailing columns of A using the upper triangle             of A and working backwards, and compute the matrix W = U12*D   

/* Subroutine */int zlarfgp_(integer *n, doublecomplex *alpha, doublecomplex *x, integer *incx, doublecomplex *tau){    /* System generated locals */    integer i__1, i__2;    doublereal d__1, d__2;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *), d_sign(doublereal *, doublereal *), z_abs( doublecomplex *);    /* Local variables */    integer j;    doublecomplex savealpha;    integer knt;    doublereal beta, alphi, alphr;    extern /* Subroutine */    int zscal_(integer *, doublecomplex *, doublecomplex *, integer *);    doublereal xnorm;    extern doublereal dlapy2_(doublereal *, doublereal *), dlapy3_(doublereal *, doublereal *, doublereal *), dznrm2_(integer *, doublecomplex * , integer *), dlamch_(char *);    extern /* Subroutine */    int zdscal_(integer *, doublereal *, doublecomplex *, integer *);    doublereal bignum;    extern /* Double Complex */    VOID zladiv_(doublecomplex *, doublecomplex *, doublecomplex *);    doublereal smlnum;    /* -- 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 .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    --x;    /* Function Body */    if (*n <= 0)    {        tau->r = 0., tau->i = 0.;        return 0;    }    i__1 = *n - 1;    xnorm = dznrm2_(&i__1, &x[1], incx);    alphr = alpha->r;    alphi = d_imag(alpha);    if (xnorm == 0.)    {        /* H = [1-alpha/abs(alpha) 0;        0 I], sign chosen so ALPHA >= 0. */        if (alphi == 0.)        {            if (alphr >= 0.)            {                /* When TAU.eq.ZERO, the vector is special-cased to be */                /* all zeros in the application routines. We do not need */                /* to clear it. */                tau->r = 0., tau->i = 0.;            }            else            {                /* However, the application routines rely on explicit */                /* zero checks when TAU.ne.ZERO, and we must clear X. */                tau->r = 2., tau->i = 0.;                i__1 = *n - 1;                for (j = 1;                        j <= i__1;                        ++j)                {                    i__2 = (j - 1) * *incx + 1;                    x[i__2].r = 0.;                    x[i__2].i = 0.; // , expr subst                }                z__1.r = -alpha->r;                z__1.i = -alpha->i; // , expr subst                alpha->r = z__1.r, alpha->i = z__1.i;            }        }        else        {            /* Only "reflecting" the diagonal entry to be real and non-negative. */            xnorm = dlapy2_(&alphr, &alphi);            d__1 = 1. - alphr / xnorm;            d__2 = -alphi / xnorm;            z__1.r = d__1;            z__1.i = d__2; // , expr subst            tau->r = z__1.r, tau->i = z__1.i;            i__1 = *n - 1;            for (j = 1;                    j <= i__1;//.........这里部分代码省略.........

/* Subroutine */ int zggbal_(char *job, integer *n, doublecomplex *a, integer 	*lda, doublecomplex *b, integer *ldb, integer *ilo, integer *ihi, 	doublereal *lscale, doublereal *rscale, doublereal *work, integer *	info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;    doublereal d__1, d__2, d__3;    /* Builtin functions */    double d_lg10(doublereal *), d_imag(doublecomplex *), z_abs(doublecomplex 	    *), d_sign(doublereal *, doublereal *), pow_di(doublereal *, 	    integer *);    /* Local variables */    integer i__, j, k, l, m;    doublereal t;    integer jc;    doublereal ta, tb, tc;    integer ir;    doublereal ew;    integer it, nr, ip1, jp1, lm1;    doublereal cab, rab, ewc, cor, sum;    integer nrp2, icab, lcab;    doublereal beta, coef;    integer irab, lrab;    doublereal basl, cmax;    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 	    integer *);    doublereal coef2, coef5, gamma, alpha;    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 	    integer *);    extern logical lsame_(char *, char *);    doublereal sfmin, sfmax;    integer iflow;    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 	    integer *, doublereal *, integer *);    integer kount;    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern doublereal dlamch_(char *);    doublereal pgamma;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    integer lsfmin;    extern integer izamax_(integer *, doublecomplex *, integer *);    integer lsfmax;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGBAL balances a pair of general complex matrices (A,B).  This *//*  involves, first, permuting A and B by similarity transformations 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 matrices, and improve the *//*  accuracy of the computed eigenvalues and/or eigenvectors in the *//*  generalized eigenvalue problem A*x = lambda*B*x. *//*  Arguments *//*  ========= *//*  JOB     (input) CHARACTER*1 *//*          Specifies the operations to be performed on A and B: *//*          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 *//*                  and RSCALE(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 matrices A and B.  N >= 0. *//*  A       (input/output) COMPLEX*16 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. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= max(1,N). *//*  B       (input/output) COMPLEX*16 array, dimension (LDB,N) *//*          On entry, the input matrix B. *//*          On exit, B is overwritten by the balanced matrix. *//*          If JOB = 'N', B is not referenced. *//*  LDB     (input) INTEGER *///.........这里部分代码省略.........

/* Subroutine */ int zhptrs_(char *uplo, integer *n, integer *nrhs, 	doublecomplex *ap, integer *ipiv, doublecomplex *b, integer *ldb, 	integer *info){    /* System generated locals */    integer b_dim1, b_offset, i__1, i__2;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(	    doublecomplex *, doublecomplex *);    /* Local variables */    integer j, k;    doublereal s;    doublecomplex ak, bk;    integer kc, kp;    doublecomplex akm1, bkm1, akm1k;    extern logical lsame_(char *, char *);    doublecomplex denom;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    logical upper;    extern /* Subroutine */ int zgeru_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), xerbla_(char *, integer *), zdscal_(integer *, doublereal *, doublecomplex *, 	    integer *), zlacgv_(integer *, doublecomplex *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZHPTRS solves a system of linear equations A*X = B with a complex *//*  Hermitian matrix A stored in packed format using the factorization *//*  A = U*D*U**H or A = L*D*L**H computed by ZHPTRF. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the details of the factorization are stored *//*          as an upper or lower triangular matrix. *//*          = 'U':  Upper triangular, form is A = U*D*U**H; *//*          = 'L':  Lower triangular, form is A = L*D*L**H. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  NRHS    (input) INTEGER *//*          The number of right hand sides, i.e., the number of columns *//*          of the matrix B.  NRHS >= 0. *//*  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          The block diagonal matrix D and the multipliers used to *//*          obtain the factor U or L as computed by ZHPTRF, stored as a *//*          packed triangular matrix. *//*  IPIV    (input) INTEGER array, dimension (N) *//*          Details of the interchanges and the block structure of D *//*          as determined by ZHPTRF. *//*  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS) *//*          On entry, the right hand side matrix B. *//*          On exit, the solution matrix X. *//*  LDB     (input) INTEGER *//*          The leading dimension of the array B.  LDB >= max(1,N). *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0: if INFO = -i, the i-th argument had an illegal value *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    --ap;    --ipiv;//.........这里部分代码省略.........

/* Subroutine */ int zhpev_(char *jobz, char *uplo, integer *n, doublecomplex 	*ap, doublereal *w, doublecomplex *z__, integer *ldz, doublecomplex *	work, doublereal *rwork, integer *info){/*  -- LAPACK driver 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       =======       ZHPEV 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*16 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) DOUBLE PRECISION array, dimension (N)               If INFO = 0, the eigenvalues in ascending order.       Z       (output) COMPLEX*16 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*16 array, dimension (max(1, 2*N-1))       RWORK   (workspace) DOUBLE PRECISION 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 */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer z_dim1, z_offset, i__1;    doublereal d__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    static integer inde;    static doublereal anrm;    static integer imax;    static doublereal rmin, rmax;    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 	    integer *);    static doublereal sigma;    extern logical lsame_(char *, char *);    static integer iinfo;    static logical wantz;    extern doublereal dlamch_(char *);    static integer iscale;    static doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    static doublereal bignum;    static integer indtau;//.........这里部分代码省略.........