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

/* Subroutine */ int zrotg_(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 */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    double sqrt(doublereal);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static doublereal norm;    extern doublereal cdabs_(doublecomplex *);    static doublecomplex alpha;    static doublereal scale;    if (cdabs_(ca) != 0.) {	goto L10;    }    *c = 0.;    s->r = 1., s->i = 0.;    ca->r = cb->r, ca->i = cb->i;    goto L20;L10:    scale = cdabs_(ca) + cdabs_(cb);    z__2.r = scale, z__2.i = 0.;    z_div(&z__1, ca, &z__2);/* Computing 2nd power */    d__1 = cdabs_(&z__1);    z__4.r = scale, z__4.i = 0.;    z_div(&z__3, cb, &z__4);/* Computing 2nd power */    d__2 = cdabs_(&z__3);    norm = scale * sqrt(d__1 * d__1 + d__2 * d__2);    d__1 = cdabs_(ca);    z__1.r = ca->r / d__1, z__1.i = ca->i / d__1;    alpha.r = z__1.r, alpha.i = z__1.i;    *c = cdabs_(ca) / norm;    d_cnjg(&z__3, cb);    z__2.r = alpha.r * z__3.r - alpha.i * z__3.i, z__2.i = alpha.r * z__3.i + 	    alpha.i * z__3.r;    z__1.r = z__2.r / norm, z__1.i = z__2.i / norm;    s->r = z__1.r, s->i = z__1.i;    z__1.r = norm * alpha.r, z__1.i = norm * alpha.i;    ca->r = z__1.r, ca->i = z__1.i;L20:    return 0;} /* zrotg_ */

voidpow_zi(dcomplex *p, dcomplex *a, long int *b) 	/* p = a**b  */{	long int n;	double t;	dcomplex x;	n = *b;	p->dreal = 1;	p->dimag = 0;	if(n == 0)		return;	if(n < 0) {		n = -n;		z_div(&x, p, a);	} else {		x.dreal = a->dreal;		x.dimag = a->dimag;	}	for( ; ; ) {		if(n & 01) {			t = p->dreal * x.dreal - p->dimag * x.dimag;			p->dimag = p->dreal * x.dimag + p->dimag * x.dreal;			p->dreal = t;		}		if(n >>= 1) {			t = x.dreal * x.dreal - x.dimag * x.dimag;			x.dimag = 2 * x.dreal * x.dimag;			x.dreal = t;		} else			break;	}

static voidVanVlietResidue(    int whichDeriv,		/* Which derivative are we evaluating */    const doublecomplex* polej,	/* Which pole are we computing for */    const doublecomplex poles[4],				/* Poles of the filter */    double gain,		/* Gain of the filter */    doublecomplex* residue	/* Output: Computed residue */) {    doublecomplex pi;    doublecomplex pj = *polej;    doublecomplex qj;    doublecomplex gz = {1.0, 0.0};    doublecomplex gp = {1.0, 0.0};    doublecomplex temp, temp2;    int i;    z_recip(&qj, &pj);    if (whichDeriv == 1) {	temp.r = (1.0 - qj.r) * gz.r + qj.i * gz.i;	temp.i = (1.0 - qj.r) * gz.i - qj.i * gz.r; /* gz * (1-qj) */	gz = temp;	temp.r = (1.0 + pj.r) * gz.r - pj.i * gz.i;	temp.i = (1.0 + pj.r) * gz.i + pj.i * gz.r; /* gz * (1+pj) */	gz = temp;	temp.r = pj.r * gz.r - pj.i * gz.i; /* gz * pj */	temp.i = pj.r * gz.i + pj.i * gz.r;	gz.r = 0.5 * temp.r; gz.i = 0.5 * temp.i;    } else if (whichDeriv == 2) {	temp.r = (1.0 - qj.r) * gz.r + qj.i * gz.i;	temp.i = (1.0 - qj.r) * gz.i - qj.i * gz.r; /* gz * (1 - qj) */	gz = temp;	temp.r = (1.0 - pj.r) * gz.r + pj.i * gz.i;	temp.i = (1.0 - pj.r) * gz.i - pj.i * gz.r; /* gz * (1 - pj) */	gz.r = -temp.r; gz.i = -temp.i;    }    for (i = 0; i < 4; ++i) {	pi = poles[i];	if ((pi.r != pj.r) || (pi.i != pj.i && pi.i != -pj.i)) {	    temp.r = 1.0 - pi.r * qj.r + pi.i * qj.i;	    temp.i = - pi.r * qj.i - pi.i * qj.r; /* 1 - pi * qj */	    temp2.r = gp.r * temp.r - gp.i * temp.i;	    temp2.i = gp.i * temp.r + gp.r * temp.i; /* gp * (1 - pi * qj) */	    gp = temp2;	}	temp.r = 1.0 - pi.r * pj.r + pi.i * pj.i;	temp.i = -pi.r * pj.i - pi.i * pj.r; /* 1 - pi * pj */	temp2.r = gp.r * temp.r - gp.i * temp.i;	temp2.i = gp.i * temp.r + gp.r * temp.i; /* gp * (1 - pi * pj) */ 	gp = temp2;    }    z_div(&temp, &gz, &gp); /* gz / gp */    residue->r = gain * temp.r;    residue->i = gain * temp.i; /* gain * gz/gp */}

void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) 	/* p = a**b  */#endif{	integer n;	unsigned long u;	double t;	doublecomplex q, x;	static doublecomplex one = {1.0, 0.0};	n = *b;	q.r = 1;	q.i = 0;	if(n == 0)		goto done;	if(n < 0)		{		n = -n;		z_div(&x, &one, a);		}	else		{		x.r = a->r;		x.i = a->i;		}	for(u = n; ; )		{		if(u & 01)			{			t = q.r * x.r - q.i * x.i;			q.i = q.r * x.i + q.i * x.r;			q.r = t;			}		if(u >>= 1)			{			t = x.r * x.r - x.i * x.i;			x.i = 2 * x.r * x.i;			x.r = t;			}		else			break;		}

void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) 	/* p = a**b  */#endif{integer n;double t;doublecomplex x;static doublecomplex one = {1.0, 0.0};n = *b;p->r = 1;p->i = 0;if(n == 0)	return;if(n < 0)	{	n = -n;	z_div(&x, &one, a);	}else	{	x.r = a->r;	x.i = a->i;	}for( ; ; )	{	if(n & 01)		{		t = p->r * x.r - p->i * x.i;		p->i = p->r * x.i + p->i * x.r;		p->r = t;		}	if(n >>= 1)		{		t = x.r * x.r - x.i * x.i;		x.i = 2 * x.r * x.i;		x.r = t;		}	else		break;	}

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);}

