C++ c_abs函数代码示例

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