C++ c_div函数代码示例

void small_conducting_Mie(double x,struct c_complex m,double*mu,
long nangles,struct c_complex*s1,
struct c_complex*s2,double*qext,double*qsca,
double*qback,double*g)

/*:31*/
#line 448 "./mie.w"


{
struct c_complex ahat1,ahat2,bhat1,bhat2;
struct c_complex ss1;
double x2,x3,x4,muj,angle;
long j;

if((s1==NULL)||(s2==NULL))nangles= 0;

m.re+= 0.0;
x2= x*x;
x3= x2*x;
x4= x2*x2;

ahat1= c_div(c_set(0.0,2.0/3.0*(1.0-0.2*x2)),c_set(1.0-0.5*x2,2.0/3.0*x3));
bhat1= c_div(c_set(0.0,(x2-10.0)/30.0),c_set(1+0.5*x2,-x3/3.0));
ahat2= c_set(0.0,x2/30.);
bhat2= c_set(0.0,-x2/45.);

*qsca= 6.0*x4*(c_norm(ahat1)+c_norm(bhat1)+
(5.0/3.0)*(c_norm(ahat2)+c_norm(bhat2)));
*qext= *qsca;
*g= 6.0*x4*(ahat1.im*(ahat2.im+bhat1.im)+
bhat2.im*(5.0/9.0*ahat2.im+bhat1.im)+
ahat1.re*bhat1.re)/(*qsca);

ss1.re= 1.5*x2*(ahat1.re-bhat1.re);
ss1.im= 1.5*x2*(ahat1.im-bhat1.im-(5.0/3.0)*(ahat2.im+bhat2.im));
*qback= 4*c_norm(ss1);

x3*= 1.5;
ahat1.re*= x3;
ahat1.im*= x3;
bhat1.re*= x3;
bhat1.im*= x3;
ahat2.im*= x3*(5.0/3.0);
bhat2.im*= x3*(5.0/3.0);
for(j= 0;j<nangles;j++){
muj= mu[j];
angle= 2.0*muj*muj-1.0;
s1[j].re= ahat1.re+(bhat1.re)*muj;
s1[j].im= ahat1.im+(bhat1.im+ahat2.im)*muj+bhat2.im*angle;;
s2[j].re= bhat1.re+(ahat1.re)*muj;
s2[j].im= bhat1.im+(ahat1.im+bhat2.im)*muj+ahat2.im*angle;
}
}

struct c_complex Lentz_Dn(struct c_complex z,long n)

/*:10*/
#line 126 "./mie.w"

{
struct c_complex alpha_j1,alpha_j2,zinv,aj;
struct c_complex alpha,result,ratio,runratio;

/*12:*/
#line 156 "./mie.w"


zinv= c_sdiv(2.0,z);
alpha= c_smul(n+0.5,zinv);
aj= c_smul(-n-1.5,zinv);
alpha_j1= c_add(aj,c_inv(alpha));
alpha_j2= aj;
ratio= c_div(alpha_j1,alpha_j2);
runratio= c_mul(alpha,ratio);


/*:12*/
#line 131 "./mie.w"


do
/*13:*/
#line 179 "./mie.w"

{
aj.re= zinv.re-aj.re;
aj.im= zinv.im-aj.im;
alpha_j1= c_add(c_inv(alpha_j1),aj);
alpha_j2= c_add(c_inv(alpha_j2),aj);
ratio= c_div(alpha_j1,alpha_j2);
zinv.re*= -1;
zinv.im*= -1;
runratio= c_mul(ratio,runratio);
}

/*:13*/
#line 134 "./mie.w"


while(fabs(c_abs(ratio)-1.0)> 1e-12);

result= c_add(c_sdiv((double)-n,z),runratio);
return result;
}