//.........这里部分代码省略.........			.r;		z__1.r = z__2.r, z__1.i = z__2.i;		tmp += (d__1 = z__1.r, abs(d__1)) + (d__2 = d_imag(&z__1), 			abs(d__2));	    }	    rwork[i__] = tmp;	    anorm = max(anorm,tmp);	}    }/*     Quick return if possible. */    if (*n == 0) {	ret_val = 1.;	return ret_val;    } else if (anorm == 0.) {	return ret_val;    }/*     Estimate the norm of inv(op(A)). */    ainvnm = 0.;    kase = 0;L10:    zlacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);    if (kase != 0) {	if (kase == 2) {/*           Multiply by R. */	    i__1 = *n;	    for (i__ = 1; i__ <= i__1; ++i__) {		i__2 = i__;		i__3 = i__;		i__4 = i__;		z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 			work[i__3].i;		work[i__2].r = z__1.r, work[i__2].i = z__1.i;	    }	    if (up) {		zhetrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[			1], n, info);	    } else {		zhetrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[			1], n, info);	    }/*           Multiply by inv(X). */	    i__1 = *n;	    for (i__ = 1; i__ <= i__1; ++i__) {		i__2 = i__;		z_div(&z__1, &work[i__], &x[i__]);		work[i__2].r = z__1.r, work[i__2].i = z__1.i;	    }	} else {/*           Multiply by inv(X'). */	    i__1 = *n;	    for (i__ = 1; i__ <= i__1; ++i__) {		i__2 = i__;		z_div(&z__1, &work[i__], &x[i__]);		work[i__2].r = z__1.r, work[i__2].i = z__1.i;	    }	    if (up) {		zhetrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[			1], n, info);	    } else {		zhetrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[			1], n, info);	    }/*           Multiply by R. */	    i__1 = *n;	    for (i__ = 1; i__ <= i__1; ++i__) {		i__2 = i__;		i__3 = i__;		i__4 = i__;		z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 			work[i__3].i;		work[i__2].r = z__1.r, work[i__2].i = z__1.i;	    }	}	goto L10;    }/*     Compute the estimate of the reciprocal condition number. */    if (ainvnm != 0.) {	ret_val = 1. / ainvnm;    }    return ret_val;} /* zla_hercond_x__ */

