C++ dvector函数代码示例

NUMERICS_EXPORT BOOL inverse(double **a, int n)
{
    double d;
    int i, j;
    BOOL ret = FALSE;
    double** ai = dmatrix(0, n - 1, 0, n - 1);
    double* col = dvector(0, n - 1);
    int* indx = ivector(0, n - 1);
    
    if(ludcmp(a, n, indx, &d)){
        for(j = 0; j < n; j++){
            for(i = 0; i < n; i++) col[i] = 0.0;
            col[j] = 1.0;
            lubksb(a, n, indx, col);
            for(i = 0; i < n; i++) ai[i][j] = col[i];
        }
        for(i = 0; i < n; i++){
            for(j = 0; j < n; j++){
                a[i][j] = ai[i][j];
            }
        }
        ret = TRUE;
    }
    
    free_dmatrix(ai, 0, n - 1, 0);
    free_dvector(col, 0);
    free_ivector(indx, 0);
    
    return ret;
}

void Tri::dump_mesh(FILE *name){
  Coord X;
  X.x = dvector(0,QGmax*QGmax-1);
  X.y = dvector(0,QGmax*QGmax-1);

  this->coord(&X);

  fprintf(name,"VARIABLES = x y z\n");
  fprintf(name,"ZONE T=\"Element %d\", I=%d, J=%d, F=POINT\n",
    id,qa,qb);

  for(int i=0;i<qa*qb;++i)
    fprintf(name,"%lf %lf %lf\n", X.x[i], X.y[i], h[0][i]);
  free(X.x);
  free(X.y);
}

void ReadHMM(FILE *fp, HMM *phmm)
{
	int i, j, k;

	fscanf(fp, "M= %d\n", &(phmm->M)); 

	fscanf(fp, "N= %d\n", &(phmm->N)); 

	fscanf(fp, "A:\n");
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);
	for (i = 1; i <= phmm->N; i++) { 
		for (j = 1; j <= phmm->N; j++) {
			fscanf(fp, "%lf", &(phmm->A[i][j])); 
		}
		fscanf(fp,"\n");
	}

	fscanf(fp, "B:\n");
	phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);
	for (j = 1; j <= phmm->N; j++) { 
		for (k = 1; k <= phmm->M; k++) {
			fscanf(fp, "%lf", &(phmm->B[j][k])); 
		}
		fscanf(fp,"\n");
	}

	fscanf(fp, "pi:\n");
	phmm->pi = (double *) dvector(1, phmm->N);
	for (i = 1; i <= phmm->N; i++) 
		fscanf(fp, "%lf", &(phmm->pi[i])); 

}

/*!
    \fn SGLPolygonObj::getCenter()
 */
SGLVektor SGLPolygonObj::getCenter()const
{
	SGLVektor ret;
	for(int i=0;i<Fl.Cnt;i++)
		ret+=Fl.Fl[i]->getCenter();
	return dvector( ret/Fl.Cnt);
}

/*
 *    Given a, b, and c as output from chebyshev_fit, and given m, the desired degree of
 *    approximation, this routine returns the array d[0..m-1], of coefficients of a polynomial
 *    expansion which is equivalent to the Chebyshev fit.
 */
void chebyshev_2_poly( double a, double b, double *c, double *d, int m ) {

    int j, k;
    double sv, *dd;

    dd = dvector( 0, m-1 );

    for ( j = 0; j < m; ++j ) d[j] = dd[j] = 0.0;
    d[0] = c[m-1];
    for ( j = m-2; j >= 1; --j ) {
        for ( k = m-j; k >= 1; --k ) {
            sv = d[k];
            d[k] = 2.0*d[k-1] - dd[k];
            dd[k] = sv;
        }
        sv = d[0];
        d[0] = -dd[0] + c[j];
        dd[0] = sv;
    }
    for ( j = m-1; j >= 1; --j ) d[j] = d[j-1] - dd[j];
    d[0] = -dd[0] + 0.5*c[0];

    free_dvector( dd, 0, m-1 );

    /*
     *    Map the interval [-1,+1] to [a,b].
     */
    poly_shift_coeff( a, b, d, m );
}