static struct polynom *pol_div(struct polynom *p,struct polynom *q,struct polynom *p1,struct polynom *p2)
        {
        int i,j;
        struct complex z,z1;

		z.x=0; z.y=0; // RS 7.2.2013
        p->n=p1->n-p2->n;
        if (p->n<0)
            {
            sur_print("\nDegree of dividend < Degree of divisor");
            WAIT; return(p);
            }
        for (i=0; i<=p->n; ++i)
            p->a[i].x=p->a[i].y=0.0;
        for (i=0; i<=p1->n; ++i)
            { q->a[i].x=p1->a[i].x; q->a[i].y=p1->a[i].y; }

        for (i=p1->n; i>=p2->n; --i)
            {
            if (c_zero(&(q->a[i]))) continue;
            c_div(&z,&(q->a[i]),&(p2->a[p2->n]));
            p->a[i-p2->n].x=z.x; p->a[i-p2->n].y=z.y;
            q->a[i].x=q->a[i].y=0.0;
            for (j=i-1; j>=i-p2->n; --j)
                c_sub(&(q->a[j]),&(q->a[j]),
                      c_mult(&z1,&z,&(p2->a[p2->n-i+j])));
            }
        i=p2->n-1;
        while (c_zero(&(q->a[i])) && i>0) --i;
        q->n=i;
        return(p);
        }

/**
 * \brief Expand pole product
 * \param c         resulting filter coefficients
 * \param poles     pole locations
 * \param K         number of poles
 * \ingroup vyv_gaussian
 *
 * This routine expands the product to obtain the filter coefficients:
 * \f[ \prod_{k=0}^{K-1}\frac{\mathrm{poles}[k]-1}{\mathrm{poles}[k]-z^{-1}}
 = \frac{c[0]}{1+\sum_{k=1}^K c[k] z^{-k}}. \f]
 */
static void expand_pole_product(double *c, const complex4c *poles, int K)
{
	complex4c denom[VYV_MAX_K + 1];
	int k, j;

	assert(K <= VYV_MAX_K);
	denom[0] = poles[0];
	denom[1] = make_complex(-1, 0);

	for (k = 1; k < K; ++k)
	{
		denom[k + 1] = c_neg(denom[k]);

		for (j = k; j > 0; --j)
			denom[j] = c_sub(c_mul(denom[j], poles[k]), denom[j - 1]);

		denom[0] = c_mul(denom[0], poles[k]);
	}

	for (k = 1; k <= K; ++k)
		c[k] = c_div(denom[k], denom[0]).real;

	for (c[0] = 1, k = 1; k <= K; ++k)
		c[0] += c[k];

	return;
}

void dd (header *hd)
{	header *st=hd,*hd1,*result;
	int c1,c2,i,j,r;
	double *m1,*m2,*mr;
	complex *mc1,*mc2,*mcr,hc1,hc2;
	interval *mi1,*mi2,*mir,hi1,hi2;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	if (c1!=c2) wrong_arg();
	if (iscomplex(hd)) /* complex values */
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c1,""); if (error) return;
		mcr=(complex *)matrixof(result);
		memmove((char *)mcr,(char *)mc2,c1*sizeof(complex));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	if (mc1[j][0]==mc1[j-i][0] &&
					mc1[j][1]==mc1[j-i][1]) wrong_arg();
				c_sub(mcr[j],mcr[j-1],hc1);
				c_sub(mc1[j],mc1[j-i],hc2);
				c_div(hc1,hc2,mcr[j]);
			}
		}
	}
	else if (isinterval(hd)) /* complex values */
	{	mi1=(complex *)m1; mi2=(complex *)m2;
		result=new_imatrix(1,c1,""); if (error) return;
		mir=(interval *)matrixof(result);
		memmove((char *)mir,(char *)mi2,c1*sizeof(interval));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	i_sub(mir[j],mir[j-1],hi1);
				if (hi1[0]<=0 && hi1[1]>=0)
				{	output("Interval points coincide\n");
					error=1; return;
				}
				i_sub(mi1[j],mi1[j-i],hi2);
				i_div(hi1,hi2,mir[j]);
			}
		}
	}
	else if (isreal(hd))
	{	result=new_matrix(1,c1,""); if (error) return;
		mr=matrixof(result);
		memmove((char *)mr,(char *)m2,c1*sizeof(double));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	if (m1[j]==m1[j-i]) wrong_arg();
				mr[j]=(mr[j]-mr[j-1])/(m1[j]-m1[j-i]);
			}
		}	
	}
	else wrong_arg();
	moveresult(st,result);
}