/* Subroutine */int zgtsv_(integer *n, integer *nrhs, doublecomplex *dl, doublecomplex *d__, doublecomplex *du, doublecomplex *b, integer *ldb, integer *info){    /* System generated locals */    integer b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4, z__5;    /* Builtin functions */    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer j, k;    doublecomplex temp, mult;    extern /* Subroutine */    int xerbla_(char *, integer *);    /* -- LAPACK driver 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 .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Statement Functions .. */    /* .. */    /* .. Statement Function definitions .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    --dl;    --d__;    --du;    b_dim1 = *ldb;    b_offset = 1 + b_dim1;    b -= b_offset;    /* Function Body */    *info = 0;    if (*n < 0)    {        *info = -1;    }    else if (*nrhs < 0)    {        *info = -2;    }    else if (*ldb < max(1,*n))    {        *info = -7;    }    if (*info != 0)    {        i__1 = -(*info);        xerbla_("ZGTSV ", &i__1);        return 0;    }    if (*n == 0)    {        return 0;    }    i__1 = *n - 1;    for (k = 1;            k <= i__1;            ++k)    {        i__2 = k;        if (dl[i__2].r == 0. && dl[i__2].i == 0.)        {            /* Subdiagonal is zero, no elimination is required. */            i__2 = k;            if (d__[i__2].r == 0. && d__[i__2].i == 0.)            {                /* Diagonal is zero: set INFO = K and return;                a unique */                /* solution can not be found. */                *info = k;                return 0;            }        }        else /* if(complicated condition) */        {            i__2 = k;            i__3 = k;            if ((d__1 = d__[i__2].r, f2c_abs(d__1)) + (d__2 = d_imag(&d__[k]), f2c_abs(d__2)) >= (d__3 = dl[i__3].r, f2c_abs(d__3)) + (d__4 = d_imag(&dl[k]), f2c_abs(d__4)))            {                /* No row interchange required */                z_div(&z__1, &dl[k], &d__[k]);                mult.r = z__1.r;                mult.i = z__1.i; // , expr subst                i__2 = k + 1;                i__3 = k + 1;                i__4 = k;//.........这里部分代码省略.........

/* 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 zdrgvx_(integer *nsize, doublereal *thresh, integer *nin, 	 integer *nout, doublecomplex *a, integer *lda, doublecomplex *b, 	doublecomplex *ai, doublecomplex *bi, doublecomplex *alpha, 	doublecomplex *beta, doublecomplex *vl, doublecomplex *vr, integer *	ilo, integer *ihi, doublereal *lscale, doublereal *rscale, doublereal 	*s, doublereal *dtru, doublereal *dif, doublereal *diftru, 	doublecomplex *work, integer *lwork, doublereal *rwork, integer *	iwork, integer *liwork, doublereal *result, logical *bwork, integer *	info){    /* Format strings */    static char fmt_9999[] = "(/002 ZDRGVX: /002,a,/002 returned INFO=/002,i"	    "6,/002./002,/9x,/002N=/002,i6,/002, JTYPE=/002,i6,/002)/002)";    static char fmt_9998[] = "(/002 ZDRGVX: /002,a,/002 Eigenvectors from"	    " /002,a,/002 incorrectly /002,/002normalized./002,//002 Bits of "	    "error=/002,0p,g10.3,/002,/002,9x,/002N=/002,i6,/002, JTYPE=/002,"	    "i6,/002, IWA=/002,i5,/002, IWB=/002,i5,/002, IWX=/002,i5,/002, I"	    "WY=/002,i5)";    static char fmt_9997[] = "(/1x,a3,/002 -- Complex Expert Eigenvalue/vect"	    "or/002,/002 problem driver/002)";    static char fmt_9995[] = "(/002 Matrix types: /002,/)";    static char fmt_9994[] = "(/002 TYPE 1: Da is diagonal, Db is identity,"	    " /002,//002     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) /002,/"	    "/002     YH and X are left and right eigenvectors. /002,/)";    static char fmt_9993[] = "(/002 TYPE 2: Da is quasi-diagonal, Db is iden"	    "tity, /002,//002     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1)"	    " /002,//002     YH and X are left and right eigenvectors. /002,/)"	    ;    static char fmt_9992[] = "(//002 Tests performed:  /002,/4x,/002 a is al"	    "pha, b is beta, l is a left eigenvector, /002,/4x,/002 r is a ri"	    "ght eigenvector and /002,a,/002 means /002,a,/002./002,//002 1 ="	    " max | ( b A - a B )/002,a,/002 l | / const./002,//002 2 = max |"	    " ( b A - a B ) r | / const./002,//002 3 = max ( Sest/Stru, Stru/"	    "Sest ) /002,/002 over all eigenvalues/002,//002 4 = max( DIFest/"	    "DIFtru, DIFtru/DIFest ) /002,/002 over the 1st and 5th eigenvect"	    "ors/002,/)";    static char fmt_9991[] = "(/002 Type=/002,i2,/002,/002,/002 IWA=/002,i2"	    ",/002, IWB=/002,i2,/002, IWX=/002,i2,/002, IWY=/002,i2,/002, res"	    "ult /002,i2,/002 is/002,0p,f8.2)";    static char fmt_9990[] = "(/002 Type=/002,i2,/002,/002,/002 IWA=/002,i2"	    ",/002, IWB=/002,i2,/002, IWX=/002,i2,/002, IWY=/002,i2,/002, res"	    "ult /002,i2,/002 is/002,1p,d10.3)";    static char fmt_9987[] = "(/002 ZDRGVX: /002,a,/002 returned INFO=/002,i"	    "6,/002./002,/9x,/002N=/002,i6,/002, Input example #/002,i2,/002"	    ")/002)";    static char fmt_9986[] = "(/002 ZDRGVX: /002,a,/002 Eigenvectors from"	    " /002,a,/002 incorrectly /002,/002normalized./002,//002 Bits of "	    "error=/002,0p,g10.3,/002,/002,9x,/002N=/002,i6,/002, Input Examp"	    "le #/002,i2,/002)/002)";    static char fmt_9996[] = "(/002Input Example/002)";    static char fmt_9989[] = "(/002 Input example #/002,i2,/002, matrix orde"	    "r=/002,i4,/002,/002,/002 result /002,i2,/002 is/002,0p,f8.2)";    static char fmt_9988[] = "(/002 Input example #/002,i2,/002, matrix orde"	    "r=/002,i4,/002,/002,/002 result /002,i2,/002 is/002,1p,d10.3)";    /* System generated locals */    integer a_dim1, a_offset, ai_dim1, ai_offset, b_dim1, b_offset, bi_dim1, 	    bi_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),	     s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_rsle(void);    /* Local variables */    integer i__, j, n, iwa, iwb;    doublereal ulp;    integer iwx, iwy, nmax, linfo;    doublereal anorm, bnorm;    extern /* Subroutine */ int zget52_(logical *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, doublecomplex *, doublereal *, 	     doublereal *);    integer nerrs;    doublereal ratio1, ratio2, thrsh2;    extern /* Subroutine */ int zlatm6_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublereal *, doublereal *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    doublereal abnorm;    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *);    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 	    *, integer *);    doublecomplex weight[5];    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    integer minwrk, maxwrk, iptype;    extern /* Subroutine */ int zggevx_(char *, char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *, 	     doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, integer *, integer *, doublereal *, //.........这里部分代码省略.........

/* Subroutine */int zsytf2_rook_(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;    doublereal d__1, d__2;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k, p;    doublecomplex t, d11, d12, d21, d22;    integer ii, kk, kp;    doublecomplex wk, wkm1, wkp1;    logical done;    integer imax, jmax;    extern /* Subroutine */    int zsyr_(char *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *);    doublereal alpha;    extern logical lsame_(char *, char *);    doublereal dtemp, sfmin;    extern /* Subroutine */    int zscal_(integer *, doublecomplex *, doublecomplex *, integer *);    integer itemp, kstep;    logical upper;    extern /* Subroutine */    int zswap_(integer *, doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    doublereal absakk;    extern /* Subroutine */    int xerbla_(char *, integer *);    doublereal colmax;    extern integer izamax_(integer *, doublecomplex *, integer *);    doublereal rowmax;    /* -- LAPACK computational routine (version 3.5.0) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* November 2013 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Statement Functions .. */    /* .. */    /* .. Statement Function definitions .. */    /* .. */    /* .. Executable Statements .. */    /* Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --ipiv;    /* 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_("ZSYTF2_ROOK", &i__1);        return 0;    }    /* Initialize ALPHA for use in choosing pivot block size. */    alpha = (sqrt(17.) + 1.) / 8.;    /* Compute machine safe minimum */    sfmin = dlamch_("S");    if (upper)    {        /* Factorize A as U*D*U**T using the upper triangle of A */        /* K is the main loop index, decreasing from N to 1 in steps of */        /* 1 or 2 */        k = *n;L10: /* If K < 1, exit from loop */        if (k < 1)        {            goto L70;//.........这里部分代码省略.........

/* Subroutine */ int ztptri_(char *uplo, char *diag, integer *n, 	doublecomplex *ap, integer *info){    /* System generated locals */    integer i__1, i__2;    doublecomplex z__1;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer j, jc, jj;    doublecomplex ajj;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    logical upper;    extern /* Subroutine */ int ztpmv_(char *, char *, char *, integer *, 	    doublecomplex *, doublecomplex *, integer *), xerbla_(char *, integer *);    integer jclast;    logical nounit;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTPTRI computes the inverse of a complex upper or lower triangular *//*  matrix A stored in packed format. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  A is upper triangular; *//*          = 'L':  A is lower triangular. *//*  DIAG    (input) CHARACTER*1 *//*          = 'N':  A is non-unit triangular; *//*          = 'U':  A is unit triangular. *//*  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 triangular matrix A, stored *//*          columnwise in a linear array.  The j-th column of A is stored *//*          in the array AP as follows: *//*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; *//*          if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n. *//*          See below for further details. *//*          On exit, the (triangular) inverse of the original matrix, in *//*          the same packed storage format. *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0:  if INFO = -i, the i-th argument had an illegal value *//*          > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular *//*                matrix is singular and its inverse can not be computed. *//*  Further Details *//*  =============== *//*  A triangular matrix A can be transferred to packed storage using one *//*  of the following program segments: *//*  UPLO = 'U':                      UPLO = 'L': *//*        JC = 1                           JC = 1 *//*        DO 2 J = 1, N                    DO 2 J = 1, N *//*           DO 1 I = 1, J                    DO 1 I = J, N *//*              AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J) *//*      1    CONTINUE                    1    CONTINUE *//*           JC = JC + J                      JC = JC + N - J + 1 *//*      2 CONTINUE                       2 CONTINUE *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    --ap;//.........这里部分代码省略.........

 int zrotg_(doublecomplex *ca, doublecomplex *cb, double *	c__, doublecomplex *s){    /* System generated locals */    double d__1, d__2;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    double sqrt(double);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    double norm;    doublecomplex alpha;    double scale;/*     .. Scalar Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*     determines a double complex Givens rotation. *//*     .. Local Scalars .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. */    if (z_abs(ca) != 0.) {	goto L10;    }    *c__ = 0.;    s->r = 1., s->i = 0.;    ca->r = cb->r, ca->i = cb->i;    goto L20;L10:    scale = z_abs(ca) + z_abs(cb);    z__2.r = scale, z__2.i = 0.;    z_div(&z__1, ca, &z__2);/* Computing 2nd power */    d__1 = z_abs(&z__1);    z__4.r = scale, z__4.i = 0.;    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.r = ca->r / d__1, z__1.i = ca->i / d__1;    alpha.r = z__1.r, alpha.i = z__1.i;    *c__ = z_abs(ca) / norm;    d_cnjg(&z__3, cb);    z__2.r = alpha.r * z__3.r - alpha.i * z__3.i, z__2.i = alpha.r * z__3.i + 	    alpha.i * z__3.r;    z__1.r = z__2.r / norm, z__1.i = z__2.i / norm;    s->r = z__1.r, s->i = z__1.i;    z__1.r = norm * alpha.r, z__1.i = norm * alpha.i;    ca->r = z__1.r, ca->i = z__1.i;L20:    return 0;} /* zrotg_ */

//.........这里部分代码省略.........	    if (ipiv[*info] > 0 && (ap[i__2].r == 0. && ap[i__2].i == 0.)) {		return 0;	    }	    kp = kp + *n - *info + 1;	}    }    *info = 0;    if (upper) {/*        Compute inv(A) from the factorization A = U*D*U'. *//*        K is the main loop index, increasing from 1 to N in steps of *//*        1 or 2, depending on the size of the diagonal blocks. */	k = 1;	kc = 1;L30:/*        If K > N, exit from loop. */	if (k > *n) {	    goto L50;	}	kcnext = kc + k;	if (ipiv[k] > 0) {/*           1 x 1 diagonal block *//*           Invert the diagonal block. */	    i__1 = kc + k - 1;	    z_div(&z__1, &c_b1, &ap[kc + k - 1]);	    ap[i__1].r = z__1.r, ap[i__1].i = z__1.i;/*           Compute column K of the inverse. */	    if (k > 1) {		i__1 = k - 1;		zcopy_(&i__1, &ap[kc], &c__1, &work[1], &c__1);		i__1 = k - 1;		z__1.r = -1., z__1.i = -0.;		zspmv_(uplo, &i__1, &z__1, &ap[1], &work[1], &c__1, &c_b2, &			ap[kc], &c__1);		i__1 = kc + k - 1;		i__2 = kc + k - 1;		i__3 = k - 1;		zdotu_(&z__2, &i__3, &work[1], &c__1, &ap[kc], &c__1);		z__1.r = ap[i__2].r - z__2.r, z__1.i = ap[i__2].i - z__2.i;		ap[i__1].r = z__1.r, ap[i__1].i = z__1.i;	    }	    kstep = 1;	} else {/*           2 x 2 diagonal block *//*           Invert the diagonal block. */	    i__1 = kcnext + k - 1;	    t.r = ap[i__1].r, t.i = ap[i__1].i;	    z_div(&z__1, &ap[kc + k - 1], &t);	    ak.r = z__1.r, ak.i = z__1.i;	    z_div(&z__1, &ap[kcnext + k], &t);	    akp1.r = z__1.r, akp1.i = z__1.i;	    z_div(&z__1, &ap[kcnext + k - 1], &t);

/* Subroutine */ int zlaghe_(integer *n, integer *k, doublereal *d, 	doublecomplex *a, integer *lda, integer *iseed, doublecomplex *work, 	integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(	    doublecomplex *, doublecomplex *);    /* Local variables */    extern /* Subroutine */ int zher2_(char *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    static integer i, j;    static doublecomplex alpha;    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    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 *), 	    zhemv_(char *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *), zaxpy_(integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static doublecomplex wa, wb;    static doublereal wn;    extern /* Subroutine */ int xerbla_(char *, integer *), zlarnv_(	    integer *, integer *, integer *, doublecomplex *);    static doublecomplex tau;/*  -- 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       =======       ZLAGHE generates a complex hermitian matrix A, by pre- and post-       multiplying a real diagonal matrix D with a random unitary matrix:       A = U*D*U'. The semi-bandwidth may then be reduced to k by additional       unitary transformations.       Arguments       =========       N       (input) INTEGER               The order of the matrix A.  N >= 0.       K       (input) INTEGER               The number of nonzero subdiagonals within the band of A.               0 <= K <= N-1.       D       (input) DOUBLE PRECISION array, dimension (N)               The diagonal elements of the diagonal matrix D.       A       (output) COMPLEX*16 array, dimension (LDA,N)               The generated n by n hermitian matrix A (the full matrix is               stored).       LDA     (input) INTEGER               The leading dimension of the array A.  LDA >= N.       ISEED   (input/output) INTEGER array, dimension (4)               On entry, the seed of the random number generator; the array               elements must be between 0 and 4095, and ISEED(4) must be               odd.               On exit, the seed is updated.       WORK    (workspace) COMPLEX*16 array, dimension (2*N)       INFO    (output) INTEGER               = 0: successful exit               < 0: if INFO = -i, the i-th argument had an illegal value       =====================================================================          Test the input arguments          Parameter adjustments */    --d;    a_dim1 = *lda;    a_offset = a_dim1 + 1;//.........这里部分代码省略.........

//.........这里部分代码省略.........            If IPIV(k) > 0, then rows and columns k and IPIV(k) were               interchanged and D(k,k) is a 1-by-1 diagonal block.               If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and               columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)               is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =               IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were               interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.       W       (workspace) COMPLEX*16 array, dimension (LDW,NB)       LDW     (input) INTEGER               The leading dimension of the array W.  LDW >= max(1,N).       INFO    (output) INTEGER               = 0: successful exit               > 0: if INFO = k, D(k,k) is exactly zero.  The factorization                    has been completed, but the block diagonal matrix D is                    exactly singular.       =====================================================================          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

doublereal zla_gbrcond_x__(char *trans, integer *n, integer *kl, integer *ku, 	doublecomplex *ab, integer *ldab, doublecomplex *afb, integer *ldafb, 	integer *ipiv, doublecomplex *x, integer *info, doublecomplex *work, 	doublereal *rwork, ftnlen trans_len){    /* System generated locals */    integer ab_dim1, ab_offset, afb_dim1, afb_offset, i__1, i__2, i__3, i__4;    doublereal ret_val, d__1, d__2;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, kd, ke;    doublereal tmp;    integer kase;    extern logical lsame_(char *, char *);    integer isave[3];    doublereal anorm;    extern /* Subroutine */ int zlacn2_(integer *, doublecomplex *, 	    doublecomplex *, doublereal *, integer *, integer *), xerbla_(	    char *, integer *);    doublereal ainvnm;    extern /* Subroutine */ int zgbtrs_(char *, integer *, integer *, integer 	    *, integer *, doublecomplex *, integer *, integer *, 	    doublecomplex *, integer *, integer *);    logical notrans;/*     -- LAPACK routine (version 3.2.1)                               -- *//*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- *//*     -- Jason Riedy of Univ. of California Berkeley.                 -- *//*     -- April 2009                                                   -- *//*     -- LAPACK is a software package provided by Univ. of Tennessee, -- *//*     -- Univ. of California Berkeley and NAG Ltd.                    -- *//*     .. *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*  Purpose *//*  ======= *//*     ZLA_GBRCOND_X Computes the infinity norm condition number of *//*     op(A) * diag(X) where X is a COMPLEX*16 vector. *//*  Arguments *//*  ========= *//*     TRANS   (input) CHARACTER*1 *//*     Specifies the form of the system of equations: *//*       = 'N':  A * X = B     (No transpose) *//*       = 'T':  A**T * X = B  (Transpose) *//*       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose) *//*     N       (input) INTEGER *//*     The number of linear equations, i.e., the order of the *//*     matrix A.  N >= 0. *//*     KL      (input) INTEGER *//*     The number of subdiagonals within the band of A.  KL >= 0. *//*     KU      (input) INTEGER *//*     The number of superdiagonals within the band of A.  KU >= 0. *//*     AB      (input) COMPLEX*16 array, dimension (LDAB,N) *//*     On entry, the matrix A in band storage, in rows 1 to KL+KU+1. *//*     The j-th column of A is stored in the j-th column of the *//*     array AB as follows: *//*     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) *//*     LDAB    (input) INTEGER *//*     The leading dimension of the array AB.  LDAB >= KL+KU+1. *//*     AFB     (input) COMPLEX*16 array, dimension (LDAFB,N) *//*     Details of the LU factorization of the band matrix A, as *//*     computed by ZGBTRF.  U is stored as an upper triangular *//*     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, *//*     and the multipliers used during the factorization are stored *//*     in rows KL+KU+2 to 2*KL+KU+1. *//*     LDAFB   (input) INTEGER *//*     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1. *//*     IPIV    (input) INTEGER array, dimension (N) *//*     The pivot indices from the factorization A = P*L*U *//*     as computed by ZGBTRF; row i of the matrix was interchanged *//*     with row IPIV(i). *//*     X       (input) COMPLEX*16 array, dimension (N) *//*     The vector X in the formula op(A) * diag(X). *//*     INFO    (output) INTEGER *//*       = 0:  Successful exit. *//*     i > 0:  The ith argument is invalid. *///.........这里部分代码省略.........

/* Subroutine */ int zlatdf_(integer *ijob, integer *n, doublecomplex *z__, 	integer *ldz, doublecomplex *rhs, doublereal *rdsum, doublereal *	rdscal, integer *ipiv, integer *jpiv){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    void z_sqrt(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k;    doublecomplex bm, bp, xm[2], xp[2];    integer info;    doublecomplex temp, work[8];    doublereal scale;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    doublecomplex pmone;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    doublereal rtemp, sminu, rwork[2];    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    doublereal splus;    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *), zgesc2_(	    integer *, doublecomplex *, integer *, doublecomplex *, integer *, 	     integer *, doublereal *), zgecon_(char *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *, 	     doublereal *, doublereal *), zlaswp_(integer *, doublecomplex *, 	    integer *, integer *, integer *, integer *, integer *);/*  -- LAPACK auxiliary routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLATDF computes the contribution to the reciprocal Dif-estimate *//*  by solving for x in Z * x = b, where b is chosen such that the norm *//*  of x is as large as possible. It is assumed that LU decomposition *//*  of Z has been computed by ZGETC2. On entry RHS = f holds the *//*  contribution from earlier solved sub-systems, and on return RHS = x. *//*  The factorization of Z returned by ZGETC2 has the form *//*  Z = P * L * U * Q, where P and Q are permutation matrices. L is lower *//*  triangular with unit diagonal elements and U is upper triangular. *//*  Arguments *//*  ========= *//*  IJOB    (input) INTEGER *//*          IJOB = 2: First compute an approximative null-vector e *//*              of Z using ZGECON, e is normalized and solve for *//*              Zx = +-e - f with the sign giving the greater value of *//*              2-norm(x).  About 5 times as expensive as Default. *//*          IJOB .ne. 2: Local look ahead strategy where *//*              all entries of the r.h.s. b is choosen as either +1 or *//*              -1.  Default. *//*  N       (input) INTEGER *//*          The number of columns of the matrix Z. *//*  Z       (input) DOUBLE PRECISION array, dimension (LDZ, N) *//*          On entry, the LU part of the factorization of the n-by-n *//*          matrix Z computed by ZGETC2:  Z = P * L * U * Q *//*  LDZ     (input) INTEGER *//*          The leading dimension of the array Z.  LDA >= max(1, N). *//*  RHS     (input/output) DOUBLE PRECISION array, dimension (N). *//*          On entry, RHS contains contributions from other subsystems. *//*          On exit, RHS contains the solution of the subsystem with *//*          entries according to the value of IJOB (see above). *//*  RDSUM   (input/output) DOUBLE PRECISION *//*          On entry, the sum of squares of computed contributions to *//*          the Dif-estimate under computation by ZTGSYL, where the *//*          scaling factor RDSCAL (see below) has been factored out. *//*          On exit, the corresponding sum of squares updated with the *//*          contributions from the current sub-system. *//*          If TRANS = 'T' RDSUM is not touched. *//*          NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL. *//*  RDSCAL  (input/output) DOUBLE PRECISION *//*          On entry, scaling factor used to prevent overflow in RDSUM. *///.........这里部分代码省略.........

/* 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   //.........这里部分代码省略.........

/* Subroutine */ int zsptrf_(char *uplo, integer *n, doublecomplex *ap, 	integer *ipiv, integer *info){    /* System generated locals */    integer 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;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k;    doublecomplex t, r1, d11, d12, d21, d22;    integer kc, kk, kp;    doublecomplex wk;    integer kx, knc, kpc, npp;    doublecomplex wkm1, wkp1;    integer imax, jmax;    extern /* Subroutine */ int zspr_(char *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *);    doublereal alpha;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    integer kstep;    logical upper;    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    doublereal absakk;    extern /* Subroutine */ int xerbla_(char *, 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 *//*  ======= *//*  ZSPTRF computes the factorization of a complex symmetric matrix A *//*  stored in packed format using the Bunch-Kaufman diagonal pivoting *//*  method: *//*     A = U*D*U**T  or  A = L*D*L**T *//*  where U (or L) is a product of permutation and unit upper (lower) *//*  triangular matrices, and D is symmetric and block diagonal with *//*  1-by-1 and 2-by-2 diagonal blocks. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  Upper triangle of A is stored; *//*          = 'L':  Lower triangle of A is stored. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          On entry, the upper or lower triangle of the symmetric matrix *//*          A, packed columnwise in a linear array.  The j-th column of A *//*          is stored in the array AP as follows: *//*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; *//*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. *//*          On exit, the block diagonal matrix D and the multipliers used *//*          to obtain the factor U or L, stored as a packed triangular *//*          matrix overwriting A (see below for further details). *//*  IPIV    (output) INTEGER array, dimension (N) *//*          Details of the interchanges and the block structure of D. *//*          If IPIV(k) > 0, then rows and columns k and IPIV(k) were *//*          interchanged and D(k,k) is a 1-by-1 diagonal block. *//*          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and *//*          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *//*          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) = *//*          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were *//*          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -i, the i-th argument had an illegal value *//*          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization *//*               has been completed, but the block diagonal matrix D is *//*               exactly singular, and division by zero will occur if it *//*               is used to solve a system of equations. *//*  Further Details *//*  =============== *///.........这里部分代码省略.........

//.........这里部分代码省略.........    nsupr      = xlsub[fsupc+1] - lptr;    lu_sup_ptr = &lusup[xlusup[fsupc]];	/* start of the current supernode */    lu_col_ptr = &lusup[xlusup[jcol]];	/* start of jcol in the supernode */    lsub_ptr   = &lsub[lptr];	/* start of row indices of the supernode */#ifdef DEBUGif ( jcol == MIN_COL ) {    printf("Before cdiv: col %d/n", jcol);    for (k = nsupc; k < nsupr; k++) 	printf("  lu[%d] %f/n", lsub_ptr[k], lu_col_ptr[k]);}#endif        /* Determine the largest abs numerical value for partial pivoting;       Also search for user-specified pivot, and diagonal element. */    if ( *usepr ) *pivrow = iperm_r[jcol];    diagind = iperm_c[jcol];    pivmax = 0.0;    pivptr = nsupc;    diag = EMPTY;    old_pivptr = nsupc;    for (isub = nsupc; isub < nsupr; ++isub) {        rtemp = z_abs1 (&lu_col_ptr[isub]);	if ( rtemp > pivmax ) {	    pivmax = rtemp;	    pivptr = isub;	}	if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;	if ( lsub_ptr[isub] == diagind ) diag = isub;    }    /* Test for singularity */    if ( pivmax == 0.0 ) {#if 1#if SCIPY_FIX	if (pivptr < nsupr) {	    *pivrow = lsub_ptr[pivptr];	}	else {	    *pivrow = diagind;	}#else	*pivrow = lsub_ptr[pivptr];#endif	perm_r[*pivrow] = jcol;#else	perm_r[diagind] = jcol;#endif	*usepr = 0;	return (jcol+1);    }    thresh = u * pivmax;        /* Choose appropriate pivotal element by our policy. */    if ( *usepr ) {        rtemp = z_abs1 (&lu_col_ptr[old_pivptr]);	if ( rtemp != 0.0 && rtemp >= thresh )	    pivptr = old_pivptr;	else	    *usepr = 0;    }    if ( *usepr == 0 ) {	/* Use diagonal pivot? */	if ( diag >= 0 ) { /* diagonal exists */            rtemp = z_abs1 (&lu_col_ptr[diag]);	    if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;        }	*pivrow = lsub_ptr[pivptr];    }        /* Record pivot row */    perm_r[*pivrow] = jcol;        /* Interchange row subscripts */    if ( pivptr != nsupc ) {	itemp = lsub_ptr[pivptr];	lsub_ptr[pivptr] = lsub_ptr[nsupc];	lsub_ptr[nsupc] = itemp;	/* Interchange numerical values as well, for the whole snode, such 	 * that L is indexed the same way as A. 	 */	for (icol = 0; icol <= nsupc; icol++) {	    itemp = pivptr + icol * nsupr;	    temp = lu_sup_ptr[itemp];	    lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];	    lu_sup_ptr[nsupc + icol*nsupr] = temp;	}    } /* if */    /* cdiv operation */    ops[FACT] += 10 * (nsupr - nsupc);    z_div(&temp, &one, &lu_col_ptr[nsupc]);    for (k = nsupc+1; k < nsupr; k++) 	zz_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);    return 0;}

/* 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 zgttrf_(integer *n, doublecomplex *dl, doublecomplex *	d__, doublecomplex *du, doublecomplex *du2, integer *ipiv, integer *	info){    /* System generated locals */    integer i__1, i__2, i__3, i__4;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    static integer i__;    static doublecomplex fact, temp;    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);/*  -- LAPACK routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     June 30, 1999 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGTTRF computes an LU factorization of a complex tridiagonal matrix A *//*  using elimination with partial pivoting and row interchanges. *//*  The factorization has the form *//*     A = L * U *//*  where L is a product of permutation and unit lower bidiagonal *//*  matrices and U is upper triangular with nonzeros in only the main *//*  diagonal and first two superdiagonals. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the matrix A. *//*  DL      (input/output) COMPLEX*16 array, dimension (N-1) *//*          On entry, DL must contain the (n-1) sub-diagonal elements of *//*          A. *//*          On exit, DL is overwritten by the (n-1) multipliers that *//*          define the matrix L from the LU factorization of A. *//*  D       (input/output) COMPLEX*16 array, dimension (N) *//*          On entry, D must contain the diagonal elements of A. *//*          On exit, D is overwritten by the n diagonal elements of the *//*          upper triangular matrix U from the LU factorization of A. *//*  DU      (input/output) COMPLEX*16 array, dimension (N-1) *//*          On entry, DU must contain the (n-1) super-diagonal elements *//*          of A. *//*          On exit, DU is overwritten by the (n-1) elements of the first *//*          super-diagonal of U. *//*  DU2     (output) COMPLEX*16 array, dimension (N-2) *//*          On exit, DU2 is overwritten by the (n-2) elements of the *//*          second super-diagonal of U. *//*  IPIV    (output) INTEGER array, dimension (N) *//*          The pivot indices; for 1 <= i <= n, row i of the matrix was *//*          interchanged with row IPIV(i).  IPIV(i) will always be either *//*          i or i+1; IPIV(i) = i indicates a row interchange was not *//*          required. *//*  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 Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Statement Functions .. *//*     .. *//*     .. Statement Function definitions .. *//*     .. *//*     .. Executable Statements .. *///.........这里部分代码省略.........

//.........这里部分代码省略.........    --iwork;    --dr;    --dl;    --d__;    --iseed;    /* Function Body */    if (*i__ < 1 || *i__ > *m || *j < 1 || *j > *n) {	 ret_val->r = 0.,  ret_val->i = 0.;	return ;    }/*     Check for banding */    if (*j > *i__ + *ku || *j < *i__ - *kl) {	 ret_val->r = 0.,  ret_val->i = 0.;	return ;    }/*     Check for sparsity */    if (*sparse > 0.) {	if (dlaran_(&iseed[1]) < *sparse) {	     ret_val->r = 0.,  ret_val->i = 0.;	    return ;	}    }/*     Compute subscripts depending on IPVTNG */    if (*ipvtng == 0) {	isub = *i__;	jsub = *j;    } else if (*ipvtng == 1) {	isub = iwork[*i__];	jsub = *j;    } else if (*ipvtng == 2) {	isub = *i__;	jsub = iwork[*j];    } else if (*ipvtng == 3) {	isub = iwork[*i__];	jsub = iwork[*j];    }/*     Compute entry and grade it according to IGRADE */    if (isub == jsub) {	i__1 = isub;	ctemp.r = d__[i__1].r, ctemp.i = d__[i__1].i;    } else {	zlarnd_(&z__1, idist, &iseed[1]);	ctemp.r = z__1.r, ctemp.i = z__1.i;    }    if (*igrade == 1) {	i__1 = isub;	z__1.r = ctemp.r * dl[i__1].r - ctemp.i * dl[i__1].i, z__1.i = 		ctemp.r * dl[i__1].i + ctemp.i * dl[i__1].r;	ctemp.r = z__1.r, ctemp.i = z__1.i;    } else if (*igrade == 2) {	i__1 = jsub;	z__1.r = ctemp.r * dr[i__1].r - ctemp.i * dr[i__1].i, z__1.i = 		ctemp.r * dr[i__1].i + ctemp.i * dr[i__1].r;	ctemp.r = z__1.r, ctemp.i = z__1.i;    } else if (*igrade == 3) {	i__1 = isub;	z__2.r = ctemp.r * dl[i__1].r - ctemp.i * dl[i__1].i, z__2.i = 		ctemp.r * dl[i__1].i + ctemp.i * dl[i__1].r;	i__2 = jsub;	z__1.r = z__2.r * dr[i__2].r - z__2.i * dr[i__2].i, z__1.i = z__2.r * 		dr[i__2].i + z__2.i * dr[i__2].r;	ctemp.r = z__1.r, ctemp.i = z__1.i;    } else if (*igrade == 4 && isub != jsub) {	i__1 = isub;	z__2.r = ctemp.r * dl[i__1].r - ctemp.i * dl[i__1].i, z__2.i = 		ctemp.r * dl[i__1].i + ctemp.i * dl[i__1].r;	z_div(&z__1, &z__2, &dl[jsub]);	ctemp.r = z__1.r, ctemp.i = z__1.i;    } else if (*igrade == 5) {	i__1 = isub;	z__2.r = ctemp.r * dl[i__1].r - ctemp.i * dl[i__1].i, z__2.i = 		ctemp.r * dl[i__1].i + ctemp.i * dl[i__1].r;	d_cnjg(&z__3, &dl[jsub]);	z__1.r = z__2.r * z__3.r - z__2.i * z__3.i, z__1.i = z__2.r * z__3.i 		+ z__2.i * z__3.r;	ctemp.r = z__1.r, ctemp.i = z__1.i;    } else if (*igrade == 6) {	i__1 = isub;	z__2.r = ctemp.r * dl[i__1].r - ctemp.i * dl[i__1].i, z__2.i = 		ctemp.r * dl[i__1].i + ctemp.i * dl[i__1].r;	i__2 = jsub;	z__1.r = z__2.r * dl[i__2].r - z__2.i * dl[i__2].i, z__1.i = z__2.r * 		dl[i__2].i + z__2.i * dl[i__2].r;	ctemp.r = z__1.r, ctemp.i = z__1.i;    }     ret_val->r = ctemp.r,  ret_val->i = ctemp.i;    return ;/*     End of ZLATM2 */} /* zlatm2_ */

/* 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;//.........这里部分代码省略.........

//.........这里部分代码省略.........		       	&nsupr, &x[fsupc], &incx, &beta, &work[0], &incy);#endif#else		    zlsolve ( nsupr, nsupc, &Lval[luptr], &x[fsupc]);				    zmatvec ( nsupr, nsupr-nsupc, nsupc, &Lval[luptr+nsupc],			&x[fsupc], &work[0] );#endif						    iptr = istart + nsupc;		    for (i = 0; i < nrow; ++i, ++iptr) {			irow = L_SUB(iptr);			z_sub(&x[irow], &x[irow], &work[i]); /* Scatter */			work[i] = comp_zero;		    }	 	}	    } /* for k ... */	    	} else {	    /* Form x := inv(U)*x */	    	    if ( U->nrow == 0 ) return 0; /* Quick return */	    	    for (k = Lstore->nsuper; k >= 0; k--) {	    	fsupc = L_FST_SUPC(k);	    	nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);	    	nsupc = L_FST_SUPC(k+1) - fsupc;	    	luptr = L_NZ_START(fsupc);		    	        solve_ops += 4 * nsupc * (nsupc + 1);		if ( nsupc == 1 ) {		    z_div(&x[fsupc], &x[fsupc], &Lval[luptr]);		    for (i = U_NZ_START(fsupc); i < U_NZ_START(fsupc+1); ++i) {			irow = U_SUB(i);			zz_mult(&comp_zero, &x[fsupc], &Uval[i]);			z_sub(&x[irow], &x[irow], &comp_zero);		    }		} else {#ifdef USE_VENDOR_BLAS#ifdef _CRAY		    CTRSV(ftcs3, ftcs2, ftcs2, &nsupc, &Lval[luptr], &nsupr,		       &x[fsupc], &incx);#else		    ztrsv_("U", "N", "N", &nsupc, &Lval[luptr], &nsupr,		       &x[fsupc], &incx);#endif#else				    zusolve ( nsupr, nsupc, &Lval[luptr], &x[fsupc] );#endif				    for (jcol = fsupc; jcol < L_FST_SUPC(k+1); jcol++) {		        solve_ops += 8*(U_NZ_START(jcol+1) - U_NZ_START(jcol));		    	for (i = U_NZ_START(jcol); i < U_NZ_START(jcol+1); 				i++) {			    irow = U_SUB(i);			zz_mult(&comp_zero, &x[jcol], &Uval[i]);			z_sub(&x[irow], &x[irow], &comp_zero);		    	}                    }		}	    } /* for k ... */	    	}    } else { /* Form x := inv(A')*x */

/* Subroutine */ int ztptri_(char *uplo, char *diag, integer *n, 	doublecomplex *ap, integer *info){/*  -- LAPACK 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       =======       ZTPTRI computes the inverse of a complex upper or lower triangular       matrix A stored in packed format.       Arguments       =========       UPLO    (input) CHARACTER*1               = 'U':  A is upper triangular;               = 'L':  A is lower triangular.       DIAG    (input) CHARACTER*1               = 'N':  A is non-unit triangular;               = 'U':  A is unit triangular.       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 triangular matrix A, stored               columnwise in a linear array.  The j-th column of A is stored               in the array AP as follows:               if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;               if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.               See below for further details.               On exit, the (triangular) inverse of the original matrix, in               the same packed storage format.       INFO    (output) INTEGER               = 0:  successful exit               < 0:  if INFO = -i, the i-th argument had an illegal value               > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular                     matrix is singular and its inverse can not be computed.       Further Details       ===============       A triangular matrix A can be transferred to packed storage using one       of the following program segments:       UPLO = 'U':                      UPLO = 'L':             JC = 1                           JC = 1             DO 2 J = 1, N                    DO 2 J = 1, N                DO 1 I = 1, J                    DO 1 I = J, N                   AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)           1    CONTINUE                    1    CONTINUE                JC = JC + J                      JC = JC + N - J + 1           2 CONTINUE                       2 CONTINUE       =====================================================================          Test the input parameters.          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer i__1, i__2;    doublecomplex z__1;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    static integer j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static logical upper;    extern /* Subroutine */ int ztpmv_(char *, char *, char *, integer *, 	    doublecomplex *, doublecomplex *, integer *);    static integer jc, jj;    extern /* Subroutine */ int xerbla_(char *, integer *);    static integer jclast;    static logical nounit;    static doublecomplex ajj;//.........这里部分代码省略.........

//.........这里部分代码省略.........		    iptr = istart + nsupc;		    for (i = 0; i < nrow; i++) {			irow = L_SUB(iptr);			z_sub(&rhs_work[irow], &rhs_work[irow], &work[i]);			work[i].r = 0.;	                work[i].i = 0.;			iptr++;		    }		}#endif		    	    } /* else ... */	} /* for L-solve */#ifdef DEBUG  	printf("After L-solve: y=/n");	zprint_soln(n, nrhs, Bmat);#endif	/*	 * Back solve Ux=y.	 */	for (k = Lstore->nsuper; k >= 0; k--) {	    fsupc = L_FST_SUPC(k);	    istart = L_SUB_START(fsupc);	    nsupr = L_SUB_START(fsupc+1) - istart;	    nsupc = L_FST_SUPC(k+1) - fsupc;	    luptr = L_NZ_START(fsupc);	    solve_ops += 4 * nsupc * (nsupc + 1) * nrhs;	    if ( nsupc == 1 ) {		rhs_work = &Bmat[0];		for (j = 0; j < nrhs; j++) {		    z_div(&rhs_work[fsupc], &rhs_work[fsupc], &Lval[luptr]);		    rhs_work += ldb;		}	    } else {#ifdef USE_VENDOR_BLAS#ifdef _CRAY		ftcs1 = _cptofcd("L", strlen("L"));		ftcs2 = _cptofcd("U", strlen("U"));		ftcs3 = _cptofcd("N", strlen("N"));		CTRSM( ftcs1, ftcs2, ftcs3, ftcs3, &nsupc, &nrhs, &alpha,		       &Lval[luptr], &nsupr, &Bmat[fsupc], &ldb);#else		ztrsm_("L", "U", "N", "N", &nsupc, &nrhs, &alpha,		       &Lval[luptr], &nsupr, &Bmat[fsupc], &ldb);#endif#else				for (j = 0; j < nrhs; j++)		    zusolve ( nsupr, nsupc, &Lval[luptr], &Bmat[fsupc+j*ldb] );#endif			    }	    for (j = 0; j < nrhs; ++j) {		rhs_work = &Bmat[j*ldb];		for (jcol = fsupc; jcol < fsupc + nsupc; jcol++) {		    solve_ops += 8*(U_NZ_START(jcol+1) - U_NZ_START(jcol));		    for (i = U_NZ_START(jcol); i < U_NZ_START(jcol+1); i++ ){			irow = U_SUB(i);			zz_mult(&temp_comp, &rhs_work[jcol], &Uval[i]);			z_sub(&rhs_work[irow], &rhs_work[irow], &temp_comp);		    }		}	    }	    

/* Double Complex */ VOID zlatm2_(doublecomplex * ret_val, integer *m, 	integer *n, integer *i, integer *j, integer *kl, integer *ku, integer 	*idist, integer *iseed, doublecomplex *d, integer *igrade, 	doublecomplex *dl, doublecomplex *dr, integer *ipvtng, integer *iwork,	 doublereal *sparse){    /* System generated locals */    integer 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 */    static integer isub, jsub;    static doublecomplex ctemp;    extern doublereal dlaran_(integer *);    extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 	    integer *);/*  -- LAPACK auxiliary test routine (version 2.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          February 29, 1992       Purpose       =======          ZLATM2 returns the (I,J) entry of a random matrix of dimension          (M, N) described by the other paramters. It is called by the          ZLATMR routine in order to build random test matrices. No error          checking on parameters is done, because this routine is called in          a tight loop by ZLATMR which has already checked the parameters.          Use of ZLATM2 differs from CLATM3 in the order in which the random          number generator is called to fill in random matrix entries.          With ZLATM2, the generator is called to fill in the pivoted matrix          columnwise. With ZLATM3, the generator is called to fill in the          matrix columnwise, after which it is pivoted. Thus, ZLATM3 can          be used to construct random matrices which differ only in their          order of rows and/or columns. ZLATM2 is used to construct band          matrices while avoiding calling the random number generator for          entries outside the band (and therefore generating random numbers          The matrix whose (I,J) entry is returned is constructed as          follows (this routine only computes one entry):            If I is outside (1..M) or J is outside (1..N), return zero               (this is convenient for generating matrices in band format).            Generate a matrix A with random entries of distribution IDIST.            Set the diagonal to D.            Grade the matrix, if desired, from the left (by DL) and/or               from the right (by DR or DL) as specified by IGRADE.            Permute, if desired, the rows and/or columns as specified by               IPVTNG and IWORK.            Band the matrix to have lower bandwidth KL and upper               bandwidth KU.            Set random entries to zero as specified by SPARSE.       Arguments       =========       M      - INTEGER                Number of rows of matrix. Not modified.       N      - INTEGER                Number of columns of matrix. Not modified.       I      - INTEGER                Row of entry to be returned. Not modified.       J      - INTEGER                Column of entry to be returned. Not modified.       KL     - INTEGER                Lower bandwidth. Not modified.       KU     - INTEGER                Upper bandwidth. Not modified.       IDIST  - INTEGER                On entry, IDIST specifies the type of distribution to be                used to generate a random matrix .   //.........这里部分代码省略.........

/* Subroutine */ int ztrsv_(char *uplo, char *trans, char *diag, integer *n, 	doublecomplex *a, integer *lda, doublecomplex *x, integer *incx){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(	    doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, ix, jx, kx, info;    doublecomplex temp;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    logical noconj, nounit;/*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTRSV  solves one of the systems of equations *//*     A*x = b,   or   A'*x = b,   or   conjg( A' )*x = b, *//*  where b and x are n element vectors and A is an n by n unit, or *//*  non-unit, upper or lower triangular matrix. *//*  No test for singularity or near-singularity is included in this *//*  routine. Such tests must be performed before calling this routine. *//*  Arguments *//*  ========== *//*  UPLO   - CHARACTER*1. *//*           On entry, UPLO specifies whether the matrix is an upper or *//*           lower triangular matrix as follows: *//*              UPLO = 'U' or 'u'   A is an upper triangular matrix. *//*              UPLO = 'L' or 'l'   A is a lower triangular matrix. *//*           Unchanged on exit. *//*  TRANS  - CHARACTER*1. *//*           On entry, TRANS specifies the equations to be solved as *//*           follows: *//*              TRANS = 'N' or 'n'   A*x = b. *//*              TRANS = 'T' or 't'   A'*x = b. *//*              TRANS = 'C' or 'c'   conjg( A' )*x = b. *//*           Unchanged on exit. *//*  DIAG   - CHARACTER*1. *//*           On entry, DIAG specifies whether or not A is unit *//*           triangular as follows: *//*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. *//*              DIAG = 'N' or 'n'   A is not assumed to be unit *//*                                  triangular. *//*           Unchanged on exit. *//*  N      - INTEGER. *//*           On entry, N specifies the order of the matrix A. *//*           N must be at least zero. *//*           Unchanged on exit. *//*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). *//*           Before entry with  UPLO = 'U' or 'u', the leading n by n *//*           upper triangular part of the array A must contain the upper *//*           triangular matrix and the strictly lower triangular part of *//*           A is not referenced. *//*           Before entry with UPLO = 'L' or 'l', the leading n by n *//*           lower triangular part of the array A must contain the lower *//*           triangular matrix and the strictly upper triangular part of *//*           A is not referenced. *//*           Note that when  DIAG = 'U' or 'u', the diagonal elements of *//*           A are not referenced either, but are assumed to be unity. *//*           Unchanged on exit. *//*  LDA    - INTEGER. *//*           On entry, LDA specifies the first dimension of A as declared *//*           in the calling (sub) program. LDA must be at least *//*           max( 1, n ). *//*           Unchanged on exit. *//*  X      - COMPLEX*16       array of dimension at least *//*           ( 1 + ( n - 1 )*abs( INCX ) ). *//*           Before entry, the incremented array X must contain the n *//*           element right-hand side vector b. On exit, X is overwritten *///.........这里部分代码省略.........