fastgauss.cpp

快速高斯变换的程序

// Implementation for Fast Gauss Transform
#include "FastGauss.h"
//#include "mkl.h"
//#include <iostream.h>
#include <math.h>
#include <limits.h>
//#include "mex.h"



////////////////////////////////////////////////////////////////////
// Constructor 1
//
// PURPOSE                                                    
// -------   
// Initialize the class without the parameter tol. 
// Allocate memory for the fast Gauss algorithm
// Determine the memory space needed. 
// Compute the coefficients of the Taylor expansions.
//
// PARAMETERS                                                      
// ----------
// dim		--> Dimension of problem
// u		--> Weight matrix with du rows and ns columns
// du		--> Row dimension of matrix u (default = 1)
// x		--> Source matrix with dim rows and ns columns
// ns		--> Column dimension of matrix x (no. of sources)
// h		--> Bandwidth
// p		--> Order of truncation
// K		--> Number of centers
// e		--> Ratio of far field
//
// OUTPUTS
// -------
// xc		<--	center of the sources
// A_k		<-- Coefficients of the Taylor expansion with wts u
////////////////////////////////////////////////////////////////////
FastGauss::FastGauss(int dim, double *u, double *x, int ns, 
					double h, int p, int K, double e)
{	
	// interface:
	Dim = dim;
	NSources = ns;	
	pWeights = u;
	pSources = x;
	bandwid = h;
	pterms = p;
	K_N = K;
	ratio_far = e;
		
	pd = nchoosek(pterms+Dim-1, Dim); //pd = C_dim^(dim+pterms-1)
	
	A_k = new double[pd*K_N];

	xc = new double[K_N*Dim];
	xcind = new int[K_N];
	indx = new int[NSources];
	xheads = new int[K_N];
	xboxsz = new int[K_N];
	
	// Divide the source space into K_N parts using K-center algorithm
	rx2 = dkcenter(Dim,NSources,K_N,pSources,xcind,indx,xboxsz);


	// computer the center of the sources
	Compute_Centers(xc,pSources,Dim,NSources,indx,K_N,xboxsz);

	// Build hashing tables for the nearest neighbor searching
	 
	TaylorExpansion(Dim, pterms);

}


//////////////////////////////////////////////////////////////
// Destructor
//
// PURPOSE                                                    
// -------   
// Free the memory.
//////////////////////////////////////////////////////////////
FastGauss::~FastGauss()
{
	delete []A_k;
	delete []xc;
	delete []xcind;
	delete []indx;
	delete []xheads;
	delete []xboxsz;
}


//////////////////////////////////////////////////////////////
// n choose k
//
// PURPOSE                                                    
// -------   
// Compute the combinatorial number $C_k^n$.
//////////////////////////////////////////////////////////////
int FastGauss::nchoosek(int n, int k)
{
	int n_k = n - k;
	
	if (k < n_k)
	{
		k = n_k;
		n_k = n - k;
	}

	int  nchsk = 1; 
	for ( int i = 1; i <= n_k; i++)
	{
		nchsk *= (++k);
		nchsk /= i;
	}

	return nchsk;
}



////////////////////////////////////////////////////////////////////
// Taylor Expansion
//
// PURPOSE                                                    
// -------   
// Compute the Taylor expansions from the sources.
//
// PARAMETERS                                                      
// ----------
// d		--> Dimension of problem
// p		--> Order of truncation
//
// OUTPUTS
// -------
// xc		<--	center of the sources
// A_k		<-- Coefficients of the Taylor expansion with wts u
// B_k		<-- Coefficients of the Taylor expansion without wts
////////////////////////////////////////////////////////////////////
void FastGauss::TaylorExpansion(int d, int p)
{

	// compute C_k(p, p)
	double *C_k = new double[pd];
	Compute_C_k(d, p, C_k);

	// compute A_k(p, p) and B_k(p, p)	
	Compute_A_k(xc, C_k);

	delete[] C_k;
	
	return;
	
}


void FastGauss::Compute_Centers(double *xc, double *x, 
								int d, int n, int *ind, int K, int *bxsz)
{

	for (int i = 0; i < d*K; i++)
		xc[i] = 0.0;

	
	for (int i = 0, nd = 0; i < n; i++,nd+=d)
	{
		int ibase = ind[i]*d;
		for (int j = 0; j < d; j++)
			xc[ibase+j] += x[nd+j];
	}

	for (int i = 0, ibase = 0; i < K; i++,ibase+=d)
		for (int j = 0; j < d; j++)
			xc[ibase+j] /= bxsz[i];

}