/**
 * \brief Compute the variance of the impulse response
 * \param poles0    unscaled pole locations
 * \param q         rescaling parameter
 * \param K         number of poles
 * \return variance achieved by poles = poles0^(1/q)
 * \ingroup vyv_gaussian
 */
static double variance(const complex4c *poles0, int K, double q)
{
	complex4c sum = { 0, 0 };
	int k;

	for (k = 0; k < K; ++k)
	{
		complex4c z = c_real_pow(poles0[k], 1 / q), denom = z;
		denom.real -= 1;
		/* Compute sum += z / (z - 1)^2. */
		sum = c_add(sum, c_div(z, c_mul(denom, denom)));
	}

	return 2 * sum.real;
}

/**
 * \brief Derivative of variance with respect to q
 * \param poles0    unscaled pole locations
 * \param q         rescaling parameter
 * \param K         number of poles
 * \return derivative of variance with respect to q
 * \ingroup vyv_gaussian
 *
 * This function is used by compute_q() in solving for q.
 */
static double dq_variance(const complex4c *poles0, int K, double q)
{
	complex4c sum = { 0, 0 };
	int k;

	for (k = 0; k < K; ++k)
	{
		complex4c z = c_real_pow(poles0[k], 1 / q), w = z, denom = z;
		w.real += 1;
		denom.real -= 1;
		/* Compute sum += z log(z) (z + 1) / (z - 1)^3 */
		sum = c_add(sum, c_div(c_mul(c_mul(z, c_log(z)), w),
			c_real_pow(denom, 3)));
	}

	return (2 / q) * sum.real;
}

//.........这里部分代码省略.........
    int          *lsub, *xlsub;
    complex       *lusup;
    int          *xlusup;
    flops_t      *ops = stat->ops;

    /* Initialize pointers */
    lsub       = Glu->lsub;
    xlsub      = Glu->xlsub;
    lusup      = Glu->lusup;
    xlusup     = Glu->xlusup;
    fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
    nsupc      = jcol - fsupc;	        /* excluding jcol; nsupc >= 0 */
    lptr       = xlsub[fsupc];
    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 DEBUG
if ( jcol == MIN_COL ) {
    printf(_T("Before cdiv: col %d\n"), jcol);
    for (k = nsupc; k < nsupr; k++) 
	printf(_T("  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 = c_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 ) {
	*pivrow = lsub_ptr[pivptr];
	perm_r[*pivrow] = jcol;
	*usepr = 0;
	return (jcol+1);
    }

    thresh = u * pivmax;
    
    /* Choose appropriate pivotal element by our policy. */
    if ( *usepr ) {
        rtemp = c_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 = c_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);

    c_div(&temp, &one, &lu_col_ptr[nsupc]);
    for (k = nsupc+1; k < nsupr; k++) 
	cc_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);

    return 0;
}