void Newton_Solver(  double * x1,  double ** A,  double * b, int vn)
{
	int i, j;
	double ger;

	int *p = ivector( 1, vn);

	double  **LU = dmatrix( 1, vn, 1, vn);

		for( i=0; i<=vn-1; i++)
		    for( j=0; j<=vn-1; j++)
			LU[i+1][j+1] = A[i][j];

	double *X = dvector( 1, vn);
		for( i=0; i<=vn-1; i++)
		    X[i+1] = b[i];

 	ludcmp( LU, vn, p, &ger);

	lubksb( LU, vn, p, X );

	for( i=0; i<=vn-1; i++)
	    x1[i] = X[i+1];

}

void Matrix_Inverse( double **invMat, double **Mat, int nn)
{
        int i, j;
        double ger;

        int *p= ivector( 1, nn);


        double **LU = dmatrix( 1, nn, 1, nn);
                for( i=0; i<=nn-1; i++)
                    for( j=0; j<=nn-1; j++)
                        LU[i+1][j+1] = Mat[i][j];

        double *col = dvector( 1, nn);


        ludcmp( LU, nn, p, &ger);

	for( j=1; j<=nn; j++) {
	    for( i=1; i<=nn; i++) col[i] = 0.0;
 	    col[j] = 1.0;
	    lubksb( LU, nn, p, col );
	    for( i=1; i<=nn; i++) invMat[i-1][j-1] = (double) col[i];
			     };

	free_dvector( col, 1, nn);
	free_ivector( p, 1, nn);
	free_dmatrix( LU, 1, nn, 1, nn);

}

double  *fit_poly_norm (double **x, double *y, int ndata, int npol)
{
  double chisq, *sig, *w, **u, **v, *a;
  double *dvector(), **dmatrix();
  int  i;
  void svdfit();

  a = dvector(1,npol);
  w = dvector(1,npol);
  sig = dvector(1,ndata);
  u = dmatrix(1,ndata,1,npol);
  v = dmatrix(1,npol,1,npol);
  for (i=1;i<=ndata;i++) sig[i] = 1.0;
  svdfit(x, y, sig, ndata, a, npol, u, v, w, &chisq, fpoly);
  return a;
}

/* Invert a double matrix, 1-indexed of size dim
 * from Numerical Recipes.  Input matrix is destroyed.
 */
int
invert_matrix(double **in, double **out,
              int dim)
{
  extern void ludcmp(double **a, int n, int *indx, double *d);
  extern void ludcmp(double **a, int n, int *indx, double *d);

  int   *indx,i,j;
  double *tvec,det;

  tvec = dvector(1,dim);
  indx = ivector(1,dim);
  ludcmp(in,dim,indx,&det);

  for (j=1; j<=dim; j++) {
    for (i=1; i<=dim; i++) tvec[i]=0.;
    tvec[j] = 1.0;
    lubksb(in,dim,indx,tvec);
    for (i=1; i<=dim; i++) out[i][j]=tvec[i];
  }

  free_ivector(indx,1,6);
  free_dvector(tvec,1,6);
  return(0);
}

/* compute the slope vector dy for the transient equation
 * dy + cy = p. useful in the transient solver
 */
void slope_fn_block(block_model_t *model, double *y, double *p, double *dy)
{
	/* shortcuts	*/
	int n = model->n_nodes;
	double **c = model->c;

	/* for our equation, dy = p - cy */
	#if (MATHACCEL == MA_INTEL || MATHACCEL == MA_APPLE)
	/* dy = p	*/
	cblas_dcopy(n, p, 1, dy, 1);
	/* dy = dy - c*y = p - c*y */
	cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n, -1, c[0],
				n, y, 1, 1, dy, 1);
	#elif (MATHACCEL == MA_AMD || MATHACCEL == MA_SUN)
	/* dy = p	*/
	dcopy(n, p, 1, dy, 1);
	/* dy = dy - c*y = p - c*y */
	dgemv('T', n, n, -1, c[0], n, y, 1, 1, dy, 1);
	#else
	int i;
	double *t = dvector(n);
	matvectmult(t, c, y, n);
	for (i = 0; i < n; i++)
		dy[i] = p[i]-t[i];
	free_dvector(t);
	#endif
}

