float clangs(char *norm, SuperMatrix *A)
{
    
    /* Local variables */
    NCformat *Astore;
    complex   *Aval;
    int      i, j, irow;
    float   value, sum;
    float   *rwork;

    Astore = A->Store;
    Aval   = Astore->nzval;
    
    if ( SUPERLU_MIN(A->nrow, A->ncol) == 0) {
	value = 0.;
	
    } else if (strncmp(norm, "M", 1)==0) {
	/* Find max(abs(A(i,j))). */
	value = 0.;
	for (j = 0; j < A->ncol; ++j)
	    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++)
		value = SUPERLU_MAX( value, c_abs( &Aval[i]) );
	
    } else if (strncmp(norm, "O", 1)==0 || *(unsigned char *)norm == '1') {
	/* Find norm1(A). */
	value = 0.;
	for (j = 0; j < A->ncol; ++j) {
	    sum = 0.;
	    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++) 
		sum += c_abs( &Aval[i] );
	    value = SUPERLU_MAX(value,sum);
	}
	
    } else if (strncmp(norm, "I", 1)==0) {
	/* Find normI(A). */
	if ( !(rwork = (float *) SUPERLU_MALLOC(A->nrow * sizeof(float))) )
	    ABORT("SUPERLU_MALLOC fails for rwork.");
	for (i = 0; i < A->nrow; ++i) rwork[i] = 0.;
	for (j = 0; j < A->ncol; ++j)
	    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++) {
		irow = Astore->rowind[i];
		rwork[irow] += c_abs( &Aval[i] );
	    }
	value = 0.;
	for (i = 0; i < A->nrow; ++i)
	    value = SUPERLU_MAX(value, rwork[i]);
	
	SUPERLU_FREE (rwork);
	
    } else if (strncmp(norm, "F", 1)==0 || strncmp(norm, "E", 1)==0) {
	/* Find normF(A). */
	ABORT("Not implemented.");
    } else
	ABORT("Illegal norm specified.");

    return (value);

} /* clangs */

/*! \brief

<pre>
    Purpose
    =======

    SCSUM1 takes the sum of the absolute values of a complex
    vector and returns a single precision result.

    Based on SCASUM from the Level 1 BLAS.
    The change is to use the 'genuine' absolute value.

    Contributed by Nick Higham for use with CLACON.

    Arguments
    =========

    N       (input) INT
            The number of elements in the vector CX.

    CX      (input) COMPLEX array, dimension (N)
            The vector whose elements will be summed.

    INCX    (input) INT
            The spacing between successive values of CX.  INCX > 0.

    =====================================================================
</pre>
*/
double scsum1_(int *n, complex *cx, int *incx)
{
    /* System generated locals */
    int i__1, i__2;
    float ret_val;
    /* Builtin functions */
    double c_abs(complex *);
    /* Local variables */
    static int i, nincx;
    static float stemp;


#define CX(I) cx[(I)-1]


    ret_val = 0.f;
    stemp = 0.f;
    if (*n <= 0) {
        return ret_val;
    }
    if (*incx == 1) {
        goto L20;
    }

/*     CODE FOR INCREMENT NOT EQUAL TO 1 */

    nincx = *n * *incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i = 1; *incx < 0 ? i >= nincx : i <= nincx; i += *incx) {

/*        NEXT LINE MODIFIED. */

        stemp += c_abs(&CX(i));
/* L10: */
    }
    ret_val = stemp;
    return ret_val;

/*     CODE FOR INCREMENT EQUAL TO 1 */

L20:
    i__2 = *n;
    for (i = 1; i <= *n; ++i) {

/*        NEXT LINE MODIFIED. */

        stemp += c_abs(&CX(i));
/* L30: */
    }
    ret_val = stemp;
    return ret_val;

/*     End of SCSUM1 */

} /* scsum1_ */

long RoundDoubleLong(double val)
{
    long ival = c_abs((long)val);
    if((c_abs(val) - ival) >= 0.5)
        ival++;

	if(val > 0)
	{
		return ival;
	}
	else
	{
		return 0 - ival;
	}
}

