Memory error in c++ (armadillo)

Question

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.

You solve the problem by debugging your code. Writing programs isn't just about writing code, running code, and hoping it works. Debugging is part and parcel of writing programs. — PaulMcKenzie, Jun 10 '15 at 16:50
To be honest, I didn't understand what you mean at all. Would you please clarify it more? — Haj Nasser, Jun 10 '15 at 16:52
You are supposed to debug your program first if there is an error. Did you expect to write programs and never need to debug them, i.e. they work every time? No programmer writes programs that work the first time, unless they are lucky or the program is very simple. Also, one of the reasons for having questions closed in SO is the lack of effort of the poster in debugging their code, — PaulMcKenzie, Jun 10 '15 at 16:53
Of course, this program is for 21 particles but it runs well for 3 particles. On the oher hand, I debuged it for 3 particles before and I just expanded it to 21 particles and then the problem apeared. — Haj Nasser, Jun 10 '15 at 16:58
Well use the debugger when the input is 21. What's the problem in doing that? And what is this? `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);}` You need to do this differently, as one reason may be that you're blowing up the stack with such code. — PaulMcKenzie, Jun 10 '15 at 16:58
Of course, I asked this question before in http://stackoverflow.com/questions/30528648/tensor-product-of-n-matrices-in-c-and-armadillo but I didn't recieve any responses. — Haj Nasser, Jun 10 '15 at 17:04
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/80207/discussion-between-user3109979-and-paulmckenzie). — Haj Nasser, Jun 10 '15 at 17:14
You are using global variables and initializing them outside of a function. These are big problems, as the order of initialization of global variables is undefined in C++. The statement `cx_mat ee = ii.col(0);` may not work properly, as `ii` may not be initialized (ie. full of garbage). Similarly for `cx_mat sink=a21*a20.t();`. Rewrite your code so you don't have global variables. — mtall, Jun 11 '15 at 02:59
the main problem is sink=a21*a20.t(), because when I delete all the lines after this comment and this comment, the code works well, but when I only add this lin to the code, I recieve the memory error. — Haj Nasser, Jun 11 '15 at 07:09

kfsone · Accepted Answer · 2015-06-10T17:49:41.363

Step 1: Identify where it happens.

Compile with

$ g++ -std=c++0x -Wall -O0 -g3 psinkt.cpp -o ./psinkt.out

and debug with

$ gdb ./psinkt.out
gdb> run

Or use valgrind

$ yum install valgrind
# or
$ apt-get install valgrind
$ valgrind --tool=memcheck ./psinkt.out

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:

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);

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:

cx_mat H11;
for (auto it = as.begin() + 1; it != as.end() - 1; ++it) {
    H11 += H1(a1, *it); // (a1,a2) ... (a1, a20)
}

You could refactor your kron nesting like this:

cx_mat tens(const cx_mat& background, size_t posn, const cx_mat& foreground, size_t width)
{
    std::vector<cx_mat*> mats;
    mats.resize(width);
    std::fill(mats.begin(), mats.end(), &background);
    mats[posn] = &foreground;
    cx_mat accum = *mats[0];
    for (size_t i = 1; i < width; ++i) {
        accum = kron(accum, *mats[i]);
    }
    return accum;
}

cx_mat a1 = tens(gg, 0, ee, 21);
cx_mat a2 = tens(gg, 1, ee, 21);
cx_mat a3 = tens(gg, 2, ee, 21);
...

const cx_mat& says to pass the value by reference, I don't know whether cx_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:

void tens(cx_mat& into, cx_mat& background, size_t posn, cx_mat& foreground, size_t width)
{
    std::vector<cx_mat*> mats;
    mats.resize(width);
    std::fill(mats.begin(), mats.end(), &background);
    mats[posn] = &foreground;
    into = *mats[0];
    for (size_t i = 1; i < width; ++i) {
        accum = kron(accum, *mats[i]);
    }
}

enum { Width = 21 };
std::vector<mat> amats;
amats.reserve(Width);
for (size_t i = 0; i < Width; ++i) {
    tens(amats[i], gg, i, ee, Width);
}

Or you could use C++11 templates:

cx_mat tens(const cx_mat& lhs, const cx_mat& rhs)
{
    return kron(lhs, rhs);
}

template<typename Args...>
cx_mat tens(const cx_mat& lhs, const cx_mat& rhs, Args&&... rest)
{
    return tens(kron(lhs, rhs), std::forward<Args>(rest)...);
}

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);

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.

Several of my friends who have recently moved from Python to C++ have found that - contrary to any of our expectations - MS Visual Studio (one of the free versions, e.g. VS 2015) has helped them hugely. There are even some really nice perf/profiling tools for Python available for MSVC (https://pytools.codeplex.com/). Alternatively, consider using such as Eclipse or Code::Blocks. Learning C++ from the command line is hardcore and - frankly - the most painful route. — kfsone, Jun 10 '15 at 17:47
Thank you very much dear kfsone. Yes, I have just moved from Python to c++ and I am really strange in c++. Would you please let me know name of a good c++11 book? — Haj Nasser, Jun 10 '15 at 19:25
There's http://stackoverflow.com/questions/388242/the-definitive-c-book-guide-and-list and more appropriately http://stackoverflow.com/questions/9822291/c-c-for-python-programmer but the take away is - it varies. I encourage you to see those threads to get some suggestions then find a book store or library so you can spend a while with the books that sound good. These books are often expensive and nothing is worse than getting 5 pages in and realizing the book is going to use examples you don't relate to in any way and that it is basically a heavy sleeping pill for you. — kfsone, Jun 10 '15 at 19:44
Excuse me kfsone. I did all you said about compiling and debuging, but I received tha same error (memory error). — Haj Nasser, Jun 10 '15 at 20:12
In gdb you can use the `bt` command to get a `back trace` from the current point of execution, which should tell you *where* the program is crashing and help you debug. You may want to consult http://www.cprogramming.com/gdbtutorial.html or similar for more help using gdb. — kfsone, Jun 10 '15 at 20:16

Memory error in c++ (armadillo)

1 Answers1