/*
CALCULATES GEOMETRIC DISTRIBUTION. 
NEEDS:
    - M (slope in log-linear space)
    - S (number of taxa)
RETURNS:
	- abundance: proportional abundances
***********************************************************************/
double *proportional_gs_distribution(double M, int S)
{
int i = 0;							/* loop variable */
double sum = 0.0f;					/* number of "occurrences" */
double *A;							/* array to return */

if (M<1.0f) {
	printf("\nproportional_gs_distribution, illegal slope = %f\n",M);
	exit(1);
	}
if (S<=0)	{
	printf("\nproportional_gs_distribution, illegal number of taxa = %d\n",S);
	exit(1);
	}

A = dvector(S);						/* allocate array */

sum = A[0] = 100;					/* taxon 0 = 100 occurrences */
for (i=1; i<S; i++) {
	A[i] = A[i-1] / M;				/* taxon i = taxon (i-1) / slope */
	sum += A[i];					/* sum number of occurrences */
	}

for (i=0; i<S; i++) A[i] /=  sum;	/* make proportional */

return A;
}

int Discount_Factors_opt(FTYPE *pdDiscountFactors, 
		     int iN, 
		     FTYPE dYears, 
		     FTYPE *pdRatePath)
{
	int i,j;				//looping variables
	int iSuccess;			//return variable
	
	FTYPE ddelt;			//HJM time-step length
	ddelt = (FTYPE) (dYears/iN);

	FTYPE *pdexpRes;
	pdexpRes = dvector(0,iN-2);

	//initializing the discount factor vector
	for (i=0; i<=iN-1; ++i)
	  pdDiscountFactors[i] = 1.0;

	//precompute the exponientials
	for (j=0; j<=(i-2); ++j){ pdexpRes[j] = -pdRatePath[j]*ddelt; }
	for (j=0; j<=(i-2); ++j){ pdexpRes[j] = exp(pdexpRes[j]);  }
	
	for (i=1; i<=iN-1; ++i)
	  for (j=0; j<=i-1; ++j)
	    pdDiscountFactors[i] *= pdexpRes[j];
		
	free_dvector(pdexpRes, 0, iN-2);
	iSuccess = 1;
	return iSuccess;
}

/* Deserializes a matrix from R into a vector.
   Note that the "top" offset needs to already be applied to "va".
 */
double *Runpack_dvectors(double *va, unsigned int n, double *a, unsigned int sample_size){
  if(!a) a=dvector(n);
  unsigned int i;
  for(i=0;i<n;i++)
    a[i]=va[sample_size*i];
  return a;
}

void pzextr(int iest, double xest, double yest[], double yz[], double dy[], int nv)
{
    int k1,j;
    double q,f2,f1,delta,*c;

    c=dvector(1,nv);
    x[iest]=xest;
    for (j=1;j<=nv;j++) dy[j]=yz[j]=yest[j];
    if (iest == 1) {
        for (j=1;j<=nv;j++) d[j][1]=yest[j];
    } else {
        for (j=1;j<=nv;j++) c[j]=yest[j];
        for (k1=1;k1<iest;k1++) {
            delta=1.0/(x[iest-k1]-xest);
            f1=xest*delta;
            f2=x[iest-k1]*delta;
            for (j=1;j<=nv;j++) {
                q=d[j][k1];
                d[j][k1]=dy[j];
                delta=c[j]-q;
                dy[j]=f1*delta;
                c[j]=f2*delta;
                yz[j] += dy[j];
            }
        }
        for (j=1;j<=nv;j++) d[j][iest]=dy[j];
    }
    free_dvector(c,1,nv);
}

void vander(double x[], double w[], double q[], int n)
{
	int i,j,k,k1;
	double b,s,t,xx;
	double *c;

	c=dvector(1,n);
	if (n == 1) w[1]=q[1];
	else {
		for (i=1;i<=n;i++) c[i]=0.0;
		c[n] = -x[1];
		for (i=2;i<=n;i++) {
			xx = -x[i];
			for (j=(n+1-i);j<=(n-1);j++) c[j] += xx*c[j+1];
			c[n] += xx;
		}
		for (i=1;i<=n;i++) {
			xx=x[i];
			t=b=1.0;
			s=q[n];
			k=n;
			for (j=2;j<=n;j++) {
				k1=k-1;
				b=c[k]+xx*b;
				s += q[k1]*b;
				t=xx*t+b;
				k=k1;
			}
			w[i]=s/t;
		}
	}
	free_dvector(c,1,n);
}