/* Subroutine */ int ctrti2_(char *uplo, char *diag, integer *n, complex *a, 
	integer *lda, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;
    complex q__1;

    /* Local variables */
    integer j;
    complex ajj;
    logical upper;
    logical nounit;

/*  -- LAPACK routine (version 3.2) -- */
/*     November 2006 */

/*  Purpose */
/*  ======= */

/*  CTRTI2 computes the inverse of a complex upper or lower triangular */
/*  matrix. */

/*  This is the Level 2 BLAS version of the algorithm. */

/*  Arguments */
/*  ========= */

/*  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. */

/*  A       (input/output) COMPLEX array, dimension (LDA,N) */
/*          On entry, the triangular matrix A.  If UPLO = 'U', the */
/*          leading n by n upper triangular part of the array A contains */
/*          the upper triangular matrix, and the strictly lower */
/*          triangular part of A is not referenced.  If UPLO = 'L', the */
/*          leading n by n lower triangular part of the array A contains */
/*          the lower triangular matrix, and the strictly upper */
/*          triangular part of A is not referenced.  If DIAG = 'U', the */
/*          diagonal elements of A are also not referenced and are */
/*          assumed to be 1. */

/*          On exit, the (triangular) inverse of the original matrix, in */
/*          the same storage format. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,N). */

/*  INFO    (output) INTEGER */
/*          = 0: successful exit */
/*          < 0: if INFO = -k, the k-th argument had an illegal value */

/*  ===================================================================== */

/*     Test the input parameters. */

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    nounit = lsame_(diag, "N");
    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (! nounit && ! lsame_(diag, "U")) {
	*info = -2;
    } else if (*n < 0) {
	*info = -3;
    } else if (*lda < max(1,*n)) {
	*info = -5;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("CTRTI2", &i__1);
	return 0;
    }

    if (upper) {

/*        Compute inverse of upper triangular matrix. */

	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    if (nounit) {
		i__2 = j + j * a_dim1;
		c_div(&q__1, &c_b1, &a[j + j * a_dim1]);
		a[i__2].r = q__1.r, a[i__2].i = q__1.i;
		i__2 = j + j * a_dim1;
//.........这里部分代码省略.........

/* Subroutine */ int claic1_(integer *job, integer *j, complex *x, real *sest,
	 complex *w, complex *gamma, real *sestpr, complex *s, complex *c__)
{
/*  -- 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   
    =======   

    CLAIC1 applies one step of incremental condition estimation in   
    its simplest version:   

    Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j   
    lower triangular matrix L, such that   
             twonorm(L*x) = sest   
    Then CLAIC1 computes sestpr, s, c such that   
    the vector   
                    [ s*x ]   
             xhat = [  c  ]   
    is an approximate singular vector of   
                    [ L     0  ]   
             Lhat = [ w' gamma ]   
    in the sense that   
             twonorm(Lhat*xhat) = sestpr.   

    Depending on JOB, an estimate for the largest or smallest singular   
    value is computed.   

    Note that [s c]' and sestpr**2 is an eigenpair of the system   

        diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]   
                                              [ conjg(gamma) ]   

    where  alpha =  conjg(x)'*w.   

    Arguments   
    =========   

    JOB     (input) INTEGER   
            = 1: an estimate for the largest singular value is computed.   
            = 2: an estimate for the smallest singular value is computed.   

    J       (input) INTEGER   
            Length of X and W   

    X       (input) COMPLEX array, dimension (J)   
            The j-vector x.   

    SEST    (input) REAL   
            Estimated singular value of j by j matrix L   

    W       (input) COMPLEX array, dimension (J)   
            The j-vector w.   

    GAMMA   (input) COMPLEX   
            The diagonal element gamma.   

    SESTPR  (output) REAL   
            Estimated singular value of (j+1) by (j+1) matrix Lhat.   

    S       (output) COMPLEX   
            Sine needed in forming xhat.   

    C       (output) COMPLEX   
            Cosine needed in forming xhat.   

    =====================================================================   


       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    
    /* System generated locals */
    real r__1, r__2;
    complex q__1, q__2, q__3, q__4, q__5, q__6;
    /* Builtin functions */
    double c_abs(complex *);
    void r_cnjg(complex *, complex *), c_sqrt(complex *, complex *);
    double sqrt(doublereal);
    void c_div(complex *, complex *, complex *);
    /* Local variables */
    static complex sine;
    static real test, zeta1, zeta2, b, t;
    static complex alpha;
    extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 
	    *, complex *, integer *);
    static real norma, s1, s2, absgam, absalp;
    extern doublereal slamch_(char *);
    static complex cosine;
    static real absest, scl, eps, tmp;


    --w;
    --x;

    /* Function Body */
//.........这里部分代码省略.........

/* Subroutine */ int cgttrf_(integer *n, complex *dl, complex *d__, complex *
                             du, complex *du2, integer *ipiv, integer *info)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    real r__1, r__2, r__3, r__4;
    complex q__1, q__2;

    /* Builtin functions */
    double r_imag(complex *);
    void c_div(complex *, complex *, complex *);

    /* Local variables */
    integer i__;
    complex fact, temp;
    extern /* Subroutine */ int xerbla_(char *, integer *);


    /*  -- LAPACK routine (version 3.2) -- */
    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
    /*     November 2006 */

    /*     .. Scalar Arguments .. */
    /*     .. */
    /*     .. Array Arguments .. */
    /*     .. */

    /*  Purpose */
    /*  ======= */

    /*  CGTTRF 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 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 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 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 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 .. */

    /* Parameter adjustments */
//.........这里部分代码省略.........

/*<       subroutine csvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) >*/
/* Subroutine */ int csvdc_(complex *x, integer *ldx, integer *n, integer *p, 
        complex *s, complex *e, complex *u, integer *ldu, complex *v, integer 
        *ldv, complex *work, integer *job, integer *info)
{
    /* System generated locals */
    integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 
            i__3, i__4;
    real r__1, r__2, r__3, r__4;
    complex q__1, q__2, q__3;

    /* Builtin functions */
    double r_imag(complex *), c_abs(complex *);
    void c_div(complex *, complex *, complex *), r_cnjg(complex *, complex *);
    double sqrt(doublereal);

    /* Local variables */
    real b, c__, f, g;
    integer i__, j, k, l=0, m;
    complex r__, t;
    real t1, el;
    integer kk;
    real cs;
    integer ll, mm, ls=0;
    real sl;
    integer lu;
    real sm, sn;
    integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
    real emm1, smm1;
    integer kase, jobu, iter;
    real test;
    integer nctp1, nrtp1;
    extern /* Subroutine */ int cscal_(integer *, complex *, complex *, 
            integer *);
    real scale;
    extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 
            *, complex *, integer *);
    real shift;
    extern /* Subroutine */ int cswap_(integer *, complex *, integer *, 
            complex *, integer *);
    integer maxit;
    extern /* Subroutine */ int caxpy_(integer *, complex *, complex *, 
            integer *, complex *, integer *), csrot_(integer *, complex *, 
            integer *, complex *, integer *, real *, real *);
    logical wantu, wantv;
    extern /* Subroutine */ int srotg_(real *, real *, real *, real *);
    real ztest;
    extern doublereal scnrm2_(integer *, complex *, integer *);

/*<       integer ldx,n,p,ldu,ldv,job,info >*/
/*<       complex x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) >*/


/*     csvdc is a subroutine to reduce a complex nxp matrix x by */
/*     unitary transformations u and v to diagonal form.  the */
/*     diagonal elements s(i) are the singular values of x.  the */
/*     columns of u are the corresponding left singular vectors, */
/*     and the columns of v the right singular vectors. */

/*     on entry */

/*         x         complex(ldx,p), where ldx.ge.n. */
/*                   x contains the matrix whose singular value */
/*                   decomposition is to be computed.  x is */
/*                   destroyed by csvdc. */

/*         ldx       integer. */
/*                   ldx is the leading dimension of the array x. */

/*         n         integer. */
/*                   n is the number of rows of the matrix x. */

/*         p         integer. */
/*                   p is the number of columns of the matrix x. */

/*         ldu       integer. */
/*                   ldu is the leading dimension of the array u */
/*                   (see below). */

/*         ldv       integer. */
/*                   ldv is the leading dimension of the array v */
/*                   (see below). */

/*         work      complex(n). */
/*                   work is a scratch array. */

/*         job       integer. */
/*                   job controls the computation of the singular */
/*                   vectors.  it has the decimal expansion ab */
/*                   with the following meaning */

/*                        a.eq.0    do not compute the left singular */
/*                                  vectors. */
/*                        a.eq.1    return the n left singular vectors */
/*                                  in u. */
/*                        a.ge.2    returns the first min(n,p) */
/*                                  left singular vectors in u. */
/*                        b.eq.0    do not compute the right singular */
/*                                  vectors. */
/*                        b.eq.1    return the right singular vectors */
//.........这里部分代码省略.........

/** CHETRF_ROOK_REC2 computes a partial factorization of a complex Hermitian indefinite matrix using the boun ded Bunch-Kaufman ("rook") diagonal pivoting method
 *
 * This routine is a minor modification of LAPACK's clahef_rook.
 * It serves as an unblocked kernel in the recursive algorithms.
 * The blocked BLAS Level 3 updates were removed and moved to the
 * recursive algorithm.
 * */
/* Subroutine */ void RELAPACK_chetrf_rook_rec2(char *uplo, int *n,
	int *nb, int *kb, complex *a, int *lda, int *ipiv,
	complex *w, int *ldw, int *info, ftnlen uplo_len)
{
    /* System generated locals */
    int a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4;
    float r__1, r__2;
    complex q__1, q__2, q__3, q__4, q__5;

    /* Builtin functions */
    double sqrt(double), r_imag(complex *);
    void r_cnjg(complex *, complex *), c_div(complex *, complex *, complex *);

    /* Local variables */
    static int j, k, p;
    static float t, r1;
    static complex d11, d21, d22;
    static int ii, jj, kk, kp, kw, jp1, jp2, kkw;
    static logical done;
    static int imax, jmax;
    static float alpha;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    extern /* Subroutine */ int cgemv_(char *, int *, int *, complex *
	    , complex *, int *, complex *, int *, complex *, complex *
	    , int *, ftnlen);
    static float sfmin;
    extern /* Subroutine */ int ccopy_(int *, complex *, int *,
	    complex *, int *);
    static int itemp;
    extern /* Subroutine */ int cswap_(int *, complex *, int *,
	    complex *, int *);
    static int kstep;
    static float stemp, absakk;
    extern /* Subroutine */ int clacgv_(int *, complex *, int *);
    extern int icamax_(int *, complex *, int *);
    extern double slamch_(char *, ftnlen);
    extern /* Subroutine */ int csscal_(int *, float *, complex *, int
	    *);
    static float colmax, rowmax;

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --ipiv;
    w_dim1 = *ldw;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    *info = 0;
    alpha = (sqrt(17.f) + 1.f) / 8.f;
    sfmin = slamch_("S", (ftnlen)1);
    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
	k = *n;
L10:
	kw = *nb + k - *n;
	if ((k <= *n - *nb + 1 && *nb < *n) || k < 1) {
	    goto L30;
	}
	kstep = 1;
	p = k;
	if (k > 1) {
	    i__1 = k - 1;
	    ccopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &w[kw * w_dim1 + 1], &
		    c__1);
	}
	i__1 = k + kw * w_dim1;
	i__2 = k + k * a_dim1;
	r__1 = a[i__2].r;
	w[i__1].r = r__1, w[i__1].i = 0.f;
	if (k < *n) {
	    i__1 = *n - k;
	    q__1.r = -1.f, q__1.i = -0.f;
	    cgemv_("No transpose", &k, &i__1, &q__1, &a[(k + 1) * a_dim1 + 1],
		     lda, &w[k + (kw + 1) * w_dim1], ldw, &c_b1, &w[kw *
		    w_dim1 + 1], &c__1, (ftnlen)12);
	    i__1 = k + kw * w_dim1;
	    i__2 = k + kw * w_dim1;
	    r__1 = w[i__2].r;
	    w[i__1].r = r__1, w[i__1].i = 0.f;
	}
	i__1 = k + kw * w_dim1;
	absakk = (r__1 = w[i__1].r, dabs(r__1));
	if (k > 1) {
	    i__1 = k - 1;
	    imax = icamax_(&i__1, &w[kw * w_dim1 + 1], &c__1);
	    i__1 = imax + kw * w_dim1;
	    colmax = (r__1 = w[i__1].r, dabs(r__1)) + (r__2 = r_imag(&w[imax
		    + kw * w_dim1]), dabs(r__2));
	} else {
	    colmax = 0.f;
	}
//.........这里部分代码省略.........

static int polroot(
struct polynom *p,
struct complex *pz,  /* pointer to root */
struct complex *pz0  /* pointer to initial value */
)
        {
        int i;
        int n_iter;
        struct polynom d;
        struct complex v,v0,vd,delta,zmin;
        double y,ymin;
//        extern struct polynom *pol_der();
//        extern struct complex *pol_value();
		delta.x=0; delta.y=0; // RS 7.2.2013
//printf("\np->n=%d",p->n); getch();
        if (p->n==1)
            {
            c_div(pz,&(p->a[0]),&(p->a[1]));
            pz->x=-(pz->x); pz->y=-(pz->y);
            return(1);
            }
//printf("\nr2");
        pol_value(p,pz0,&v0);
//printf("\nr3");
        if (c_zero(&v0))
            {
            pz->x=pz0->x;
            pz->y=pz0->y;
            return(1);
            }
//printf("\nr4");
        zmin.x=pz0->x; zmin.y=pz0->y;
        ymin=v0.x*v0.x+v0.y*v0.y;
        pol_der(&d,p);
//printf("\nr5");
        n_iter=0;
        while (1)
            {
            pol_value(&d,pz0,&vd);

            c_div(&delta,&v0,&vd);
            c_sub(pz,pz0,&delta);
            pol_value(p,pz,&v);
            v0.x=v.x; v0.y=v.y;
            pz0->x=pz->x; pz0->y=pz->y;

            ++n_iter;
            PR_UP;
            sprintf(sbuf,"\npolroot: N=%d  Re=%e Im=%e",n_iter,pz->x,pz->y); sur_print(sbuf); // RS CHA Rprintf

/* Rprintf("zero: %d\n",c_zero(pz)); getch();  */
  /*        if (c_zero(&pz)) break;     */
            if (c_zero(pz)) break;
            y=v.x*v.x+v.y*v.y;
            if (y<ymin)
                { ymin=y; zmin.x=pz->x; zmin.y=pz->y; }

            if (fabs(delta.x)<roots_eps && fabs(delta.y)<roots_eps) break;

            if (n_iter>roots_max_iter)
                {
                pz->x=zmin.x; pz->y=zmin.y;
                break;
                }
            if (sur_kbhit())
                {
                pz->x=zmin.x; pz->y=zmin.y;
                i=sur_getch(); if (i=='.') break;
                }
            }
        return (n_iter);
        }

/* Subroutine */ int clarge_(integer *n, complex *a, integer *lda, integer *
                             iseed, complex *work, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;
    real r__1;
    complex q__1;

    /* Builtin functions */
    double c_abs(complex *);
    void c_div(complex *, complex *, complex *);

    /* Local variables */
    static integer i__;
    extern /* Subroutine */ int cgerc_(integer *, integer *, complex *,
                                       complex *, integer *, complex *, integer *, complex *, integer *),
                                               cscal_(integer *, complex *, complex *, integer *), cgemv_(char *
                                                       , integer *, integer *, complex *, complex *, integer *, complex *
                                                       , integer *, complex *, complex *, integer *);
    extern doublereal scnrm2_(integer *, complex *, integer *);
    static complex wa, wb;
    static real wn;
    extern /* Subroutine */ int xerbla_(char *, integer *), clarnv_(
        integer *, integer *, integer *, complex *);
    static complex tau;


#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)]


    /*  -- LAPACK auxiliary test 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
        =======

        CLARGE pre- and post-multiplies a complex general n by n matrix A
        with a random unitary matrix: A = U*D*U'.

        Arguments
        =========

        N       (input) INTEGER
                The order of the matrix A.  N >= 0.

        A       (input/output) COMPLEX array, dimension (LDA,N)
                On entry, the original n by n matrix A.
                On exit, A is overwritten by U*A*U' for some random
                unitary matrix U.

        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 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 */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    --iseed;
    --work;

    /* Function Body */
    *info = 0;
    if (*n < 0) {
        *info = -1;
    } else if (*lda < max(1,*n)) {
        *info = -3;
    }
    if (*info < 0) {
        i__1 = -(*info);
        xerbla_("CLARGE", &i__1);
        return 0;
    }

    /*     pre- and post-multiply A by random unitary matrix */

    for (i__ = *n; i__ >= 1; --i__) {

        /*        generate random reflection */

//.........这里部分代码省略.........