I wanted to solve a kind of ordinary differential equation (master equation) and I wrote the following program in c++ by help of armadillo:
#include <iostream>
#include <armadillo>
#include <iomanip>
using namespace std;
using namespace arma;
cx_mat tens( cx_mat a1,cx_mat a2,cx_mat a3,cx_mat a4,cx_mat a5,cx_mat
a6,cx_mat a7,cx_mat a8,cx_mat a9,cx_mat a10,cx_mat a11,cx_mata12,cx_mat a13,cx_mat a14,cx_mat a15,cx_mat a16,cx_mat a17,cx_mat a18,cx_mat a19,cx_mat a20,cx_mat a21)
{return kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(a1,a2),a3),a4),a5),a6),a7),a8),a9),a10),a11),a12),a13),a14),a15),a16),a17),a18),a19),a20),a21);}
cx_mat ii(2,2,fill::eye);// make a 2*2 identify cx_matrix
cx_mat ee = ii.col(0); // extract a column vector
cx_mat gg = ii.col(1);
cx_mat a1 =tens(ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a2 =tens(gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a3 =tens(gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a4 =tens(gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a5 =tens(gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a6 =tens(gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a7 =tens(gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a8 =tens(gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a9 =tens(gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a10=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a11=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a12=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a13=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg,gg);
cx_mat a14=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg,gg);
cx_mat a15=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg,gg);
cx_mat a16=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg,gg);
cx_mat a17=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg,gg);
cx_mat a18=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg,gg);
cx_mat a19=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg,gg);
cx_mat a20=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee,gg);
cx_mat a21=tens(gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,ee);
cx_mat sink=a21*a20.t();
cx_mat H0(cx_mat a){
return a*a.t();}
cx_mat H1(cx_mat a,cx_mat b){
return a*b.t()+b*a.t();}
cx_mat H00=H0(a1)+H0(a2)+H0(a3)+H0(a4)+H0(a5)+H0(a6)+H0(a7);
cx_mat H11=H1(a1,a2)+H1(a1,a3)+H1(a1,a4)+H1(a1,a5)+H1(a1,a6)+H1(a1,a7)+H1(a1,a8)+H1(a1,a9)+H1(a1,a10)+H1(a1,a11)+H1(a1,a12)+H1(a1,a13)+H1(a1,a14)+H1(a1,a15)+H1(a1,a16)+H1(a1,a17)+H1(a1,a18)+H1(a1,a19)+H1(a1,a20);
cx_mat H=H00+H11;//system Hamiltonian
cx_mat rhot(float t,cx_mat y){
return complex<double>(0, 1)*(-H*y+y*H)+0.5*(2*sink*y*sink.t()-sink.t()*sink*y-y*sink.t()*sink);}//Master equation
int rk4(cx_mat y,float dt,float tmax)//Runge kutta 4th order
{float t = 0.;cx_mat ydot1, ydot2, ydot3, ydot4;
while (t < tmax)
{
ydot1 = rhot(t, y);
ydot2 = rhot(t+0.5*dt, y+0.5*dt*ydot1);
ydot3 = rhot(t+0.5*dt, y+0.5*dt*ydot2);
ydot4 = rhot(t+dt, y+dt*ydot3);
cout<<t<< real(a21.t()*y*a21) ;
y=y+ (dt/6.0)*(ydot1 + 2.0*ydot2 + 2.0*ydot3 + ydot4);
t=t+ dt;
}
return 0;
}
int main()
{
rk4(a1*a1.t(),0.01,40.);
return 0;
}
I run this program by typying the following comment on the Ubuntu terminal:
g++ -std=c++0x psinkt.cpp -o ./psinkt.out -O3 -march=native -larmadillo
but I encountered to the following memory error:
error: arma::memory::acquire(): out of memory
terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc
Aborted (core dumped)
Generally, is there a way to solve such a problem? if yes, please tell me even the key word ! I really need to solve the problem.
Step 1: Identify where it happens.
Compile with
and debug with
Or use valgrind
Step 2: Improve your C++
I'd strongly encourage you to visit a book store and skim some C++ books and find a more readable style. Your current code is very hard to penetrate and coming from Python to C++ is going to introduce a level of hurt like this - C++ is a language full of stupidly sharp edges and loaded hand cannons; it's C++, not you :)
For example, that horrible kron thing you're doing... Rethink what it is you're trying to do, and things like this:
Here's the bad news: You're going to have to write more code in C++ and use less direct paths compared to how you write in Python. When you are more familiar with the language, its idioms, etc, the paths seem more obvious and natural, but coming from Python, stuff is often going to seem backwards.
You basically have to provide a lot more detail to the system, specificity is the cost of performance.
For example, you might want to consider putting your "a"s into a
std::vector
so that you can say things like:You could refactor your kron nesting like this:
const cx_mat&
says to pass the value by reference, I don't know whethercx_mat
is a trivial object or whether passing it by value as you are is expensive (pass by value requires a deep copy).Or you could write it like this:
Or you could use C++11 templates:
But that's my least favorite of the options. Note: the elipsis (...) are not me leaving stuff out, that's actual C++11 variadic template syntax: see http://ideone.com/VgfmVB
I don't know that the above code examples solve your problem or which is preferable in your case, but I'm hoping that between the two I've equipped you to make progress.