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

/* * Check the inf-norm of the error vector */void pzinf_norm_error(int iam, int_t n, int_t nrhs, doublecomplex x[], int_t ldx,                      doublecomplex xtrue[], int_t ldxtrue, gridinfo_t *grid){    double err, xnorm, temperr, tempxnorm;    doublecomplex *x_work, *xtrue_work;    doublecomplex temp;    int i, j;    for (j = 0; j < nrhs; j++) {        x_work = &x[j*ldx];        xtrue_work = &xtrue[j*ldxtrue];        err = xnorm = 0.0;        for (i = 0; i < n; i++) {            z_sub(&temp, &x_work[i], &xtrue_work[i]);            err = SUPERLU_MAX(err, z_abs(&temp));            xnorm = SUPERLU_MAX(xnorm, z_abs(&x_work[i]));        }        /* get the golbal max err & xnrom */        temperr = err;        tempxnorm = xnorm;        MPI_Allreduce( &temperr, &err, 1, MPI_DOUBLE, MPI_MAX, grid->comm);        MPI_Allreduce( &tempxnorm, &xnorm, 1, MPI_DOUBLE, MPI_MAX, grid->comm);        err = err / xnorm;        if ( !iam ) printf("/tSol %2d: ||X-Xtrue||/||X|| = %e/n", j, err);    }}