struct c_float mult_cc(struct c_float a, struct c_float b)
{
    struct c_float c;
    int sign_c = 1;

    if (a.mantisse < 0)
    {
        sign_c = -1;
        a.mantisse = -a.mantisse;
    }

    if (b.mantisse < 0)
    {
        sign_c = -sign_c;
        b.mantisse = -b.mantisse;
    }

    a.mantisse = a.mantisse >> (INT_MAXE / 2);
    b.mantisse = b.mantisse >> (INT_MAXE / 2);

    c.exponent = a.exponent + b.exponent;
    c.mantisse = a.mantisse * b.mantisse * sign_c;

    if (c_abs(c.mantisse) < INT_MAXX)
    {
        c.exponent--;
        c.mantisse = c.mantisse << 1;
    };

    return c;
}

void MainContentComponent::calc_tuning() {								//calculate tuning

	//initialise fft i/o-buffers:									
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * file_len);
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * file_len);
	p = fftw_plan_dft_1d(file_len, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
	for(long n=0; n<file_len; n++) {									// fill fft input buffer:
		in[n][0] = file_buf[n];		//	real part
		in[n][1] = 0;				//	imag part
	}
	fftw_execute(p);													// execute fft

	file_mX = new float[file_len];										//get mag and phase:
	file_pX = new float[file_len];
	for(long n=0; n<file_len/2; n++) {
		file_mX[n] = c_abs(out[n][0], out[n][1]);
		file_pX[n] = c_arg(out[n][0], out[n][1]);
/*
		std::cout<<"magnitude: "<<file_mX[n];
		std::cout<<"\tphase: "<<file_pX[n]<<std::endl;
*/
	}
		fftw_destroy_plan(p);
	fftw_free(in); fftw_free(out); 
	ladder = new tuning(file_mX, file_len/2, key_count, mirror_tuning, peak_tresh, file_fs);	//(spectrum, length of spectrum (N/2 ivm symmetrie),
																								//	amount of keys, mirror intervals round 3/2?, 
																								//		peak detection treshold)  
	

}

//只能接收长度为128的数组
//若不足128补全为0
double* fftMe(double list[],int length) {
    complex theList[128];
    int i;
    for (i = 0; i < length; i++) {
        complex co;
        co.real=list[i];
        co.imag=0;
        theList[i] = co;
    }
    for (i = length; i < 128; i++) {
        complex co;
        co.imag=0;
        co.real=0;
        theList[i] = co;
    }

    // fft
    fft(length,theList);

    double alpha=1.0/(double)length;
    for (i = 0; i < length; i++) {
        complex co;
        co.real = alpha*theList[i].real;
        co.imag = alpha*theList[i].imag;
        theList[i] = co;
    }
    float fftSeries[length];
    double out[length/2];
    c_abs(theList,fftSeries,length);
    for(i=0;i<length/2;i++){
       out[i] = fftSeries[i+1]*2;
       //printf("%.10f\n",out[i]);
    }
    return out;
}

/*! \brief SIGN functions for complex number. Returns z/abs(z) */
complex c_sgn(complex *z)
{
    register float t = c_abs(z);
    register complex retval;

    if (t == 0.0) {
	retval.r = 1.0, retval.i = 0.0;
    } else {
	retval.r = z->r / t, retval.i = z->i / t;
    }

    return retval;
}

/*! \brief Check the inf-norm of the error vector 
 */
void cinf_norm_error(int nrhs, SuperMatrix *X, complex *xtrue)
{
    DNformat *Xstore;
    float err, xnorm;
    complex *Xmat, *soln_work;
    complex temp;
    int i, j;

    Xstore = X->Store;
    Xmat = Xstore->nzval;

    for (j = 0; j < nrhs; j++) {
      soln_work = &Xmat[j*Xstore->lda];
      err = xnorm = 0.0;
      for (i = 0; i < X->nrow; i++) {
        c_sub(&temp, &soln_work[i], &xtrue[i]);
	err = SUPERLU_MAX(err, c_abs(&temp));
	xnorm = SUPERLU_MAX(xnorm, c_abs(&soln_work[i]));
      }
      err = err / xnorm;
      printf("||X - Xtrue||/||X|| = %e\n", err);
    }
}

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

