File: compile2.cpp

package info (click to toggle)
ginac 1.8.9-1
  • links: PTS
  • area: main
  • in suites: forky, sid
  • size: 6,640 kB
  • sloc: cpp: 49,195; sh: 5,402; makefile: 448; python: 193
file content (67 lines) | stat: -rw-r--r-- 1,686 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include <iostream>
#ifdef IN_GINAC
#include "ginac.h"
#else
#include <ginac/ginac.h>
#endif
// Yes, we are using CUBA (should be installed on the system!)
#include <cuba.h>

using namespace std;
using namespace GiNaC;

/*
 * Demonstrates the use of compile_ex with the CUBA library.
 *
 * The user can enter an expression on the command line. This expression is
 * compiled via compile_ex and integrated over the region 0 <= x,y <= 1 with
 * the help of the CUBA library (http://www.feynarts.de/cuba).
 *
 */

int main()
{
	// Let the user enter a expression
	symbol x("x"), y("y");
	string s;
	cout << "Enter an expression containing 'x' and/or 'y': ";
	cin >> s;
	// Expression now in expr
	ex expr(s, lst{x,y});
 
	cout << "start integration of " << expr << " ..." << endl;

	// Some definitions for VEGAS 
	#define NDIM 2
	#define NCOMP 1
	#define EPSREL 1e-3
	#define EPSABS 1e-12
	#define VERBOSE 0
	#define MINEVAL 0
	#define MAXEVAL 50000
	#define NSTART 1000
	#define NINCREASE 500

	// Some variables for VEGAS
	int comp, nregions, neval, fail;
	double integral[NCOMP], error[NCOMP], prob[NCOMP];

	// Our function pointer that points to the compiled ex
	FUNCP_CUBA fp;

	// Optionally, compile with custom compiler flags:
	// setenv("CXXFLAGS", "-O3 -fomit-frame-pointer -ffast-math", 1);
	compile_ex(lst{expr}, lst{x,y}, fp);

	// Starting VEGAS
	// By invocation of compile() the expression in expr is converted into the
	// appropriate function pointer
	Vegas(NDIM, NCOMP, fp,
	      EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL, NSTART, NINCREASE,
	      &neval, &fail, integral, error, prob);

	// Show the result
	cout << "result: " << integral[0] << endl;

	return 0;
}