/* Subroutine */int zlaev2_(doublecomplex *a, doublecomplex *b, doublecomplex *c__, doublereal *rt1, doublereal *rt2, doublereal *cs1, doublecomplex *sn1){    /* System generated locals */    doublereal d__1, d__2, d__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal t;    doublecomplex w;    extern /* Subroutine */    int dlaev2_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);    /* -- 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 .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Executable Statements .. */    if (z_abs(b) == 0.)    {        w.r = 1.;        w.i = 0.; // , expr subst    }    else    {        d_cnjg(&z__2, b);        d__1 = z_abs(b);        z__1.r = z__2.r / d__1;        z__1.i = z__2.i / d__1; // , expr subst        w.r = z__1.r;        w.i = z__1.i; // , expr subst    }    d__1 = a->r;    d__2 = z_abs(b);    d__3 = c__->r;    dlaev2_(&d__1, &d__2, &d__3, rt1, rt2, cs1, &t);    z__1.r = t * w.r;    z__1.i = t * w.i; // , expr subst    sn1->r = z__1.r, sn1->i = z__1.i;    return 0;    /* End of ZLAEV2 */}

/*! /brief <pre>    Purpose       =======       DZSUM1 takes the sum of the absolute values of a complex       vector and returns a double precision result.       Based on DZASUM from the Level 1 BLAS.       The change is to use the 'genuine' absolute value.       Contributed by Nick Higham for use with ZLACON.       Arguments       =========       N       (input) INT               The number of elements in the vector CX.       CX      (input) COMPLEX*16 array, dimension (N)               The vector whose elements will be summed.       INCX    (input) INT               The spacing between successive values of CX.  INCX > 0.       ===================================================================== </pre>*/  double dzsum1_slu(int *n, doublecomplex *cx, int *incx){    /* Builtin functions */    double z_abs(doublecomplex *);        /* Local variables */    int i, nincx;    double stemp;#define CX(I) cx[(I)-1]    stemp = 0.;    if (*n <= 0) {	return stemp;    }    if (*incx == 1) {	goto L20;    }    /*     CODE FOR INCREMENT NOT EQUAL TO 1 */    nincx = *n * *incx;    for (i = 1; *incx < 0 ? i >= nincx : i <= nincx; i += *incx) {	/*        NEXT LINE MODIFIED. */	stemp += z_abs(&CX(i));/* L10: */    }        return stemp;    /*     CODE FOR INCREMENT EQUAL TO 1 */L20:    for (i = 1; i <= *n; ++i) {	/*        NEXT LINE MODIFIED. */	stemp += z_abs(&CX(i));/* L30: */    }        return stemp;    /*     End of DZSUM1 */} /* dzsum1_slu */

static voidVanVlietAdjustPoles(    const doublecomplex unscaled[4],				/* Unscaled pole locations from Table 1 or				 * 2 of the paper. */    double sigma,		/* Standard deviation */    doublecomplex scaled[4]	/* Scaled pole locations */) {    double q = sigma;				/* Scale factor */    int iter;			/* Iteration number */    double s = VanVlietComputeSigma(q, unscaled);    int i;    /* Search for a scale factor q that yields the desiged sigma. */    for (iter = 0; fabs(sigma-s) > sigma * 1.0e-8; ++iter) {	s = VanVlietComputeSigma(q, unscaled);	q *= sigma/s;    }    /* Adjust poles */    for (i = 0; i < 4; ++i) {	doublecomplex pi = unscaled[i];	double a = pow(z_abs(&pi), 2.0 / q);	double t = atan2(pi.i, pi.r) * 2.0 / q;	scaled[i].r = 1.0/a * cos(t); scaled[i].i = -1.0/a * sin(t);    }}

/*  * Check the inf-norm of the error vector  */void zinf_norm_error_dist(int_t n, int_t nrhs, doublecomplex *x, int_t ldx,			  doublecomplex *xtrue, int_t ldxtrue,                          gridinfo_t *grid){    double err, xnorm;    doublecomplex *x_work, *xtrue_work;    doublecomplex temp;    int i, j;    for (j = 0; j < nrhs; j++) {      x_work = &x[j*ldx];      xtrue_work = &xtrue[j*ldxtrue];      err = xnorm = 0.0;      for (i = 0; i < n; i++) {        z_sub(&temp, &x_work[i], &xtrue_work[i]);	err = SUPERLU_MAX(err, z_abs(&temp));	xnorm = SUPERLU_MAX(xnorm, z_abs(&x_work[i]));      }      err = err / xnorm;      printf("/tRHS %2d: ||X-Xtrue||/||X|| = %e/n", j, err);    }}

/*! /brief Check the inf-norm of the error vector  */void zinf_norm_error(int nrhs, SuperMatrix *X, doublecomplex *xtrue){    DNformat *Xstore;    double err, xnorm;    doublecomplex *Xmat, *soln_work;    doublecomplex temp;    int i, j;    Xstore = X->Store;    Xmat = Xstore->nzval;    for (j = 0; j < nrhs; j++) {      soln_work = &Xmat[j*Xstore->lda];      err = xnorm = 0.0;      for (i = 0; i < X->nrow; i++) {        z_sub(&temp, &soln_work[i], &xtrue[i]);	err = SUPERLU_MAX(err, z_abs(&temp));	xnorm = SUPERLU_MAX(xnorm, z_abs(&soln_work[i]));      }      err = err / xnorm;      printf("||X - Xtrue||/||X|| = %e/n", err);    }}

/* Subroutine */ int PASTEF77(z,rotg)(doublecomplex *ca, doublecomplex *cb, doublereal *c__, doublecomplex *s){    /* System generated locals */    doublereal d__1, d__2;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double z_abs(doublecomplex *);    void bla_z_div(doublecomplex *, doublecomplex *, doublecomplex *);    double sqrt(doublereal);    void bla_d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal norm;    doublecomplex alpha;    doublereal scale;    if (z_abs(ca) != 0.) {	goto L10;    }    *c__ = 0.;    s->real = 1., s->imag = 0.;    ca->real = cb->real, ca->imag = cb->imag;    goto L20;L10:    scale = z_abs(ca) + z_abs(cb);    z__2.real = scale, z__2.imag = 0.;    bla_z_div(&z__1, ca, &z__2);/* Computing 2nd power */    d__1 = z_abs(&z__1);    z__4.real = scale, z__4.imag = 0.;    bla_z_div(&z__3, cb, &z__4);/* Computing 2nd power */    d__2 = z_abs(&z__3);    norm = scale * sqrt(d__1 * d__1 + d__2 * d__2);    d__1 = z_abs(ca);    z__1.real = ca->real / d__1, z__1.imag = ca->imag / d__1;    alpha.real = z__1.real, alpha.imag = z__1.imag;    *c__ = z_abs(ca) / norm;    bla_d_cnjg(&z__3, cb);    z__2.real = alpha.real * z__3.real - alpha.imag * z__3.imag, z__2.imag = alpha.real * z__3.imag + 	    alpha.imag * z__3.real;    z__1.real = z__2.real / norm, z__1.imag = z__2.imag / norm;    s->real = z__1.real, s->imag = z__1.imag;    z__1.real = norm * alpha.real, z__1.imag = norm * alpha.imag;    ca->real = z__1.real, ca->imag = z__1.imag;L20:    return 0;} /* zrotg_ */

static doubleVanVlietComputeSigma(    double sigma,		/* Scale factor */    const doublecomplex poles[4]				/* Poles of the filter */) {    double q = sigma / 2.0;    doublecomplex cs = {0.0, 0.0};    doublecomplex b, c, d, temp;    int i;    for (i = 0; i < 4; ++i) {	doublecomplex pi = poles[i];	double a = pow(z_abs(&pi), -1.0 / q);	double t = atan2(pi.i, pi.r) / q;	b.r = a * cos(t); b.i = a * sin(t);	c.r = 1.0 - b.r; c.i = - b.i;	d.r = c.r * c.r - c.i * c.i; d.i = 2.0 * c.r * c.i;	b.r *= 2.0; b.i *= 2.0;	z_div(&temp, &b, &d);	cs.r += temp.r; cs.i += temp.i;    }    return sqrt(cs.r);}

/* Subroutine */ int ztgex2_(logical *wantq, logical *wantz, integer *n, 	doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, 	doublecomplex *q, integer *ldq, doublecomplex *z__, integer *ldz, 	integer *j1, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 	    z_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    double sqrt(doublereal), z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static doublecomplex f, g;    static integer i__, m;    static doublecomplex s[4]	/* was [2][2] */, t[4]	/* was [2][2] */;    static doublereal cq, sa, sb, cz;    static doublecomplex sq;    static doublereal ss, ws;    static doublecomplex sz;    static doublereal eps, sum;    static logical weak;    static doublecomplex cdum, work[8];    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    static doublereal scale;    extern doublereal dlamch_(char *, ftnlen);    static logical dtrong;    static doublereal thresh;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, ftnlen), 	    zlartg_(doublecomplex *, doublecomplex *, doublereal *, 	    doublecomplex *, doublecomplex *);    static doublereal smlnum;    extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *,	     doublereal *, doublereal *);/*  -- LAPACK auxiliary 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 *//*  ======= *//*  ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22) *//*  in an upper triangular matrix pair (A, B) by an unitary equivalence *//*  transformation. *//*  (A, B) must be in generalized Schur canonical form, that is, A and *//*  B are both upper triangular. *//*  Optionally, the matrices Q and Z of generalized Schur vectors are *//*  updated. *//*         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' *//*         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' *//*  Arguments *//*  ========= *//*  WANTQ   (input) LOGICAL *//*          .TRUE. : update the left transformation matrix Q; *//*          .FALSE.: do not update Q. *//*  WANTZ   (input) LOGICAL *//*          .TRUE. : update the right transformation matrix Z; *//*          .FALSE.: do not update Z. *//*  N       (input) INTEGER *//*          The order of the matrices A and B. N >= 0. *//*  A       (input/output) COMPLEX*16 arrays, dimensions (LDA,N) *//*          On entry, the matrix A in the pair (A, B). *//*          On exit, the updated matrix A. *//*  LDA     (input)  INTEGER *//*          The leading dimension of the array A. LDA >= max(1,N). *//*  B       (input/output) COMPLEX*16 arrays, dimensions (LDB,N) *//*          On entry, the matrix B in the pair (A, B). *//*          On exit, the updated matrix B. *//*  LDB     (input)  INTEGER *//*          The leading dimension of the array B. LDB >= max(1,N). *//*  Q       (input/output) COMPLEX*16 array, dimension (LDZ,N) *//*          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit, *//*          the updated matrix Q. *//*          Not referenced if WANTQ = .FALSE.. *///.........这里部分代码省略.........

/* Subroutine */ int zgetf2_(integer *m, integer *n, doublecomplex *a, 	integer *lda, integer *ipiv, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, jp;    doublereal sfmin;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zgeru_(integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zswap_(integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGETF2 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) COMPLEX*16 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 (itemp < mn) {/*        Initialize partial column norms. The first n elements of *//*        work store the exact column norms. */	i__1 = *n;	for (i__ = itemp + 1; i__ <= i__1; ++i__) {	    i__2 = *m - itemp;	    rwork[i__] = dznrm2_(&i__2, &a[itemp + 1 + i__ * a_dim1], &c__1);	    rwork[*n + i__] = rwork[i__];	}/*        Compute factorization */	i__1 = mn;	for (i__ = itemp + 1; i__ <= i__1; ++i__) {/*           Determine ith pivot column and swap if necessary */	    i__2 = *n - i__ + 1;	    pvt = i__ - 1 + idamax_(&i__2, &rwork[i__], &c__1);	    if (pvt != i__) {		zswap_(m, &a[pvt * a_dim1 + 1], &c__1, &a[i__ * a_dim1 + 1], &			c__1);		itemp = jpvt[pvt];		jpvt[pvt] = jpvt[i__];		jpvt[i__] = itemp;		rwork[pvt] = rwork[i__];		rwork[*n + pvt] = rwork[*n + i__];	    }/*           Generate elementary reflector H(i) */	    i__2 = i__ + i__ * a_dim1;	    aii.r = a[i__2].r, aii.i = a[i__2].i;	    i__2 = *m - i__ + 1;/* Computing MIN */	    i__3 = i__ + 1;	    zlarfp_(&i__2, &aii, &a[min(i__3, *m)+ i__ * a_dim1], &c__1, &tau[		    i__]);	    i__2 = i__ + i__ * a_dim1;	    a[i__2].r = aii.r, a[i__2].i = aii.i;	    if (i__ < *n) {/*              Apply H(i) to A(i:m,i+1:n) from the left */		i__2 = i__ + i__ * a_dim1;		aii.r = a[i__2].r, aii.i = a[i__2].i;		i__2 = i__ + i__ * a_dim1;		a[i__2].r = 1., a[i__2].i = 0.;		i__2 = *m - i__ + 1;		i__3 = *n - i__;		d_cnjg(&z__1, &tau[i__]);		zlarf_("Left", &i__2, &i__3, &a[i__ + i__ * a_dim1], &c__1, &			z__1, &a[i__ + (i__ + 1) * a_dim1], lda, &work[1]);		i__2 = i__ + i__ * a_dim1;		a[i__2].r = aii.r, a[i__2].i = aii.i;	    }/*           Update partial column norms */	    i__2 = *n;	    for (j = i__ + 1; j <= i__2; ++j) {		if (rwork[j] != 0.) {/*                 NOTE: The following 4 lines follow from the analysis in *//*                 Lapack Working Note 176. */		    temp = z_abs(&a[i__ + j * a_dim1]) / rwork[j];/* Computing MAX */		    d__1 = 0., d__2 = (temp + 1.) * (1. - temp);		    temp = max(d__1,d__2);/* Computing 2nd power */		    d__1 = rwork[j] / rwork[*n + j];		    temp2 = temp * (d__1 * d__1);		    if (temp2 <= tol3z) {			if (*m - i__ > 0) {			    i__3 = *m - i__;			    rwork[j] = dznrm2_(&i__3, &a[i__ + 1 + j * a_dim1], &c__1);			    rwork[*n + j] = rwork[j];			} else {			    rwork[j] = 0.;			    rwork[*n + j] = 0.;			}		    } else {			rwork[j] *= sqrt(temp);		    }		}	    }	}    }    return 0;/*     End of ZGEQPF */} /* zgeqpf_ */

/* 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 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 zlaqps_(integer *m, integer *n, integer *offset, integer 	*nb, integer *kb, doublecomplex *a, integer *lda, integer *jpvt, 	doublecomplex *tau, doublereal *vn1, doublereal *vn2, doublecomplex *	auxv, doublecomplex *f, integer *ldf){    /* System generated locals */    integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal);    void d_cnjg(doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    integer i_dnnt(doublereal *);    /* Local variables */    integer j, k, rk;    doublecomplex akk;    integer pvt;    doublereal temp, temp2, tol3z;    integer itemp;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    extern integer idamax_(integer *, doublereal *, integer *);    integer lsticc;    extern /* Subroutine */ int zlarfp_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *);    integer lastrk;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLAQPS computes a step of QR factorization with column pivoting *//*  of a complex M-by-N matrix A by using Blas-3.  It tries to factorize *//*  NB columns from A starting from the row OFFSET+1, and updates all *//*  of the matrix with Blas-3 xGEMM. *//*  In some cases, due to catastrophic cancellations, it cannot *//*  factorize NB columns.  Hence, the actual number of factorized *//*  columns is returned in KB. *//*  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. *//*  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 *//*  OFFSET  (input) INTEGER *//*          The number of rows of A that have been factorized in *//*          previous steps. *//*  NB      (input) INTEGER *//*          The number of columns to factorize. *//*  KB      (output) INTEGER *//*          The number of columns actually factorized. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the M-by-N matrix A. *//*          On exit, block A(OFFSET+1:M,1:KB) is the triangular *//*          factor obtained and block A(1:OFFSET,1:N) has been *//*          accordingly pivoted, but no factorized. *//*          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has *//*          been updated. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= max(1,M). *//*  JPVT    (input/output) INTEGER array, dimension (N) *//*          JPVT(I) = K <==> Column K of the full matrix A has been *//*          permuted into position I in AP. *//*  TAU     (output) COMPLEX*16 array, dimension (KB) *//*          The scalar factors of the elementary reflectors. *//*  VN1     (input/output) DOUBLE PRECISION array, dimension (N) *///.........这里部分代码省略.........

/* Subroutine */ int zhgeqz_(char *job, char *compq, char *compz, integer *n, 	integer *ilo, integer *ihi, doublecomplex *a, integer *lda, 	doublecomplex *b, integer *ldb, doublecomplex *alpha, doublecomplex *	beta, doublecomplex *q, integer *ldq, doublecomplex *z__, integer *	ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *	info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 	    z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2, d__3, d__4, d__5, d__6;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zi(	    doublecomplex *, doublecomplex *, integer *), z_sqrt(	    doublecomplex *, doublecomplex *);    /* Local variables */    static doublereal absb, atol, btol, temp, opst;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    static doublereal temp2, c__;    static integer j;    static doublecomplex s, t;    extern logical lsame_(char *, char *);    static doublecomplex ctemp;    static integer iiter, ilast, jiter;    static doublereal anorm;    static integer maxit;    static doublereal bnorm;    static doublecomplex shift;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal tempr;    static doublecomplex ctemp2, ctemp3;    static logical ilazr2;    static integer jc, in;    static doublereal ascale, bscale;    static doublecomplex u12;    extern doublereal dlamch_(char *);    static integer jr, nq;    static doublecomplex signbc;    static integer nz;    static doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    static doublecomplex eshift;    static logical ilschr;    static integer icompq, ilastm;    static doublecomplex rtdisc;    static integer ischur;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublereal *);    static logical ilazro;    static integer icompz, ifirst;    extern /* Subroutine */ int zlartg_(doublecomplex *, doublecomplex *, 	    doublereal *, doublecomplex *, doublecomplex *);    static integer ifrstm;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static integer istart;    static logical lquery;    static doublecomplex ad11, ad12, ad21, ad22;    static integer jch;    static logical ilq, ilz;    static doublereal ulp;    static doublecomplex abi22;#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 b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]#define q_subscr(a_1,a_2) (a_2)*q_dim1 + a_1#define q_ref(a_1,a_2) q[q_subscr(a_1,a_2)]#define z___subscr(a_1,a_2) (a_2)*z_dim1 + a_1#define z___ref(a_1,a_2) z__[z___subscr(a_1,a_2)]/*  -- LAPACK routine (instrumented to count operations, version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          June 30, 1999          ----------------------- Begin Timing Code ------------------------          Common block to return operation count and iteration count          ITCNT is initialized to 0, OPS is only incremented          OPST is used to accumulate small contributions to OPS          to avoid roundoff error          ------------------------ End Timing Code -------------------------       Purpose       =======       ZHGEQZ implements a single-shift version of the QZ   //.........这里部分代码省略.........

/* ===================================================================== */doublereal zlantb_(char *norm, char *uplo, char *diag, integer *n, integer *k, doublecomplex *ab, integer *ldab, doublereal *work){    /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5;    doublereal ret_val;    /* Builtin functions */    double z_abs(doublecomplex *), sqrt(doublereal);    /* Local variables */    integer i__, j, l;    doublereal sum, scale;    logical udiag;    extern logical lsame_(char *, char *);    doublereal value;    extern logical disnan_(doublereal *);    extern /* Subroutine */    int zlassq_(integer *, doublecomplex *, integer *, doublereal *, doublereal *);    /* -- 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 .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    ab_dim1 = *ldab;    ab_offset = 1 + ab_dim1;    ab -= ab_offset;    --work;    /* Function Body */    if (*n == 0)    {        value = 0.;    }    else if (lsame_(norm, "M"))    {        /* Find max(f2c_abs(A(i,j))). */        if (lsame_(diag, "U"))        {            value = 1.;            if (lsame_(uplo, "U"))            {                i__1 = *n;                for (j = 1;                        j <= i__1;                        ++j)                {                    /* Computing MAX */                    i__2 = *k + 2 - j;                    i__3 = *k;                    for (i__ = max(i__2,1);                            i__ <= i__3;                            ++i__)                    {                        sum = z_abs(&ab[i__ + j * ab_dim1]);                        if (value < sum || disnan_(&sum))                        {                            value = sum;                        }                        /* L10: */                    }                    /* L20: */                }            }            else            {                i__1 = *n;                for (j = 1;                        j <= i__1;                        ++j)                {                    /* Computing MIN */                    i__2 = *n + 1 - j;                    i__4 = *k + 1; // , expr subst                    i__3 = min(i__2,i__4);                    for (i__ = 2;                            i__ <= i__3;                            ++i__)                    {                        sum = z_abs(&ab[i__ + j * ab_dim1]);                        if (value < sum || disnan_(&sum))                        {                            value = sum;                        }                        /* L30: */                    }                    /* L40: *///.........这里部分代码省略.........

//.........这里部分代码省略.........	}L220:/*        Scale by AMAGN */	i__1 = kend;	for (jd = kbeg; jd <= i__1; ++jd) {	    i__2 = jd + jd * a_dim1;	    i__3 = jd + jd * a_dim1;	    d__1 = *amagn * a[i__3].r;	    a[i__2].r = d__1, a[i__2].i = 0.;/* L230: */	}	i__1 = isde;	for (jd = isdb; jd <= i__1; ++jd) {	    i__2 = jd + 1 + jd * a_dim1;	    i__3 = jd + 1 + jd * a_dim1;	    d__1 = *amagn * a[i__3].r;	    a[i__2].r = d__1, a[i__2].i = 0.;/* L240: */	}/*        If RSIGN = .TRUE., assign random signs to diagonal and *//*        subdiagonal */	if (*rsign) {	    i__1 = kend;	    for (jd = kbeg; jd <= i__1; ++jd) {		i__2 = jd + jd * a_dim1;		if (a[i__2].r != 0.) {		    zlarnd_(&z__1, &c__3, &iseed[1]);		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    d__1 = z_abs(&ctemp);		    z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    i__2 = jd + jd * a_dim1;		    i__3 = jd + jd * a_dim1;		    d__1 = a[i__3].r;		    z__1.r = d__1 * ctemp.r, z__1.i = d__1 * ctemp.i;		    a[i__2].r = z__1.r, a[i__2].i = z__1.i;		}/* L250: */	    }	    i__1 = isde;	    for (jd = isdb; jd <= i__1; ++jd) {		i__2 = jd + 1 + jd * a_dim1;		if (a[i__2].r != 0.) {		    zlarnd_(&z__1, &c__3, &iseed[1]);		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    d__1 = z_abs(&ctemp);		    z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    i__2 = jd + 1 + jd * a_dim1;		    i__3 = jd + 1 + jd * a_dim1;		    d__1 = a[i__3].r;		    z__1.r = d__1 * ctemp.r, z__1.i = d__1 * ctemp.i;		    a[i__2].r = z__1.r, a[i__2].i = z__1.i;		}/* L260: */	    }	}/*        Reverse if ITYPE < 0 */	if (*itype < 0) {

/* Subroutine */ int zptcon_(integer *n, doublereal *d__, doublecomplex *e, 	doublereal *anorm, doublereal *rcond, doublereal *rwork, 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       =======       ZPTCON computes the reciprocal of the condition number (in the       1-norm) of a complex Hermitian positive definite tridiagonal matrix       using the factorization A = L*D*L**H or A = U**H*D*U computed by       ZPTTRF.       Norm(inv(A)) is computed by a direct method, and the reciprocal of       the condition number is computed as                        RCOND = 1 / (ANORM * norm(inv(A))).       Arguments       =========       N       (input) INTEGER               The order of the matrix A.  N >= 0.       D       (input) DOUBLE PRECISION array, dimension (N)               The n diagonal elements of the diagonal matrix D from the               factorization of A, as computed by ZPTTRF.       E       (input) COMPLEX*16 array, dimension (N-1)               The (n-1) off-diagonal elements of the unit bidiagonal factor               U or L from the factorization of A, as computed by ZPTTRF.       ANORM   (input) DOUBLE PRECISION               The 1-norm of the original matrix A.       RCOND   (output) DOUBLE PRECISION               The reciprocal of the condition number of the matrix A,               computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the               1-norm of inv(A) computed in this routine.       RWORK   (workspace) DOUBLE PRECISION array, dimension (N)       INFO    (output) INTEGER               = 0:  successful exit               < 0:  if INFO = -i, the i-th argument had an illegal value       Further Details       ===============       The method used is described in Nicholas J. Higham, "Efficient       Algorithms for Computing the Condition Number of a Tridiagonal       Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.       =====================================================================          Test the input arguments.          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer i__1;    doublereal d__1;    /* Builtin functions */    double z_abs(doublecomplex *);    /* Local variables */    static integer i__, ix;    extern integer idamax_(integer *, doublereal *, integer *);    extern /* Subroutine */ int xerbla_(char *, integer *);    static doublereal ainvnm;    --rwork;    --e;    --d__;    /* Function Body */    *info = 0;    if (*n < 0) {	*info = -1;    } else if (*anorm < 0.) {	*info = -4;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZPTCON", &i__1);	return 0;    }/*     Quick return if possible */    *rcond = 0.;    if (*n == 0) {	*rcond = 1.;//.........这里部分代码省略.........

doublereal zqrt14_(char *trans, integer *m, integer *n, integer *nrhs, 	doublecomplex *a, integer *lda, doublecomplex *x, integer *ldx, 	doublecomplex *work, integer *lwork){    /* System generated locals */    integer a_dim1, a_offset, x_dim1, x_offset, i__1, i__2, i__3;    doublereal ret_val, d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static integer info;    static doublereal anrm;    static logical tpsd;    static doublereal xnrm;    static integer i__, j;    extern logical lsame_(char *, char *);    static doublereal rwork[1];    extern /* Subroutine */ int zgelq2_(integer *, integer *, doublecomplex *,	     integer *, doublecomplex *, doublecomplex *, integer *), zgeqr2_(	    integer *, integer *, doublecomplex *, integer *, doublecomplex *,	     doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int zlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublecomplex *,	     integer *, integer *);    static integer ldwork;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal err;#define x_subscr(a_1,a_2) (a_2)*x_dim1 + a_1#define x_ref(a_1,a_2) x[x_subscr(a_1,a_2)]/*  -- LAPACK test 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       =======       ZQRT14 checks whether X is in the row space of A or A'.  It does so       by scaling both X and A such that their norms are in the range       [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X]       (if TRANS = 'C') or an LQ factorization of [A',X]' (if TRANS = 'N'),       and returning the norm of the trailing triangle, scaled by       MAX(M,N,NRHS)*eps.       Arguments       =========       TRANS   (input) CHARACTER*1               = 'N':  No transpose, check for X in the row space of A               = 'C':  Conjugate transpose, check for X in row space of A'.       M       (input) INTEGER               The number of rows of the matrix A.       N       (input) INTEGER               The number of columns of the matrix A.       NRHS    (input) INTEGER               The number of right hand sides, i.e., the number of columns               of X.       A       (input) COMPLEX*16 array, dimension (LDA,N)               The M-by-N matrix A.       LDA     (input) INTEGER               The leading dimension of the array A.       X       (input) COMPLEX*16 array, dimension (LDX,NRHS)               If TRANS = 'N', the N-by-NRHS matrix X.               IF TRANS = 'C', the M-by-NRHS matrix X.       LDX     (input) INTEGER               The leading dimension of the array X.       WORK    (workspace) COMPLEX*16 array dimension (LWORK)       LWORK   (input) INTEGER               length of workspace array required               If TRANS = 'N', LWORK >= (M+NRHS)*(N+2);               if TRANS = 'C', LWORK >= (N+NRHS)*(M+2).       =====================================================================          Parameter adjustments */    a_dim1 = *lda;//.........这里部分代码省略.........

/* Subroutine */ int zgetc2_(integer *n, doublecomplex *a, integer *lda,                             integer *ipiv, integer *jpiv, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    static integer i__, j, ip, jp;    static doublereal eps;    static integer ipv, jpv;    static doublereal smin, xmax;    extern /* Subroutine */ int zgeru_(integer *, integer *, doublecomplex *,                                       doublecomplex *, integer *, doublecomplex *, integer *,                                       doublecomplex *, integer *), zswap_(integer *, doublecomplex *,                                               integer *, doublecomplex *, integer *), dlabad_(doublereal *,                                                       doublereal *);    extern doublereal dlamch_(char *, ftnlen);    static doublereal bignum, smlnum;    /*  -- LAPACK auxiliary 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 */    /*  ======= */    /*  ZGETC2 computes an LU factorization, using complete pivoting, of the */    /*  n-by-n matrix A. The factorization has the form A = P * L * U * Q, */    /*  where P and Q are permutation matrices, L is lower triangular with */    /*  unit diagonal elements and U is upper triangular. */    /*  This is a level 1 BLAS version of the algorithm. */    /*  Arguments */    /*  ========= */    /*  N       (input) INTEGER */    /*          The order of the matrix A. N >= 0. */    /*  A       (input/output) COMPLEX*16 array, dimension (LDA, N) */    /*          On entry, the n-by-n matrix to be factored. */    /*          On exit, the factors L and U from the factorization */    /*          A = P*L*U*Q; the unit diagonal elements of L are not stored. */    /*          If U(k, k) appears to be less than SMIN, U(k, k) is given the */    /*          value of SMIN, giving a nonsingular perturbed system. */    /*  LDA     (input) INTEGER */    /*          The leading dimension of the array A.  LDA >= max(1, N). */    /*  IPIV    (output) INTEGER array, dimension (N). */    /*          The pivot indices; for 1 <= i <= N, row i of the */    /*          matrix has been interchanged with row IPIV(i). */    /*  JPIV    (output) INTEGER array, dimension (N). */    /*          The pivot indices; for 1 <= j <= N, column j of the */    /*          matrix has been interchanged with column JPIV(j). */    /*  INFO    (output) INTEGER */    /*           = 0: successful exit */    /*           > 0: if INFO = k, U(k, k) is likely to produce overflow if */    /*                one tries to solve for x in Ax = b. So U is perturbed */    /*                to avoid the overflow. */    /*  Further Details */    /*  =============== */    /*  Based on contributions by */    /*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */    /*     Umea University, S-901 87 Umea, Sweden. */    /*  ===================================================================== */    /*     .. Parameters .. */    /*     .. */    /*     .. Local Scalars .. */    /*     .. */    /*     .. External Subroutines .. */    /*     .. */    /*     .. External Functions .. */    /*     .. */    /*     .. Intrinsic Functions .. */    /*     .. */    /*     .. Executable Statements .. */    /*     Set constants to control overflow */    /* Parameter adjustments *///.........这里部分代码省略.........

/* Subroutine */ int zlacn2_(integer *n, doublecomplex *v, doublecomplex *x, 	doublereal *est, integer *kase, integer *isave){    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *), d_imag(doublecomplex *);    /* Local variables */    integer i__;    doublereal temp, absxi;    integer jlast;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern integer izmax1_(integer *, doublecomplex *, integer *);    extern doublereal dzsum1_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    doublereal safmin, altsgn, estold;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLACN2 estimates the 1-norm of a square, complex matrix A. *//*  Reverse communication is used for evaluating matrix-vector products. *//*  Arguments *//*  ========= *//*  N      (input) INTEGER *//*         The order of the matrix.  N >= 1. *//*  V      (workspace) COMPLEX*16 array, dimension (N) *//*         On the final return, V = A*W,  where  EST = norm(V)/norm(W) *//*         (W is not returned). *//*  X      (input/output) COMPLEX*16 array, dimension (N) *//*         On an intermediate return, X should be overwritten by *//*               A * X,   if KASE=1, *//*               A' * X,  if KASE=2, *//*         where A' is the conjugate transpose of A, and ZLACN2 must be *//*         re-called with all the other parameters unchanged. *//*  EST    (input/output) DOUBLE PRECISION *//*         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be *//*         unchanged from the previous call to ZLACN2. *//*         On exit, EST is an estimate (a lower bound) for norm(A). *//*  KASE   (input/output) INTEGER *//*         On the initial call to ZLACN2, KASE should be 0. *//*         On an intermediate return, KASE will be 1 or 2, indicating *//*         whether X should be overwritten by A * X  or A' * X. *//*         On the final return from ZLACN2, KASE will again be 0. *//*  ISAVE  (input/output) INTEGER array, dimension (3) *//*         ISAVE is used to save variables between calls to ZLACN2 *//*  Further Details *//*  ======= ======= *//*  Contributed by Nick Higham, University of Manchester. *//*  Originally named CONEST, dated March 16, 1988. *//*  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of *//*  a real or complex matrix, with applications to condition estimation", *//*  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. *//*  Last modified:  April, 1999 *//*  This is a thread safe version of ZLACON, which uses the array ISAVE *//*  in place of a SAVE statement, as follows: *//*     ZLACON     ZLACN2 *//*      JUMP     ISAVE(1) *//*      J        ISAVE(2) *//*      ITER     ISAVE(3) *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *///.........这里部分代码省略.........

//.........这里部分代码省略.........	if (*ku > 0 && *m < *n) {/*           Annihilate a(m,m+1) by applying plane rotations from the *//*           right */	    i__1 = *ku + (*m + 1) * ab_dim1;	    rb.r = ab[i__1].r, rb.i = ab[i__1].i;	    for (i__ = *m; i__ >= 1; --i__) {		zlartg_(&ab[*ku + 1 + i__ * ab_dim1], &rb, &rc, &rs, &ra);		i__1 = *ku + 1 + i__ * ab_dim1;		ab[i__1].r = ra.r, ab[i__1].i = ra.i;		if (i__ > 1) {		    d_cnjg(&z__3, &rs);		    z__2.r = -z__3.r, z__2.i = -z__3.i;		    i__1 = *ku + i__ * ab_dim1;		    z__1.r = z__2.r * ab[i__1].r - z__2.i * ab[i__1].i, 			    z__1.i = z__2.r * ab[i__1].i + z__2.i * ab[i__1]			    .r;		    rb.r = z__1.r, rb.i = z__1.i;		    i__1 = *ku + i__ * ab_dim1;		    i__2 = *ku + i__ * ab_dim1;		    z__1.r = rc * ab[i__2].r, z__1.i = rc * ab[i__2].i;		    ab[i__1].r = z__1.r, ab[i__1].i = z__1.i;		}		if (wantpt) {		    d_cnjg(&z__1, &rs);		    zrot_(n, &pt[i__ + pt_dim1], ldpt, &pt[*m + 1 + pt_dim1], 			    ldpt, &rc, &z__1);		}	    }	}    }/*     Make diagonal and superdiagonal elements real, storing them in D *//*     and E */    i__1 = *ku + 1 + ab_dim1;    t.r = ab[i__1].r, t.i = ab[i__1].i;    i__1 = minmn;    for (i__ = 1; i__ <= i__1; ++i__) {	abst = z_abs(&t);	d__[i__] = abst;	if (abst != 0.) {	    z__1.r = t.r / abst, z__1.i = t.i / abst;	    t.r = z__1.r, t.i = z__1.i;	} else {	    t.r = 1., t.i = 0.;	}	if (wantq) {	    zscal_(m, &t, &q[i__ * q_dim1 + 1], &c__1);	}	if (wantc) {	    d_cnjg(&z__1, &t);	    zscal_(ncc, &z__1, &c__[i__ + c_dim1], ldc);	}	if (i__ < minmn) {	    if (*ku == 0 && *kl == 0) {		e[i__] = 0.;		i__2 = (i__ + 1) * ab_dim1 + 1;		t.r = ab[i__2].r, t.i = ab[i__2].i;	    } else {		if (*ku == 0) {		    i__2 = i__ * ab_dim1 + 2;		    d_cnjg(&z__2, &t);		    z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, 			    z__1.i = ab[i__2].r * z__2.i + ab[i__2].i * 			    z__2.r;		    t.r = z__1.r, t.i = z__1.i;		} else {		    i__2 = *ku + (i__ + 1) * ab_dim1;		    d_cnjg(&z__2, &t);		    z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, 			    z__1.i = ab[i__2].r * z__2.i + ab[i__2].i * 			    z__2.r;		    t.r = z__1.r, t.i = z__1.i;		}		abst = z_abs(&t);		e[i__] = abst;		if (abst != 0.) {		    z__1.r = t.r / abst, z__1.i = t.i / abst;		    t.r = z__1.r, t.i = z__1.i;		} else {		    t.r = 1., t.i = 0.;		}		if (wantpt) {		    zscal_(n, &t, &pt[i__ + 1 + pt_dim1], ldpt);		}		i__2 = *ku + 1 + (i__ + 1) * ab_dim1;		d_cnjg(&z__2, &t);		z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, z__1.i = 			ab[i__2].r * z__2.i + ab[i__2].i * z__2.r;		t.r = z__1.r, t.i = z__1.i;	    }	}    }    return 0;/*     End of ZGBBRD */} /* zgbbrd_ */

/* Subroutine */ int zptrfs_(char *uplo, integer *n, integer *nrhs, 	doublereal *d__, doublecomplex *e, doublereal *df, doublecomplex *ef, 	doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, 	doublereal *ferr, doublereal *berr, doublecomplex *work, doublereal *	rwork, integer *info){    /* System generated locals */    integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5, 	    i__6;    doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8, d__9, d__10, 	    d__11, d__12;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    /* Local variables */    integer i__, j;    doublereal s;    doublecomplex bi, cx, dx, ex;    integer ix, nz;    doublereal eps, safe1, safe2;    extern logical lsame_(char *, char *);    integer count;    logical upper;    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern integer idamax_(integer *, doublereal *, integer *);    doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    doublereal lstres;    extern /* Subroutine */ int zpttrs_(char *, integer *, integer *, 	    doublereal *, doublecomplex *, doublecomplex *, integer *, 	    integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZPTRFS improves the computed solution to a system of linear *//*  equations when the coefficient matrix is Hermitian positive definite *//*  and tridiagonal, and provides error bounds and backward error *//*  estimates for the solution. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the superdiagonal or the subdiagonal of the *//*          tridiagonal matrix A is stored and the form of the *//*          factorization: *//*          = 'U':  E is the superdiagonal of A, and A = U**H*D*U; *//*          = 'L':  E is the subdiagonal of A, and A = L*D*L**H. *//*          (The two forms are equivalent if A is real.) *//*  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. *//*  D       (input) DOUBLE PRECISION array, dimension (N) *//*          The n real diagonal elements of the tridiagonal matrix A. *//*  E       (input) COMPLEX*16 array, dimension (N-1) *//*          The (n-1) off-diagonal elements of the tridiagonal matrix A *//*          (see UPLO). *//*  DF      (input) DOUBLE PRECISION array, dimension (N) *//*          The n diagonal elements of the diagonal matrix D from *//*          the factorization computed by ZPTTRF. *//*  EF      (input) COMPLEX*16 array, dimension (N-1) *//*          The (n-1) off-diagonal elements of the unit bidiagonal *//*          factor U or L from the factorization computed by ZPTTRF *//*          (see UPLO). *//*  B       (input) COMPLEX*16 array, dimension (LDB,NRHS) *//*          The right hand side matrix B. *//*  LDB     (input) INTEGER *//*          The leading dimension of the array B.  LDB >= max(1,N). *//*  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS) *//*          On entry, the solution matrix X, as computed by ZPTTRS. *//*          On exit, the improved solution matrix X. *///.........这里部分代码省略.........

doublereal zlanhp_(char *norm, char *uplo, integer *n, doublecomplex *ap, 	doublereal *work){    /* System generated locals */    integer i__1, i__2;    doublereal ret_val, d__1, d__2, d__3;    /* Builtin functions */    double z_abs(doublecomplex *), sqrt(doublereal);    /* Local variables */    integer i__, j, k;    doublereal sum, absa, scale;    extern logical lsame_(char *, char *);    doublereal value;    extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *, 	     doublereal *, doublereal *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLANHP  returns the value of the one norm,  or the Frobenius norm, or *//*  the  infinity norm,  or the  element of  largest absolute value  of a *//*  complex hermitian matrix A,  supplied in packed form. *//*  Description *//*  =========== *//*  ZLANHP returns the value *//*     ZLANHP = ( max(abs(A(i,j))), NORM = 'M' or 'm' *//*              ( *//*              ( norm1(A),         NORM = '1', 'O' or 'o' *//*              ( *//*              ( normI(A),         NORM = 'I' or 'i' *//*              ( *//*              ( normF(A),         NORM = 'F', 'f', 'E' or 'e' *//*  where  norm1  denotes the  one norm of a matrix (maximum column sum), *//*  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and *//*  normF  denotes the  Frobenius norm of a matrix (square root of sum of *//*  squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm. *//*  Arguments *//*  ========= *//*  NORM    (input) CHARACTER*1 *//*          Specifies the value to be returned in ZLANHP as described *//*          above. *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the upper or lower triangular part of the *//*          hermitian matrix A is supplied. *//*          = 'U':  Upper triangular part of A is supplied *//*          = 'L':  Lower triangular part of A is supplied *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0.  When N = 0, ZLANHP is *//*          set to zero. *//*  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          The upper or lower triangle of the hermitian matrix A, packed *//*          columnwise in a linear array.  The j-th column of A is stored *//*          in the array AP as follows: *//*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; *//*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. *//*          Note that the  imaginary parts of the diagonal elements need *//*          not be set and are assumed to be zero. *//*  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), *//*          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, *//*          WORK is not referenced. *//* ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    --work;    --ap;//.........这里部分代码省略.........

/* Subroutine */int zhptri_(char *uplo, integer *n, doublecomplex *ap, integer *ipiv, doublecomplex *work, integer *info){    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1, z__2;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal d__;    integer j, k;    doublereal t, ak;    integer kc, kp, kx, kpc, npp;    doublereal akp1;    doublecomplex temp, akkp1;    extern logical lsame_(char *, char *);    extern /* Double Complex */    VOID zdotc_f2c_(doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *);    integer kstep;    logical upper;    extern /* Subroutine */    int zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, integer *), zhpmv_(char *, integer *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *), zswap_( integer *, doublecomplex *, integer *, doublecomplex *, integer *) , xerbla_(char *, integer *);    integer kcnext;    /* -- LAPACK computational 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 */    --work;    --ipiv;    --ap;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    if (! upper && ! lsame_(uplo, "L"))    {        *info = -1;    }    else if (*n < 0)    {        *info = -2;    }    if (*info != 0)    {        i__1 = -(*info);        xerbla_("ZHPTRI", &i__1);        return 0;    }    /* Quick return if possible */    if (*n == 0)    {        return 0;    }    /* Check that the diagonal matrix D is nonsingular. */    if (upper)    {        /* Upper triangular storage: examine D from bottom to top */        kp = *n * (*n + 1) / 2;        for (*info = *n;                *info >= 1;                --(*info))        {            i__1 = kp;            if (ipiv[*info] > 0 && (ap[i__1].r == 0. && ap[i__1].i == 0.))            {                return 0;            }            kp -= *info;            /* L10: */        }    }    else    {        /* Lower triangular storage: examine D from top to bottom. */        kp = 1;        i__1 = *n;        for (*info = 1;                *info <= i__1;                ++(*info))        {            i__2 = kp;//.........这里部分代码省略.........

/* Subroutine */ int zlatm6_(integer *type__, integer *n, doublecomplex *a, 	integer *lda, doublecomplex *b, doublecomplex *x, integer *ldx, 	doublecomplex *y, integer *ldy, doublecomplex *alpha, doublecomplex *	beta, doublecomplex *wx, doublecomplex *wy, doublereal *s, doublereal 	*dif){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, y_dim1, 	    y_offset, i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *), sqrt(doublereal);    /* Local variables */    integer i__, j;    doublecomplex z__[64]	/* was [8][8] */;    integer info;    doublecomplex work[26];    doublereal rwork[50];    extern /* Subroutine */ int zlakf2_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *), zgesvd_(char *, char *, integer *, 	    integer *, doublecomplex *, integer *, doublereal *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLATM6 generates test matrices for the generalized eigenvalue *//*  problem, their corresponding right and left eigenvector matrices, *//*  and also reciprocal condition numbers for all eigenvalues and *//*  the reciprocal condition numbers of eigenvectors corresponding to *//*  the 1th and 5th eigenvalues. *//*  Test Matrices *//*  ============= *//*  Two kinds of test matrix pairs *//*           (A, B) = inverse(YH) * (Da, Db) * inverse(X) *//*  are used in the tests: *//*  Type 1: *//*     Da = 1+a   0    0    0    0    Db = 1   0   0   0   0 *//*           0   2+a   0    0    0         0   1   0   0   0 *//*           0    0   3+a   0    0         0   0   1   0   0 *//*           0    0    0   4+a   0         0   0   0   1   0 *//*           0    0    0    0   5+a ,      0   0   0   0   1 *//*  and Type 2: *//*     Da = 1+i   0    0       0       0    Db = 1   0   0   0   0 *//*           0   1-i   0       0       0         0   1   0   0   0 *//*           0    0    1       0       0         0   0   1   0   0 *//*           0    0    0 (1+a)+(1+b)i  0         0   0   0   1   0 *//*           0    0    0       0 (1+a)-(1+b)i,   0   0   0   0   1 . *//*  In both cases the same inverse(YH) and inverse(X) are used to compute *//*  (A, B), giving the exact eigenvectors to (A,B) as (YH, X): *//*  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x *//*          0    1   -y    y   -y         0   1   x  -x  -x *//*          0    0    1    0    0         0   0   1   0   0 *//*          0    0    0    1    0         0   0   0   1   0 *//*          0    0    0    0    1,        0   0   0   0   1 , where *//*  a, b, x and y will have all values independently of each other. *//*  Arguments *//*  ========= *//*  TYPE    (input) INTEGER *//*          Specifies the problem type (see futher details). *//*  N       (input) INTEGER *//*          Size of the matrices A and B. *//*  A       (output) COMPLEX*16 array, dimension (LDA, N). *//*          On exit A N-by-N is initialized according to TYPE. *//*  LDA     (input) INTEGER *//*          The leading dimension of A and of B. *//*  B       (output) COMPLEX*16 array, dimension (LDA, N). *//*          On exit B N-by-N is initialized according to TYPE. *//*  X       (output) COMPLEX*16 array, dimension (LDX, N). *//*          On exit X is the N-by-N matrix of right eigenvectors. *///.........这里部分代码省略.........

doublereal zlantb_(char *norm, char *uplo, char *diag, integer *n, integer *k,	 doublecomplex *ab, integer *ldab, doublereal *work){/*  -- LAPACK auxiliary routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          October 31, 1992       Purpose       =======       ZLANTB  returns the value of the one norm,  or the Frobenius norm, or       the  infinity norm,  or the element of  largest absolute value  of an       n by n triangular band matrix A,  with ( k + 1 ) diagonals.       Description       ===========       ZLANTB returns the value          ZLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm'                   (                   ( norm1(A),         NORM = '1', 'O' or 'o'                   (                   ( normI(A),         NORM = 'I' or 'i'                   (                   ( normF(A),         NORM = 'F', 'f', 'E' or 'e'       where  norm1  denotes the  one norm of a matrix (maximum column sum),       normI  denotes the  infinity norm  of a matrix  (maximum row sum) and       normF  denotes the  Frobenius norm of a matrix (square root of sum of       squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.       Arguments       =========       NORM    (input) CHARACTER*1               Specifies the value to be returned in ZLANTB as described               above.       UPLO    (input) CHARACTER*1               Specifies whether the matrix A is upper or lower triangular.               = 'U':  Upper triangular               = 'L':  Lower triangular       DIAG    (input) CHARACTER*1               Specifies whether or not the matrix A is unit triangular.               = 'N':  Non-unit triangular               = 'U':  Unit triangular       N       (input) INTEGER               The order of the matrix A.  N >= 0.  When N = 0, ZLANTB is               set to zero.       K       (input) INTEGER               The number of super-diagonals of the matrix A if UPLO = 'U',               or the number of sub-diagonals of the matrix A if UPLO = 'L'.               K >= 0.       AB      (input) COMPLEX*16 array, dimension (LDAB,N)               The upper or lower triangular band matrix A, stored in the               first k+1 rows of AB.  The j-th column of A is stored               in the j-th column of the array AB as follows:               if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;               if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).               Note that when DIAG = 'U', the elements of the array AB               corresponding to the diagonal elements of the matrix A are               not referenced, but are assumed to be one.       LDAB    (input) INTEGER               The leading dimension of the array AB.  LDAB >= K+1.       WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK),               where LWORK >= N when NORM = 'I'; otherwise, WORK is not               referenced.      =====================================================================          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5;    doublereal ret_val, d__1, d__2;    /* Builtin functions */    double z_abs(doublecomplex *), sqrt(doublereal);    /* Local variables */    static integer i__, j, l;    static doublereal scale;    static logical udiag;    extern logical lsame_(char *, char *);    static doublereal value;    extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *,	     doublereal *, doublereal *);    static doublereal sum;#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)]//.........这里部分代码省略.........

doublereal dzsum1_(integer *n, doublecomplex *cx, integer *incx){/*  -- LAPACK auxiliary routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          October 31, 1992       Purpose       =======       DZSUM1 takes the sum of the absolute values of a complex       vector and returns a double precision result.       Based on DZASUM from the Level 1 BLAS.       The change is to use the 'genuine' absolute value.       Contributed by Nick Higham for use with ZLACON.       Arguments       =========       N       (input) INTEGER               The number of elements in the vector CX.       CX      (input) COMPLEX*16 array, dimension (N)               The vector whose elements will be summed.       INCX    (input) INTEGER               The spacing between successive values of CX.  INCX > 0.       =====================================================================          Parameter adjustments */    /* System generated locals */    integer i__1, i__2;    doublereal ret_val;    /* Builtin functions */    double z_abs(doublecomplex *);    /* Local variables */    static integer i__, nincx;    static doublereal stemp;    --cx;    /* Function Body */    ret_val = 0.;    stemp = 0.;    if (*n <= 0) {	return ret_val;    }    if (*incx == 1) {	goto L20;    }/*     CODE FOR INCREMENT NOT EQUAL TO 1 */    nincx = *n * *incx;    i__1 = nincx;    i__2 = *incx;    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {/*        NEXT LINE MODIFIED. */	stemp += z_abs(&cx[i__]);/* L10: */    }    ret_val = stemp;    return ret_val;/*     CODE FOR INCREMENT EQUAL TO 1 */L20:    i__2 = *n;    for (i__ = 1; i__ <= i__2; ++i__) {/*        NEXT LINE MODIFIED. */	stemp += z_abs(&cx[i__]);/* L30: */    }    ret_val = stemp;    return ret_val;/*     End of DZSUM1 */} /* dzsum1_ */

/* Subroutine */ int zlaror_slu(char *side, char *init, integer *m, integer *n, 	doublecomplex *a, integer *lda, integer *iseed, doublecomplex *x, 	integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static integer kbeg, jcol;    static doublereal xabs;    static integer irow, j;    static doublecomplex csign;    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static integer ixfrm;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static integer itype, nxfrm;    static doublereal xnorm;    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern int input_error(char *, int *);    static doublereal factor;    extern /* Subroutine */ int zlacgv_slu(integer *, doublecomplex *, integer *)	    ;    extern /* Double Complex */ VOID zlarnd_slu(doublecomplex *, integer *, 	    integer *);    extern /* Subroutine */ int zlaset_slu(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static doublecomplex 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       =======          ZLAROR pre- or post-multiplies an M by N matrix A by a random          unitary 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, pp. 403-409 ).          (BLAS-2 version)       Arguments       =========       SIDE   - CHARACTER*1                SIDE specifies whether A is multiplied on the left or right                by U.            SIDE = 'L'   Multiply A on the left (premultiply) by U            SIDE = 'R'   Multiply A on the right (postmultiply) by U*            SIDE = 'C'   Multiply A on the left by U and the right by U*            SIDE = 'T'   Multiply A on the left by U and the right by U'                Not modified.       INIT   - CHARACTER*1                INIT specifies whether or not A should be initialized to                the identity matrix.                   INIT = 'I'   Initialize A to (a section of) the                                identity matrix before applying U.                   INIT = 'N'   No initialization.  Apply U to the                                input matrix A.                INIT = 'I' may be used to generate square (i.e., unitary)                or rectangular orthogonal matrices (orthogonality being                in the sense of ZDOTC):                For square matrices, M=N, and SIDE many be either 'L' or                'R'; the rows will be orthogonal to each other, as will the                columns.                For rectangular matrices where M < N, SIDE = 'R' will                produce a dense matrix whose rows will be orthogonal and                whose columns will not, while SIDE = 'L' will produce a                matrix whose rows will be orthogonal, and whose first M                columns will be orthogonal, the remaining columns being                zero.                For matrices where M > N, just use the previous                explaination, interchanging 'L' and 'R' and "rows" and                "columns".                Not modified.       M      - INTEGER                Number of rows of A. Not modified.       N      - INTEGER   //.........这里部分代码省略.........

//.........这里部分代码省略.........			    zlarnd_(&z__1, &c__3, &iseed[1]);			    z__[i__5].r = z__1.r, z__[i__5].i = z__1.i;/* L40: */			}			i__4 = n + 1 - jc;			zlarfg_(&i__4, &q[jc + jc * q_dim1], &q[jc + 1 + jc * 				q_dim1], &c__1, &work[jc]);			i__4 = (n << 1) + jc;			i__5 = jc + jc * q_dim1;			d__2 = q[i__5].r;			d__1 = d_sign(&c_b29, &d__2);			work[i__4].r = d__1, work[i__4].i = 0.;			i__4 = jc + jc * q_dim1;			q[i__4].r = 1., q[i__4].i = 0.;			i__4 = n + 1 - jc;			zlarfg_(&i__4, &z__[jc + jc * z_dim1], &z__[jc + 1 + 				jc * z_dim1], &c__1, &work[n + jc]);			i__4 = n * 3 + jc;			i__5 = jc + jc * z_dim1;			d__2 = z__[i__5].r;			d__1 = d_sign(&c_b29, &d__2);			work[i__4].r = d__1, work[i__4].i = 0.;			i__4 = jc + jc * z_dim1;			z__[i__4].r = 1., z__[i__4].i = 0.;/* L50: */		    }		    zlarnd_(&z__1, &c__3, &iseed[1]);		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    i__3 = n + n * q_dim1;		    q[i__3].r = 1., q[i__3].i = 0.;		    i__3 = n;		    work[i__3].r = 0., work[i__3].i = 0.;		    i__3 = n * 3;		    d__1 = z_abs(&ctemp);		    z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;		    work[i__3].r = z__1.r, work[i__3].i = z__1.i;		    zlarnd_(&z__1, &c__3, &iseed[1]);		    ctemp.r = z__1.r, ctemp.i = z__1.i;		    i__3 = n + n * z_dim1;		    z__[i__3].r = 1., z__[i__3].i = 0.;		    i__3 = n << 1;		    work[i__3].r = 0., work[i__3].i = 0.;		    i__3 = n << 2;		    d__1 = z_abs(&ctemp);		    z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;		    work[i__3].r = z__1.r, work[i__3].i = z__1.i;/*                 Apply the diagonal matrices */		    i__3 = n;		    for (jc = 1; jc <= i__3; ++jc) {			i__4 = n;			for (jr = 1; jr <= i__4; ++jr) {			    i__5 = jr + jc * a_dim1;			    i__6 = (n << 1) + jr;			    d_cnjg(&z__3, &work[n * 3 + jc]);			    z__2.r = work[i__6].r * z__3.r - work[i__6].i * 				    z__3.i, z__2.i = work[i__6].r * z__3.i + 				    work[i__6].i * z__3.r;			    i__7 = jr + jc * a_dim1;			    z__1.r = z__2.r * a[i__7].r - z__2.i * a[i__7].i, 				    z__1.i = z__2.r * a[i__7].i + z__2.i * a[				    i__7].r;			    a[i__5].r = z__1.r, a[i__5].i = z__1.i;			    i__5 = jr + jc * b_dim1;			    i__6 = (n << 1) + jr;