complex c_sqrt(complex x)  /* 平方根 $\sqrt{x}$ */
{
    double r, w;

    r = c_abs(x);
    w = sqrt(r + fabs(x.re));
    if (x.re >= 0) {
        x.re = SQRT05 * w;
        x.im = SQRT05 * x.im / w;
    } else {
        x.re = SQRT05 * fabs(x.im) / w;
        x.im = (x.im >= 0) ? SQRT05 * w : -SQRT05 * w;
    }
    return x;
}

static PyObject *
complex_abs(PyComplexObject *v)
{
	double result;

	PyFPE_START_PROTECT("complex_abs", return 0)
	result = c_abs(v->cval);
	PyFPE_END_PROTECT(result)

	if (errno == ERANGE) {
		PyErr_SetString(PyExc_OverflowError,
				"absolute value too large");
		return NULL;
	}
	return PyFloat_FromDouble(result);
}

/* Subroutine */ int PASTEF77(c,rotg)(singlecomplex *ca, singlecomplex *cb, real *c__, singlecomplex *s)
{
    /* System generated locals */
    real r__1, r__2;
    singlecomplex q__1, q__2, q__3;

    /* Builtin functions */
    double c_abs(singlecomplex *), sqrt(doublereal);
    void bla_r_cnjg(singlecomplex *, singlecomplex *);

    /* Local variables */
    real norm;
    singlecomplex alpha;
    real scale;

    if (c_abs(ca) != 0.f) {
	goto L10;
    }
    *c__ = 0.f;
    s->real = 1.f, s->imag = 0.f;
    ca->real = cb->real, ca->imag = cb->imag;
    goto L20;
L10:
    scale = c_abs(ca) + c_abs(cb);
    q__1.real = ca->real / scale, q__1.imag = ca->imag / scale;
/* Computing 2nd power */
    r__1 = c_abs(&q__1);
    q__2.real = cb->real / scale, q__2.imag = cb->imag / scale;
/* Computing 2nd power */
    r__2 = c_abs(&q__2);
    norm = scale * sqrt(r__1 * r__1 + r__2 * r__2);
    r__1 = c_abs(ca);
    q__1.real = ca->real / r__1, q__1.imag = ca->imag / r__1;
    alpha.real = q__1.real, alpha.imag = q__1.imag;
    *c__ = c_abs(ca) / norm;
    bla_r_cnjg(&q__3, cb);
    q__2.real = alpha.real * q__3.real - alpha.imag * q__3.imag, q__2.imag = alpha.real * q__3.imag + 
	    alpha.imag * q__3.real;
    q__1.real = q__2.real / norm, q__1.imag = q__2.imag / norm;
    s->real = q__1.real, s->imag = q__1.imag;
    q__1.real = norm * alpha.real, q__1.imag = norm * alpha.imag;
    ca->real = q__1.real, ca->imag = q__1.imag;
L20:
    return 0;
} /* crotg_ */

int fft(float *cbd,float *dx2,int *idxcnt)
{

    static integer nbin = 64;
    /*static integer nxx = (6*64)+150;*/

    /* Local variables */
    complex chat[64];
    float magn[64],wk[534];
    integer iwk[534],nbinx=nbin;
    float a, x, dx;
    integer icbd,icnt,iok;

/*       ** perform FFT ananlysis on data */
    for (icbd = 0; icbd < nbin; ++icbd) {
	chat[icbd].r = cbd[icbd], chat[icbd].i = (float)0.;
    }
    iok=fftcc_(chat, &nbinx, iwk, wk);
    if (iok!=0) return -1;
    for (icbd = 0; icbd < nbin; ++icbd) {
	magn[icbd] = c_abs(&chat[icbd]) / (double)nbin;
    }
    dx = 0.0F;
    icnt = 0;
    for (icbd = 1; icbd < (nbin/2)-1; ++icbd) {
	a = magn[icbd-1];
	x = magn[icbd];
	dx += (x + a) / 2.0F;
	if ((magn[icbd] > magn[icbd-1])&&(magn[icbd] > magn[icbd+1])) {
	    ++icnt;
	}
    }
/* added to "normalize" fft for various satellite resolutions */
    *dx2 = dx / magn[0];
    *idxcnt=icnt;

    return 0;
} /* MAIN__ */

