Commit 85f78394 authored by Martin Robinson's avatar Martin Robinson
Browse files

add everything

parents
# location of the Boost Python include files and library
BOOST_INC = /usr/include
BOOST_LIB = /usr/lib
printRank: printRank.o
/usr/bin/mpic++ -o printRank printRank.o -L$(BOOST_LIB) -lboost_mpi -lboost_serialization
printRank.o: printRank.cpp
/usr/bin/mpic++ -I$(BOOST_INC) -c printRank.cpp
#include <boost/mpi/environment.hpp>
#include <boost/mpi/communicator.hpp>
#include <iostream>
namespace mpi = boost::mpi;
int main()
{
mpi::environment env;
mpi::communicator world;
std::cout << "I am process " << world.rank() << " of " << world.size()
<< "." << std::endl;
return 0;
}
# location of the Boost Python include files and library
BOOST_INC = /usr/include
BOOST_LIB = /usr/lib
common_factor: common_factor.o
g++ -o common_factor common_factor.o
common_factor.o: common_factor.cpp
g++ -I$(BOOST_INC) -c common_factor.cpp
#include <boost/math/common_factor.hpp>
#include <iostream>
int main()
{
using std::cout;
using std::endl;
cout << "The GCD and LCM of 6 and 15 are "
<< boost::math::gcd(6, 15) << " and "
<< boost::math::lcm(6, 15) << ", respectively."
<< endl;
}
# location of the Boost include files and library
BOOST_INC = /usr/include
BOOST_LIB = /usr/lib
harmonic_oscillator: harmonic_oscillator.o
g++ -o $@ -L$(BOOST_LIB) $^
%.o: %.cpp
g++ -I$(BOOST_INC) -c $^ -o $@
clean:
rm *.o
rm harmonic_oscillator
/*
Copyright 2010-2012 Karsten Ahnert
Copyright 2011-2013 Mario Mulansky
Copyright 2013 Pascal Germroth
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <iostream>
#include <vector>
#include <boost/numeric/odeint.hpp>
//[ rhs_function
/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;
const double gam = 0.15;
/* The rhs of x' = f(x) */
void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - gam*x[1];
}
//]
//[ rhs_class
/* The rhs of x' = f(x) defined as a class */
class harm_osc {
double m_gam;
public:
harm_osc( double gam ) : m_gam(gam) { }
void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - m_gam*x[1];
}
};
//]
//[ integrate_observer
struct push_back_state_and_time
{
std::vector< state_type >& m_states;
std::vector< double >& m_times;
push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
: m_states( states ) , m_times( times ) { }
void operator()( const state_type &x , double t )
{
m_states.push_back( x );
m_times.push_back( t );
}
};
//]
struct write_state
{
void operator()( const state_type &x ) const
{
std::cout << x[0] << "\t" << x[1] << "\n";
}
};
int main(int /* argc */ , char** /* argv */ )
{
using namespace std;
using namespace boost::numeric::odeint;
//[ state_initialization
state_type x(2);
x[0] = 1.0; // start at x=1.0, p=0.0
x[1] = 0.0;
//]
//[ integration
size_t steps = integrate( harmonic_oscillator ,
x , 0.0 , 10.0 , 0.1 );
//]
//[ integration_class
harm_osc ho(0.15);
steps = integrate( ho ,
x , 0.0 , 10.0 , 0.1 );
//]
//[ integrate_observ
vector<state_type> x_vec;
vector<double> times;
steps = integrate( harmonic_oscillator ,
x , 0.0 , 10.0 , 0.1 ,
push_back_state_and_time( x_vec , times ) );
/* output */
for( size_t i=0; i<=steps; i++ )
{
cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
}
//]
//[ define_const_stepper
runge_kutta4< state_type > stepper;
integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
//]
//[ integrate_const_loop
const double dt = 0.01;
for( double t=0.0 ; t<10.0 ; t+= dt )
stepper.do_step( harmonic_oscillator , x , t , dt );
//]
//[ define_adapt_stepper
typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
//]
//[ integrate_adapt
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
//]
{
//[integrate_adapt_full
double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
controlled_stepper_type controlled_stepper(
default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
//]
}
//[integrate_adapt_make_controlled
integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
//]
//[integrate_adapt_make_controlled_alternative
integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
//]
#ifdef BOOST_NUMERIC_ODEINT_CXX11
//[ define_const_stepper_cpp11
{
runge_kutta4< state_type > stepper;
integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) {
dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; }
, x , 0.0 , 10.0 , 0.01 );
}
//]
//[ harm_iterator_const_step]
std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
[]( std::pair< const state_type & , const double & > x ) {
cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );
//]
#endif
}
cmake_minimum_required(VERSION 2.8.3)
set(Python_ADDITIONAL_VERSIONS 2.7 2.6 2.5)
FIND_PACKAGE(PythonInterp)
FIND_PACKAGE(PythonLibs)
FIND_PACKAGE(Boost COMPONENTS python)
INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS})
LINK_LIBRARIES(${Boost_LIBRARIES} ${PYTHON_LIBRARIES})
PYTHON_ADD_MODULE(hello hello.cpp)
TARGET_LINK_LIBRARIES(hello ${Boost_LIBRARIES})
# location of the Python header files
PYTHON_VERSION = 2.7
PYTHON_INCLUDE = /usr/include/python$(PYTHON_VERSION)
# location of the Boost Python include files and library
BOOST_INC = /usr/include
BOOST_LIB = /usr/lib
# compile example
TARGET = hello
$(TARGET).so: $(TARGET).o
g++ -shared -Wl,--export-dynamic $(TARGET).o -L$(BOOST_LIB) -lboost_python -L/usr/lib/python$(PYTHON_VERSION)/config -lpython$(PYTHON_VERSION) -o $(TARGET).so
$(TARGET).o: $(TARGET).cpp
g++ -I$(PYTHON_INCLUDE) -I$(BOOST_INC) -fPIC -c $(TARGET).cpp
char const* greet()
{
return "hello, world";
}
#include <boost/python.hpp>
BOOST_PYTHON_MODULE(hello)
{
using namespace boost::python;
def("greet", greet);
}
CVODE_INC = /usr/include/sundials
CVODE_LIB = /usr/lib
cvRoberts_dns: cvRoberts_dns.o
g++ -O -o $@ $^ -L$(CVODE_LIB) -lsundials_cvode -lsundials_nvecserial
%.o: %.cpp
g++ -I$(CVODE_INC) -c $^ -o $@
clean:
rm *.o
rm cvRoberts_dns
/*
* -----------------------------------------------------------------
* $Revision: 4074 $
* $Date: 2014-04-23 14:13:52 -0700 (Wed, 23 Apr 2014) $
* -----------------------------------------------------------------
* Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
* Radu Serban @ LLNL
* -----------------------------------------------------------------
* Example problem:
*
* The following is a simple example problem, with the coding
* needed for its solution by CVODE. The problem is from
* chemical kinetics, and consists of the following three rate
* equations:
* dy1/dt = -.04*y1 + 1.e4*y2*y3
* dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
* dy3/dt = 3.e7*(y2)^2
* on the interval from t = 0.0 to t = 4.e10, with initial
* conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
* While integrating the system, we also use the rootfinding
* feature to find the points at which y1 = 1e-4 or at which
* y3 = 0.01. This program solves the problem with the BDF method,
* Newton iteration with the CVDENSE dense linear solver, and a
* user-supplied Jacobian routine.
* It uses a scalar relative tolerance and a vector absolute
* tolerance. Output is printed in decades from t = .4 to t = 4.e10.
* Run statistics (optional outputs) are printed at the end.
* -----------------------------------------------------------------
*/
#include <stdio.h>
/* Header files with a description of contents used */
#include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
#include <cvode/cvode_dense.h> /* prototype for CVDense */
#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
#include <sundials/sundials_types.h> /* definition of type realtype */
/* User-defined vector and matrix accessor macros: Ith, IJth */
/* These macros are defined in order to write code which exactly matches
the mathematical problem description given above.
Ith(v,i) references the ith component of the vector v, where i is in
the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
using the N_VIth macro in nvector.h. N_VIth numbers the components of
a vector starting from 0.
IJth(A,i,j) references the (i,j)th element of the dense matrix A, where
i and j are in the range [1..NEQ]. The IJth macro is defined using the
DENSE_ELEM macro in dense.h. DENSE_ELEM numbers rows and columns of a
dense matrix starting from 0. */
#define Ith(v,i) NV_Ith_S(v,i-1) /* Ith numbers components 1..NEQ */
#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
/* Problem Constants */
#define NEQ 3 /* number of equations */
#define Y1 RCONST(1.0) /* initial y components */
#define Y2 RCONST(0.0)
#define Y3 RCONST(0.0)
#define RTOL RCONST(1.0e-4) /* scalar relative tolerance */
#define ATOL1 RCONST(1.0e-8) /* vector absolute tolerance components */
#define ATOL2 RCONST(1.0e-14)
#define ATOL3 RCONST(1.0e-6)
#define T0 RCONST(0.0) /* initial time */
#define T1 RCONST(0.4) /* first output time */
#define TMULT RCONST(10.0) /* output time factor */
#define NOUT 12 /* number of output times */
/* Functions Called by the Solver */
static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
static int g(realtype t, N_Vector y, realtype *gout, void *user_data);
static int Jac(long int N, realtype t,
N_Vector y, N_Vector fy, DlsMat J, void *user_data,
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
/* Private functions to output results */
static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3);
static void PrintRootInfo(int root_f1, int root_f2);
/* Private function to print final statistics */
static void PrintFinalStats(void *cvode_mem);
/* Private function to check function return values */
static int check_flag(void *flagvalue, char *funcname, int opt);
/*
*-------------------------------
* Main Program
*-------------------------------
*/
int main()
{
realtype reltol, t, tout;
N_Vector y, abstol;
void *cvode_mem;
int flag, flagr, iout;
int rootsfound[2];
y = abstol = NULL;
cvode_mem = NULL;
/* Create serial vector of length NEQ for I.C. and abstol */
y = N_VNew_Serial(NEQ);
if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
abstol = N_VNew_Serial(NEQ);
if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);
/* Initialize y */
Ith(y,1) = Y1;
Ith(y,2) = Y2;
Ith(y,3) = Y3;
/* Set the scalar relative tolerance */
reltol = RTOL;
/* Set the vector absolute tolerance */
Ith(abstol,1) = ATOL1;
Ith(abstol,2) = ATOL2;
Ith(abstol,3) = ATOL3;
/* Call CVodeCreate to create the solver memory and specify the
* Backward Differentiation Formula and the use of a Newton iteration */
cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
/* Call CVodeInit to initialize the integrator memory and specify the
* user's right hand side function in y'=f(t,y), the inital time T0, and
* the initial dependent variable vector y. */
flag = CVodeInit(cvode_mem, f, T0, y);
if (check_flag(&flag, "CVodeInit", 1)) return(1);
/* Call CVodeSVtolerances to specify the scalar relative tolerance
* and vector absolute tolerances */
flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
/* Call CVodeRootInit to specify the root function g with 2 components */
flag = CVodeRootInit(cvode_mem, 2, g);
if (check_flag(&flag, "CVodeRootInit", 1)) return(1);
/* Call CVDense to specify the CVDENSE dense linear solver */
flag = CVDense(cvode_mem, NEQ);
if (check_flag(&flag, "CVDense", 1)) return(1);
/* Set the Jacobian routine to Jac (user-supplied) */
flag = CVDlsSetDenseJacFn(cvode_mem, Jac);
if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
/* In loop, call CVode, print results, and test for error.
Break out of loop when NOUT preset output times have been reached. */
printf(" \n3-species kinetics problem\n\n");
iout = 0; tout = T1;
while(1) {
flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
if (flag == CV_ROOT_RETURN) {
flagr = CVodeGetRootInfo(cvode_mem, rootsfound);
if (check_flag(&flagr, "CVodeGetRootInfo", 1)) return(1);
PrintRootInfo(rootsfound[0],rootsfound[1]);
}
if (check_flag(&flag, "CVode", 1)) break;
if (flag == CV_SUCCESS) {
iout++;
tout *= TMULT;
}
if (iout == NOUT) break;
}
/* Print some final statistics */
PrintFinalStats(cvode_mem);
/* Free y and abstol vectors */
N_VDestroy_Serial(y);
N_VDestroy_Serial(abstol);
/* Free integrator memory */
CVodeFree(&cvode_mem);
return(0);
}
/*
*-------------------------------
* Functions called by the solver
*-------------------------------
*/
/*
* f routine. Compute function f(t,y).
*/
static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
{
realtype y1, y2, y3, yd1, yd3;
y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
yd1 = Ith(ydot,1) = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
yd3 = Ith(ydot,3) = RCONST(3.0e7)*y2*y2;
Ith(ydot,2) = -yd1 - yd3;
return(0);
}
/*
* g routine. Compute functions g_i(t,y) for i = 0,1.
*/
static int g(realtype t, N_Vector y, realtype *gout, void *user_data)
{
realtype y1, y3;
y1 = Ith(y,1); y3 = Ith(y,3);
gout[0] = y1 - RCONST(0.0001);
gout[1] = y3 - RCONST(0.01);
return(0);
}
/*
* Jacobian routine. Compute J(t,y) = df/dy. *
*/
static int Jac(long int N, realtype t,
N_Vector y, N_Vector fy, DlsMat J, void *user_data,
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
realtype y1, y2, y3;
y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
IJth(J,1,1) = RCONST(-0.04);
IJth(J,1,2) = RCONST(1.0e4)*y3;
IJth(J,1,3) = RCONST(1.0e4)*y2;
IJth(J,2,1) = RCONST(0.04);
IJth(J,2,2) = RCONST(-1.0e4)*y3-RCONST(6.0e7)*y2;
IJth(J,2,3) = RCONST(-1.0e4)*y2;
IJth(J,3,2) = RCONST(6.0e7)*y2;
return(0);
}
/*
*-------------------------------
* Private helper functions
*-------------------------------
*/
static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3)
{
#if defined(SUNDIALS_EXTENDED_PRECISION)
printf("At t = %0.4Le y =%14.6Le %14.6Le %14.6Le\n", t, y1, y2, y3);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
printf("At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
#else