我想求解一种常微分方程(主方程),并在犰狳的帮助下用C ++编写了以下程序:

#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;
  }


我通过在Ubuntu终端上键入以下注释来运行该程序:

 g++ -std=c++0x psinkt.cpp -o ./psinkt.out -O3 -march=native -larmadillo


但是我遇到了以下内存错误:

error: arma::memory::acquire(): out of memory

terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
Aborted (core dumped)


通常,有没有办法解决这个问题?如果是,请甚至告诉我关键词!
我真的需要解决这个问题。

最佳答案

步骤1:找出发生的地方。

编译

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


和调试

$ gdb ./psinkt.out
gdb> run


或使用valgrind

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


步骤2:改善您的C ++

我强烈建议您访问书店并浏览一些C ++书籍,并找到更具可读性的样式。您当前的代码很难渗透,并且从Python到C ++都会带来类似的伤害。它是C ++,不是你:)

例如,您正在做的可怕的kron事情...重新考虑您想做的事情,像这样的事情:

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


这是一个坏消息:与用Python编写相比,您将不得不用C ++编写更多代码并且使用更少的直接路径。当您对语言,它的惯用法等更加熟悉时,这些路径似乎更加明显和自然,但是来自Python时,东西通常看起来会倒退。

您基本上必须为系统提供更多细节,特异性就是性能成本。

例如,您可能需要考虑将您的“ a”放入std::vector中,这样您可以说出以下内容:

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


您可以像这样重构您的kron嵌套:

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&表示要通过引用传递值,我不知道cx_mat是琐碎的对象还是因为按价传递值比较昂贵(按值传递需要深层复制)。

或者您可以这样写:

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


或者,您可以使用C ++ 11模板:

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


但这是我最不喜欢的选择。注意:省略号(...)不是我遗漏的东西,这是实际的C ++ 11可变参数模板语法:请参见http://ideone.com/VgfmVB

我不知道上面的代码示例可以解决您的问题,或者哪种情况更适合您,但我希望在这两者之间,我已经为您提供了帮助。

关于c++ - C++中的内存错误(armadillo),我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/30762396/

10-12 17:52