//.........这里部分代码省略.........
		for(i=0; i<training_data_n; i++)
		{
		trn_datapair_error_sorted[0][i]=trn_datapair_error[i];
		trn_datapair_error_sorted[1][i]= i+1;
		}

		for(j=1; j<training_data_n; j++)
		{		
		for(i=0; i<training_data_n-j; i++)
		{
		if(trn_datapair_error_sorted[0][i]>trn_datapair_error_sorted[0][i+1])
		{	
		sorting=trn_datapair_error_sorted[0][i+1];
		trn_datapair_error_sorted[0][i+1]=trn_datapair_error_sorted[0][i];
		trn_datapair_error_sorted[0][i]=sorting;
		sortingindex = sorting=trn_datapair_error_sorted[1][i+1];
		trn_datapair_error_sorted[1][i+1]=trn_datapair_error_sorted[1][i];
		trn_datapair_error_sorted[1][i]=sortingindex;
		}
		}
		}

		for(j=0; j<training_data_n; j++)
		{
		fprintf(fppp, "\n");		
		fprintf(fppp, "training data pair sorted number \t %d \n", j+1);
		fprintf(fppp, "training data pair original number \t %d \n", (int)(trn_datapair_error_sorted[1][j]));
		fprintf(fppp, "training data pair sorted error in RMSE is \t %lf \n",trn_datapair_error_sorted[0][j]);
		fprintf(fpppp, "%d \t", (int)(trn_datapair_error_sorted[1][j]));
		complexsum = complex(0.0, 0.0);
		fprintf(fppp,"Normalized layer 3 outputs are as follows \n");
		for(k= In_n*Mf_n + Rule_n; k< In_n*Mf_n + 2*Rule_n; k++)
		{
		fprintf(fppp, "%d.\t %lf + i%lf \t %lf < %lf \n", k, (layer_1_to_4_output[j][k]).real,(layer_1_to_4_output[j][k]).imag, c_abs(layer_1_to_4_output[j][k]), c_phase(layer_1_to_4_output[j][k])*180/PI);
		complexsum = c_add(complexsum, layer_1_to_4_output[j][k]);
		}
		
		
		fprintf(fppp, "Sum of the outputs of layer 3 is \t %lf+i%lf \t %lf<%lf \n", complexsum.real, complexsum.imag, c_abs(complexsum), c_phase(complexsum)*180/PI);
		complexsum = complex(0.0, 0.0);
		fprintf(fppp,"dot producted layer 4 outputs are as follows \n");
		for(k=In_n*Mf_n + 2*Rule_n; k< In_n*Mf_n + 3*Rule_n; k++)
		{
		
		fprintf(fppp, "%d.\t %lf + i%lf \t %lf < %lf \n", k, (layer_1_to_4_output[j][k]).real,(layer_1_to_4_output[j][k]).imag, c_abs(layer_1_to_4_output[j][k]), c_phase(layer_1_to_4_output[j][k])*180/PI);
		complexsum = c_add(complexsum, layer_1_to_4_output[j][k]);
		}

		fprintf(fppp, "sum of the outputs of layer 4 is \t %lf +i%lf \t %lf<%lf \n", complexsum.real, complexsum.imag, c_abs(complexsum), c_phase(complexsum)*180/PI);
		if(j> training_data_n -6 )
		{
		trnnumcheck[(int)(trn_datapair_error_sorted[1][j])]= trnnumcheck[(int)(trn_datapair_error_sorted[1][j])] +1;
		}
		if(j<5 )
		{
		trnnumchecku[(int)(trn_datapair_error_sorted[1][j])]= trnnumchecku[(int)(trn_datapair_error_sorted[1][j])] +1;
		}

		}
		fprintf(fpppp, "\n");
		
		//Find RMSE of testing error
