Commit aab5e808 authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

GslWrapper additions 

integration over semi-infinite interval and
Log of Gamma Function as given by the GSL
parent 39656742

lib/configure.pl

+2 −1
Original line number Diff line number Diff line
@@ -40,7 +40,8 @@ sub createMakefile 


	local *FH = $fh;

	my @units = ("MersenneTwister","Matrix","Mpi","ApplicationInfo","Concurrency",

	"ProgressIndicator","Tokenizer","MemResolv","PsimagLite","PsiBase64");

	"ProgressIndicator","Tokenizer","MemResolv","PsimagLite","PsiBase64",

	"GammaFunction");

	my $combinedUnits = combine("",\@units,".o ");

	my $combinedUnits2 = combine("../src/",\@units,".cpp ");

src/GslWrapper.h

+56 −10
Original line number Diff line number Diff line
@@ -82,6 +82,8 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS. 
#ifdef USE_GSL

#include <gsl/gsl_integration.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_result.h>

#include <gsl/gsl_sf_gamma.h>

#endif



#include <stdexcept>
@@ -95,6 +97,7 @@ public: 
	typedef int DummyType;

	typedef DummyType gsl_integration_workspace;

	typedef double (* GslWrapperFunctionType) (double, void * );

	typedef double gsl_sf_result;



	typedef void (* gsl_error_handler_t) (const char *,const char *,int,int);
@@ -139,6 +142,19 @@ public: 
		return 0;

	}



	int gsl_integration_qagiu(gsl_function*,

	                          double,

	                          double,

	                          double,

	                          size_t ,

	                          gsl_integration_workspace*,

	                          double*,

	                          double*) const

	{

		thereSnoGsl();

		return 0;

	}



	int gsl_integration_qagp (const gsl_function*,

	                          double*,

	                          SizeType,
@@ -170,6 +186,15 @@ public: 
		thereSnoGsl();

	}



	int gsl_sf_lngamma_complex_e(double zr,

	                             double zi,

	                             gsl_sf_result* lnr,

	                             gsl_sf_result* arg) const

	{

		thereSnoGsl();

		return 0;

	}



private:

	void thereSnoGsl() const

	{
@@ -217,6 +242,18 @@ public: 
		return ::gsl_integration_qagi(f,epsabs,epsrel,limit,workspace,result,abserr);

	}



	int gsl_integration_qagiu(gsl_function* f,

	                          double a,

	                          double epsabs,

	                          double epsrel,

	                          size_t limit,

	                          gsl_integration_workspace* workspace,

	                          double* result,

	                          double* abserr) const

	{

		return ::gsl_integration_qagiu(f,a,epsabs,epsrel,limit,workspace,result,abserr);

	}



	int gsl_integration_qagp (const gsl_function * f,

	                          double * pts,

	                          SizeType npts,
@@ -260,6 +297,15 @@ public: 
		                             result,

		                             abserr);

	}



	int gsl_sf_lngamma_complex_e(double zr,

	                             double zi,

	                             gsl_sf_result* lnr,

	                             gsl_sf_result* arg) const

	{

		return ::gsl_sf_lngamma_complex_e(zr, zi, lnr, arg);

	}



}; // class GslWrapper



#endif

src/Integrator.h

+24 −6
Original line number Diff line number Diff line
@@ -51,6 +51,24 @@ public: 
		}

	}



	RealType toInfinity(RealType a,

	                    IntegrationEnum integ = INTEG_QAG,

	                    int key = 4)

	{

		int status = gslWrapper_.gsl_integration_qagiu(&f_,

		                                               a,

		                                               epsabs_,

		                                               epsrel_,

		                                               limit_,

		                                               workspace_,

		                                               &result_,

		                                               &abserr_);

		if (status)

			gslWrapper_.printError(status);



		return result_;

	}



	RealType operator()()

	{

		int status = gslWrapper_.gsl_integration_qagi(&f_,