我想求解一种常微分方程(主方程),并在犰狳的帮助下用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/