/********************************************************************************
if(checking_data_n != 0) 
		{
			printf("testing 3 \n");	

 int chgeqz_(char *job, char *compq, char *compz, int *n, 
	int *ilo, int *ihi, complex *h__, int *ldh, complex *t, 
	int *ldt, complex *alpha, complex *beta, complex *q, int *ldq, 
	 complex *z__, int *ldz, complex *work, int *lwork, float *
	rwork, int *info)
{
    /* System generated locals */
    int h_dim1, h_offset, q_dim1, q_offset, t_dim1, t_offset, z_dim1, 
	    z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
    float r__1, r__2, r__3, r__4, r__5, r__6;
    complex q__1, q__2, q__3, q__4, q__5, q__6;

    /* Builtin functions */
    double c_abs(complex *);
    void r_cnjg(complex *, complex *);
    double r_imag(complex *);
    void c_div(complex *, complex *, complex *), pow_ci(complex *, complex *, 
	    int *), c_sqrt(complex *, complex *);

    /* Local variables */
    float c__;
    int j;
    complex s, t1;
    int jc, in;
    complex u12;
    int jr;
    complex ad11, ad12, ad21, ad22;
    int jch;
    int ilq, ilz;
    float ulp;
    complex abi22;
    float absb, atol, btol, temp;
    extern  int crot_(int *, complex *, int *, 
	    complex *, int *, float *, complex *);
    float temp2;
    extern  int cscal_(int *, complex *, complex *, 
	    int *);
    extern int lsame_(char *, char *);
    complex ctemp;
    int iiter, ilast, jiter;
    float anorm, bnorm;
    int maxit;
    complex shift;
    float tempr;
    complex ctemp2, ctemp3;
    int ilazr2;
    float ascale, bscale;
    complex signbc;
    extern double slamch_(char *), clanhs_(char *, int *, 
	    complex *, int *, float *);
    extern  int claset_(char *, int *, int *, complex 
	    *, complex *, complex *, int *), clartg_(complex *, 
	    complex *, float *, complex *, complex *);
    float safmin;
    extern  int xerbla_(char *, int *);
    complex eshift;
    int ilschr;
    int icompq, ilastm;
    complex rtdisc;
    int ischur;
    int ilazro;
    int icompz, ifirst, ifrstm, istart;
    int lquery;


/*  -- LAPACK routine (version 3.2) -- */
/*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
/*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/*     November 2006 */

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

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

/*  CHGEQZ computes the eigenvalues of a complex matrix pair (H,T), */
/*  where H is an upper Hessenberg matrix and T is upper triangular, */
/*  using the single-shift QZ method. */
/*  Matrix pairs of this type are produced by the reduction to */
/*  generalized upper Hessenberg form of a complex matrix pair (A,B): */

/*     A = Q1*H*Z1**H,  B = Q1*T*Z1**H, */

/*  as computed by CGGHRD. */

/*  If JOB='S', then the Hessenberg-triangular pair (H,T) is */
/*  also reduced to generalized Schur form, */

/*     H = Q*S*Z**H,  T = Q*P*Z**H, */

/*  where Q and Z are unitary matrices and S and P are upper triangular. */

