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

void get_hmm_from_file(FILE *fp, HMM *hmm_ptr){

	int i, j, num_trans, num_state, from, to;
	double pr;

	fscanf(fp, "Symbol= %d\n", &(hmm_ptr->M));
	fscanf(fp, "State= %d\n", &(hmm_ptr->N));

	/* default transition */
	hmm_ptr->A = (double **)dmatrix(hmm_ptr->N, hmm_ptr->N);
	for (i=0; i<hmm_ptr->N; i++){
		for (j=0; j<hmm_ptr->N; j++){
			hmm_ptr->A[i][j] = 0;
		}
	}

	/* transition */
	fscanf(fp, "Transition= %d\n", &num_trans);
	for (i=0; i<num_trans; i++){
		fscanf(fp, "%d %d %lf\n", &from, &to, &pr);
		hmm_ptr->A[from][to] = pr;
	}

	/* start state*/
	fscanf(fp, "Pi= %d\n", &num_state);
	hmm_ptr->pi = (double *)dvector(hmm_ptr->N);
	for (i=0; i<hmm_ptr->N; i++){
		fscanf(fp, "%lf\n", &pr);
		hmm_ptr->pi[i] = pr;
	}
}

void CopyHMM(HMM *phmm1, HMM *phmm2)
{
        int i, j, k;
 
        phmm2->M = phmm1->M;

 
        phmm2->N = phmm1->N;
 
        phmm2->A = (double **) dmatrix(1, phmm2->N, 1, phmm2->N);
 
        for (i = 1; i <= phmm2->N; i++)
                for (j = 1; j <= phmm2->N; j++)
                        phmm2->A[i][j] = phmm1->A[i][j];
 
        phmm2->B = (double **) dmatrix(1, phmm2->N, 1, phmm2->M);
        for (j = 1; j <= phmm2->N; j++)
                for (k = 1; k <= phmm2->M; k++)
                        phmm2->B[j][k] = phmm1->B[j][k];
 
        phmm2->pi = (double *) dvector(1, phmm2->N);
        for (i = 1; i <= phmm2->N; i++)
                phmm2->pi[i] = phmm1->pi[i]; 
 
}

void svbksb(double **u, double *w, double **v, int m, int n, double *b, double *x)
{
  int jj,j,i;
  double s,*tmp,*dvector();
  void free_dvector();

  tmp=dvector(1,n);
  for (j=1;j<=n;j++)
    {  /* calculate U(transpose)B */
      s=0.0;
      if (w[j])
	{
	  for (i=1;i<=m;i++) s += u[i][j]*b[i];
	  s /= w[j];
	}
      tmp[j]=s;
    }
  for (j=1;j<=n;j++)
    {
      s=0.0;
      for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
      x[j]=s;
    }
  free_dvector(tmp,1,n);
}

/*
CALCULATES NORMAL DISTRIBUTION HISTOGRAM WITH SYMMETRICAL BARS AROUND MEAN,
	THIS MEANS THAT THE HEIGHT 
NEEDS:
    - num_stdev - Number of divisions along normal curve to be used
    - oct_per_stdev - Octaves per SD (1 makes 1 Octave = 1 SD; 2 makes 2 Octaves = 1 SD
    - modal_oct - Octave with "mean" (=Octaves to the left of the "mean" on the normal curve.)
		NOTE: I am assuming that if someone enters "5," then they mean the fifth element, which is element 4;
			hence, we used modal_oct-1 now.  
RETURNS:
    - ARRAY GIVING NORMAL DISTRIBUTION BEGINING [START] SD's BEFORE THE MEAN
*************************************************************/
double* normdistevn(int num_stdev, int oct_per_stdev, int modal_oct) {
	int	a;
	int length = oct_per_stdev * num_stdev;
	double y;
	double Oct = 0.5+(-1.0 * (double) (modal_oct-1) / (double) oct_per_stdev);	/* changed to modal_oct-1 to accomodate start at zero */
	double Area = 0.0f;
	double *NormA;
	double p = pow(2*PI,0.5);

	NormA=dvector(length);

	/* Find the height of the histogram for each octave 		*/
	/* This will be used to determine the proportion of species */
	/* that fall into a category								*/
	for (a=0 ; a<length ; a++)	{
		y = exp(-(Oct*Oct)/2);
		y/= p;
		NormA[a] = y;
		Area+= y;
		Oct+= 1.0/ (double) oct_per_stdev;
		}
	
	for (a=0; a<(length); a++)
		NormA[a] /= Area;
	
	return NormA;
}

void Tri::genFile (Curve *cur, double *x, double *y){
  register int i;
  Geometry *g;
  Point    p1, p2, a;
  double   *z, *w, *eta, xoff, yoff;
  int      fac;

  fac = cur->face;

  p1.x = vert[vnum(fac,0)].x;  p1.y = vert[vnum(fac,0)].y;
  p2.x = vert[vnum(fac,1)].x;  p2.y = vert[vnum(fac,1)].y;

  getzw(qa,&z,&w,'a');

  eta    = dvector (0, qa);
  if ((g = lookupGeom (cur->info.file.name)) == (Geometry *) NULL)
       g = loadGeom   (cur->info.file.name);


  /* If the current edge has an offset, apply it now */
  xoff = cur->info.file.xoffset;
  yoff = cur->info.file.yoffset;
  if (xoff != 0.0 || yoff != 0.0) {
    dsadd (g->npts, xoff, g->x, 1, g->x, 1);
    dsadd (g->npts, yoff, g->y, 1, g->y, 1);
    if (option("verbose") > 1)
      printf ("shifting current geometry by (%g,%g)\n", xoff, yoff);
  }

  /* get the end points which are assumed to lie on the curve */
  /* set up search direction in normal to element point -- This
     assumes that vertices already lie on spline */

  a.x      = p1.x  - (p2.y - p1.y);
  a.y      = p1.y  + (p2.x - p1.x);
  eta[0]   = searchGeom (a, p1, g);
  a.x      = p2.x  - (p2.y - p1.y);
  a.y      = p2.y  + (p2.x - p1.x);
  eta[qa-1] = searchGeom (a, p2, g);

  /* Now generate the points where we'll evaluate the geometry */

  for (i = 1; i < qa-1; i++)
    eta [i] = eta[0] + 0.5 * (eta[qa-1] - eta[0]) * (z[i] + 1.);

  for (i = 0; i < qa; i++) {
    x[i] = splint (g->npts, eta[i], g->arclen, g->x, g->sx);
    y[i] = splint (g->npts, eta[i], g->arclen, g->y, g->sy);
  }

  g->pos = 0;     /* reset the geometry */
  if (xoff != 0.)
    dvsub (g->npts, g->x, 1, &xoff, 0, g->x, 1);
  if (yoff != 0.)
    dvsub (g->npts, g->y, 1, &yoff, 0, g->y, 1);

  free (eta);    /* free the workspace */

  return;
}