////////////////////////////////////////////////////////////////////
// Compute Taylor Expansion Coefficients
//
// PURPOSE                                                    
// -------   
// Compute the coefficients of the Taylor expansions.
//
// PARAMETERS                                                      
// ----------
// d		--> Dimension of problem
// p		--> Order of truncation
//
// OUTPUTS
// -------
// C_k		<-- Coefficients of the Taylor expansion
////////////////////////////////////////////////////////////////////
void FastGauss::Compute_C_k(int d, int p, double *C_k)
{
	
	int *heads = new int[d+1];
	int *cinds = new int[pd];
	
	for (int i = 0; i < d; i++)
		heads[i] = 0;
	heads[d] = INT_MAX;
	
	cinds[0] = 0;
	C_k[0] = 1.0;
	for (int k=1, t=1, tail=1; k < p; k++, tail=t)
	{
		for (int i = 0; i < d; i++)
		{
			int head = heads[i];
			heads[i] = t;
			for ( int j = head; j < tail; j++, t++)
			{
				cinds[t] = (j < heads[i+1])? cinds[j] + 1 : 1;
				C_k[t] = 2.0 * C_k[j];
				C_k[t] /= (double) cinds[t];
			}
		}
	}
	
	delete []cinds;
	delete []heads;
	return;
}


////////////////////////////////////////////////////////////////////
// Compute Taylor Expansion
//
// PURPOSE                                                    
// -------   
// Compute the Taylor expansions from the sources.
//
// PARAMETERS                                                      
// ----------
// xc		--> Center of the sources
// C_k		--> Coefficents of the Taylor expansion
//
// OUTPUTS
// -------
// A_k		<-- Coefficients of the Taylor expansion with wts u
////////////////////////////////////////////////////////////////////
void FastGauss::Compute_A_k(double *xc, double *C_k)
{
	double *dx = new double[Dim];
	double *prods = new double[pd];
	
	int *heads = new int[Dim];

	for (int i = 0; i < pd*K_N; i++)
		A_k[i] = 0.0;

		
	for (int n=0; n < NSources; n++)
	{
		int nbase = n*Dim;
		int ix2c = indx[n];
		int ix2cbase = ix2c*Dim;
		double sum = 0.0;
		for (int i = 0; i < Dim; i++)
		{
			dx[i] = (pSources[nbase+i] - xc[ix2cbase+i])/bandwid;
			sum -= dx[i] * dx[i];
			heads[i] = 0;
		}
		
		prods[0] = exp(sum);
		for (int k=1, t=1, tail=1; k < pterms; k++, tail=t)
		{
			for (int i = 0; i < Dim; i++)
			{
				int head = heads[i];
				heads[i] = t;
				for ( int j = head; j < tail; j++, t++)
					prods[t] = dx[i] * prods[j];
			} // for i
		}// for k
#ifdef _MKL_H_
		cblas_daxpy(pd, pWeights[n], prods, 1, A_k+ix2c*pd, 1);
#else
		for (int i = 0; i < pd; i++)
			A_k[ix2c*pd+i] += pWeights[n]*prods[i];
#endif		
	}// for n
	
	for (int k = 0; k < K_N; k++)
	{
		for (int i=0; i < pd; i++)
			A_k[k*pd+i] *= C_k[i];
	}

	delete []dx;
	delete []heads;
	delete []prods;

	return;
}
	



////////////////////////////////////////////////////////////////////
// Compute the summation
//
// PURPOSE                                                    
// -------   
// Compute the summatioin of the Taylor expansions from the sources.
//
// PARAMETERS                                                      
// ----------
// y		--> Target points
// mt		--> Column dimension of the target points
//
// OUTPUTS
// -------
// vi		<-- Results (1-by-mt vector)
////////////////////////////////////////////////////////////////////
void FastGauss::Compute_v_i(double *y, int mt, double *v_i)
{
	double *dy = new double[Dim];
	double *prods = new double[pd];
	int *heads = new int[Dim];

	for (int m=0; m < mt; m++)
	{	
		v_i[m] = 0.0;

		int mbase = m*Dim;
		double sumgx = 0.0;
		
		for (int kn=0; kn < K_N; kn++)
		{
			int xbase = kn*Dim;
			double sum2 = 0.0;
			for (int i = 0; i < Dim; i++)
			{
				dy[i] = (y[mbase+i] - xc[xbase+i])/bandwid;
				sum2 += dy[i] * dy[i];
			}
		
			if (sum2 > ratio_far) continue; //skip to next kn

			for (int i = 0; i < Dim; i++)
				heads[i] = 0;

			prods[0] = exp(-sum2);		
			for (int k=1, t=1, tail=1; k < pterms; k++, tail=t)
			{
				for (int i = 0; i < Dim; i++)
				{
					int head = heads[i];
					heads[i] = t;
					for ( int j = head; j < tail; j++, t++)
						prods[t] = dy[i] * prods[j];
				} // for i
			}// for k
#ifdef _MKL_H_
			v_i[m] += cblas_ddot(pd, prods, 1, A_k+kn*pd, 1);
#else
			for (int i = 0; i < pd; i++)
				v_i[m] += A_k[kn*pd+i]*prods[i];
#endif	
		}// for kn
	}//for m

	delete []dy;
	delete []prods;
	delete []heads;

	return;
}