/*  Optionally, the unitary matrix Q from the generalized Schur */
/*  factorization may be postmultiplied into an input matrix Q1, and the */
/*  unitary matrix Z may be postmultiplied into an input matrix Z1. */
/*  If Q1 and Z1 are the unitary matrices from CGGHRD that reduced */
/*  the matrix pair (A,B) to generalized Hessenberg form, then the output */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        Generalized Real Schur Form of a Regular Matrix Pair (A, B), in   
        M.S. Moonen et al (eds), Linear Algebra for Large Scale and   
        Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.   

    [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified   
        Eigenvalues of a Regular Matrix Pair (A, B) and Condition   
        Estimation: Theory, Algorithms and Software, Report   
        UMINF - 94.04, Department of Computing Science, Umea University,   
        S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.   
        To appear in Numerical Algorithms, 1996.   

    [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software   
        for Solving the Generalized Sylvester Equation and Estimating the   
        Separation between Regular Matrix Pairs, Report UMINF - 93.23,   
        Department of Computing Science, Umea University, S-901 87 Umea,   
        Sweden, December 1993, Revised April 1994, Also as LAPACK working   
        Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,   
        1996.   

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


       Decode and test the input parameters   

       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    
    /* 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;
    complex q__1, q__2;
    /* Builtin functions */
    double sqrt(doublereal), c_abs(complex *);
    void r_cnjg(complex *, complex *);
    /* Local variables */
    static integer kase, ierr;
    static real dsum;
    static logical swap;
    static integer i__, k;
    extern /* Subroutine */ int cscal_(integer *, complex *, complex *, 
	    integer *);
    static logical wantd;
    static integer lwmin;
    static logical wantp;
    static integer n1, n2;
    static logical wantd1, wantd2;
    static real dscale;
    static integer ks;
    extern /* Subroutine */ int clacon_(integer *, complex *, complex *, real 
	    *, integer *);
    extern doublereal slamch_(char *);
    static real rdscal;
    extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
	    *, integer *, complex *, integer *);
    static real safmin;
    extern /* Subroutine */ int ctgexc_(logical *, logical *, integer *, 
	    complex *, integer *, complex *, integer *, complex *, integer *, 
	    complex *, integer *, integer *, integer *, integer *), xerbla_(
	    char *, integer *), classq_(integer *, complex *, integer 
	    *, real *, real *);
    static integer liwmin;
    extern /* Subroutine */ int ctgsyl_(char *, integer *, integer *, integer 
	    *, complex *, integer *, complex *, integer *, complex *, integer 
	    *, complex *, integer *, complex *, integer *, complex *, integer 
	    *, real *, real *, complex *, integer *, integer *, integer *);

/* Subroutine */ int cchkgl_(integer *nin, integer *nout)
{
    /* Format strings */
    static char fmt_9999[] = "(\002 .. test output of CGGBAL .. \002)";
    static char fmt_9998[] = "(\002 ratio of largest test error             "
	    " = \002,e12.3)";
    static char fmt_9997[] = "(\002 example number where info is not zero   "
	    " = \002,i4)";
    static char fmt_9996[] = "(\002 example number where ILO or IHI is wrong"
	    " = \002,i4)";
    static char fmt_9995[] = "(\002 example number having largest error     "
	    " = \002,i4)";
    static char fmt_9994[] = "(\002 number of examples where info is not 0  "
	    " = \002,i4)";
    static char fmt_9993[] = "(\002 total number of examples tested         "
	    " = \002,i4)";

    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    real r__1, r__2, r__3;
    complex q__1;

    /* Builtin functions */
    integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_rsle(void);
    double c_abs(complex *);
    integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);

    /* Local variables */
    complex a[400]	/* was [20][20] */, b[400]	/* was [20][20] */;
    integer i__, j, n;
    complex ain[400]	/* was [20][20] */, bin[400]	/* was [20][20] */;
    integer ihi, ilo;
    real eps;
    integer knt, info, lmax[3];
    real rmax, vmax, work[120];
    integer ihiin, ninfo, iloin;
    real anorm, bnorm;
    extern /* Subroutine */ int cggbal_(char *, integer *, complex *, integer 
	    *, complex *, integer *, integer *, integer *, real *, real *, 
	    real *, integer *);
    extern doublereal clange_(char *, integer *, integer *, complex *, 
	    integer *, real *);
    real lscale[20];
    extern doublereal slamch_(char *);
    real rscale[20], lsclin[20], rsclin[20];

    /* Fortran I/O blocks */
    static cilist io___6 = { 0, 0, 0, 0, 0 };
    static cilist io___9 = { 0, 0, 0, 0, 0 };
    static cilist io___12 = { 0, 0, 0, 0, 0 };
    static cilist io___14 = { 0, 0, 0, 0, 0 };
    static cilist io___17 = { 0, 0, 0, 0, 0 };
    static cilist io___19 = { 0, 0, 0, 0, 0 };
    static cilist io___21 = { 0, 0, 0, 0, 0 };
    static cilist io___23 = { 0, 0, 0, 0, 0 };
    static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___35 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___36 = { 0, 0, 0, fmt_9997, 0 };
    static cilist io___37 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___38 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___39 = { 0, 0, 0, fmt_9994, 0 };
    static cilist io___40 = { 0, 0, 0, fmt_9993, 0 };



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

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

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

/*  CCHKGL tests CGGBAL, a routine for balancing a matrix pair (A, B). */

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

/*  NIN     (input) INTEGER */
/*          The logical unit number for input.  NIN > 0. */

/*  NOUT    (input) INTEGER */
/*          The logical unit number for output.  NOUT > 0. */

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

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Local Arrays .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
//.........这里部分代码省略.........

/* Subroutine */ int cdrvgb_(logical *dotype, integer *nn, integer *nval, 
	integer *nrhs, real *thresh, logical *tsterr, complex *a, integer *la, 
	 complex *afb, integer *lafb, complex *asav, complex *b, complex *
	bsav, complex *x, complex *xact, real *s, complex *work, real *rwork, 
	integer *iwork, integer *nout)
{
    /* Initialized data */

    static integer iseedy[4] = { 1988,1989,1990,1991 };
    static char transs[1*3] = "N" "T" "C";
    static char facts[1*3] = "F" "N" "E";
    static char equeds[1*4] = "N" "R" "C" "B";

    /* Format strings */
    static char fmt_9999[] = "(\002 *** In CDRVGB, LA=\002,i5,\002 is too sm"
	    "all for N=\002,i5,\002, KU=\002,i5,\002, KL=\002,i5,/\002 ==> In"
	    "crease LA to at least \002,i5)";
    static char fmt_9998[] = "(\002 *** In CDRVGB, LAFB=\002,i5,\002 is too "
	    "small for N=\002,i5,\002, KU=\002,i5,\002, KL=\002,i5,/\002 ==> "
	    "Increase LAFB to at least \002,i5)";
    static char fmt_9997[] = "(1x,a6,\002, N=\002,i5,\002, KL=\002,i5,\002, "
	    "KU=\002,i5,\002, type \002,i1,\002, test(\002,i1,\002)=\002,g12."
	    "5)";
    static char fmt_9995[] = "(1x,a6,\002( '\002,a1,\002','\002,a1,\002',"
	    "\002,i5,\002,\002,i5,\002,\002,i5,\002,...), EQUED='\002,a1,\002"
	    "', type \002,i1,\002, test(\002,i1,\002)=\002,g12.5)";
    static char fmt_9996[] = "(1x,a6,\002( '\002,a1,\002','\002,a1,\002',"
	    "\002,i5,\002,\002,i5,\002,\002,i5,\002,...), type \002,i1,\002, "
	    "test(\002,i1,\002)=\002,g12.5)";

    /* System generated locals */
    address a__1[2];
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9, i__10, 
	    i__11[2];
    real r__1, r__2;
    char ch__1[2];

    /* Builtin functions */
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
    /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
    double c_abs(complex *);

    /* Local variables */
    integer i__, j, k, n, i1, i2, k1, nb, in, kl, ku, nt, lda, ldb, ikl, nkl, 
	    iku, nku;
    char fact[1];
    integer ioff, mode;
    real amax;
    char path[3];
    integer imat, info;
    char dist[1];
    real rdum[1];
    char type__[1];
    integer nrun, ldafb;
    extern /* Subroutine */ int cgbt01_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, complex *, integer *, integer *, 
	    complex *, real *), cgbt02_(char *, integer *, integer *, integer 
	    *, integer *, integer *, complex *, integer *, complex *, integer 
	    *, complex *, integer *, real *), cgbt05_(char *, integer 
	    *, integer *, integer *, integer *, complex *, integer *, complex 
	    *, integer *, complex *, integer *, complex *, integer *, real *, 
	    real *, real *);
    integer ifact;
    extern /* Subroutine */ int cget04_(integer *, integer *, complex *, 
	    integer *, complex *, integer *, real *, real *);
    integer nfail, iseed[4], nfact;
    extern logical lsame_(char *, char *);
    char equed[1];
    integer nbmin;
    real rcond, roldc;
    extern /* Subroutine */ int cgbsv_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, integer *, complex *, integer *, 
	    integer *);
    integer nimat;
    real roldi;
    extern doublereal sget06_(real *, real *);
    real anorm;
    integer itran;
    logical equil;
    real roldo;
    char trans[1];
    integer izero, nerrs;
    logical zerot;
    char xtype[1];
    extern /* Subroutine */ int clatb4_(char *, integer *, integer *, integer 
	    *, char *, integer *, integer *, real *, integer *, real *, char *
), aladhd_(integer *, char *);
    extern doublereal clangb_(char *, integer *, integer *, integer *, 
	    complex *, integer *, real *), clange_(char *, integer *, 
	    integer *, complex *, integer *, real *);
    extern /* Subroutine */ int claqgb_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, real *, real *, real *, real *, 
	    real *, char *), alaerh_(char *, char *, integer *, 
	    integer *, char *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *);
    logical prefac;
    real colcnd;
    extern doublereal clantb_(char *, char *, char *, integer *, integer *, 
	    complex *, integer *, real *);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    Bjcol.ncol  = 1;
    Bjcol.Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
    if ( !Bjcol.Store ) SUPERLU_ABORT("SUPERLU_MALLOC fails for Bjcol.Store");
    Bjcol_store = Bjcol.Store;
    Bjcol_store->lda = ldb;
    Bjcol_store->nzval = work; /* address aliasing */
	
    /* Do for each right hand side ... */
    for (j = 0; j < nrhs; ++j) {
	count = 0;
	lstres = 3.;
	Bptr = &Bmat[j*ldb];
	Xptr = &Xmat[j*ldx];

	while (1) { /* Loop until stopping criterion is satisfied. */

	    /* Compute residual R = B - op(A) * X,   
	       where op(A) = A, A**T, or A**H, depending on TRANS. */
	    
#ifdef _CRAY
	    CCOPY(&A->nrow, Bptr, &ione, work, &ione);
#else
	    ccopy_(&A->nrow, Bptr, &ione, work, &ione);
#endif
	    sp_cgemv(transc, ndone, A, Xptr, ione, done, work, ione);

	    /* Compute componentwise relative backward error from formula 
	       max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )   
	       where abs(Z) is the componentwise absolute value of the matrix
	       or vector Z.  If the i-th component of the denominator is less
	       than SAFE2, then SAFE1 is added to the i-th component of the   
	       numerator before dividing. */

	    for (i = 0; i < A->nrow; ++i) rwork[i] = c_abs1( &Bptr[i] );
	    
	    /* Compute abs(op(A))*abs(X) + abs(B). */
	    if (notran) {
		for (k = 0; k < A->ncol; ++k) {
		    xk = c_abs1( &Xptr[k] );
		    for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
			rwork[Astore->rowind[i]] += c_abs1(&Aval[i]) * xk;
		}
	    } else {
		for (k = 0; k < A->ncol; ++k) {
		    s = 0.;
		    for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) {
			irow = Astore->rowind[i];
			s += c_abs1(&Aval[i]) * c_abs1(&Xptr[irow]);
		    }
		    rwork[k] += s;
		}
	    }
	    s = 0.;
	    for (i = 0; i < A->nrow; ++i) {
		if (rwork[i] > safe2) {
		    s = SUPERLU_MAX( s, c_abs1(&work[i]) / rwork[i] );
		} else if ( rwork[i] != 0.0 ) {
		    s = SUPERLU_MAX( s, (c_abs1(&work[i]) + safe1) / rwork[i] );
                }
                /* If rwork[i] is exactly 0.0, then we know the true 
                   residual also must be exactly 0.0. */
	    }
	    berr[j] = s;

	    /* Test stopping criterion. Continue iterating if   
	       1) The residual BERR(J) is larger than machine epsilon, and   

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