在基本的C++中，我发现这里的问题是C++会分配一个新的cx_mat对象来存储evol​​utionMatrix*stateMatrix，然后用operator=()将这个新对象复制到stateMatrix。

我认为你是对的，它正在创造临时工，这太慢了，但我认为它这样做的原因是错误的。

Armadillo 与任何优秀的 C++ 线性代数库一样，使用表达式模板来实现表达式的延迟求值。当您写下类似 的矩阵乘积A*B时，不会创建临时对象，而是 Armadillo 创建一个临时对象 ( x)，该对象保留对 A和的引用B，然后，给定类似 的表达式C = x，计算将结果直接存储在 中的矩阵乘积C，而不创建任何临时工。

它还使用这种优化来处理表达式，如A*B*C*D，其中，根据矩阵大小，某些乘法顺序比其他顺序更有效。

这不是犰狳的实施问题吗？

如果犰狳没有执行此优化，那将是犰狳中的一个错误，应报告给开发人员。

但是，在您的情况下，还有另一个更重要的问题。A=B*C在类似storage of的表达式中，A如果A没有别名B或C. 在您的情况下， A = A*B向输出矩阵写入任何内容也会修改输入矩阵之一。

即使给定您建议的功能

multiply(const cx_mat&, const cx_mat&, cx_mat&)

该函数究竟对表达式有什么帮助multiply(A, B, A)？对于该功能的大多数普通实现，这将导致错误。它需要自己使用一些临时存储，以确保其输入数据没有损坏。你的建议几乎是犰狳已经如何实现矩阵乘法，但我认为它可能需要注意避免像 multiply(A, B, A)分配临时这样的情况。

犰狳为什么不做这种优化的最可能的解释是这样做是不正确的。

最后，有一个更简单的方法来做你想做的事，像这样：

cx_mat *A, *Atemp, B;
for (;;) {
  *Atemp = (*A)*B;
  swap(A, Atemp);
}

这与

cx_mat A, B;
for (;;) {
  A = A*B;
}

但它分配一个临时矩阵，而不是每次迭代一个临时矩阵。

@BillGreene 指出“返回值优化”是解决基本问题的一种方法，但这实际上只帮助了其中的一半。假设您有这种形式的代码：

struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); };

ExpensiveObject operator+ (ExpensiveObject &obj1,
                           ExpensiveObject &obj2)
{
   ExpensiveObject tmp;
   ...compute tmp based on obj1 and obj2...
   return tmp;
}

void f() {
  ExpensiveObject o1, o2, o3;
  ...initialize o1, o2...;
  o3 = o1 + o2;
}

一个天真的编译器会

1. 创建一个槽来存储加号操作的结果（一个临时的），
2. 呼叫操作员+，
3. 在 operator+ 中创建 'tmp' 对象（第二个临时对象），
4. 计算 tmp，
5. 将 tmp 复制到结果槽中，
6. 破坏tmp，
7. 将结果对象复制到 o3
8. 销毁结果对象

返回值优化只能统一 'tmp' 对象和 'result' 槽，但不能去掉对副本的需要。因此，您仍然需要创建临时文件、复制操作和销毁临时文件。

解决这个问题的唯一方法是 operator+ 不返回对象，而是某个中间类的对象，当分配给 an 时ExpensiveObject，它会就地执行加法和复制操作。这是表达式模板库中典型使用的方法。

Stackoverflow ( https://stackoverflow.com/ ) 可能是这个问题的更好讨论论坛。但是，这是一个简短的答案。

我怀疑 C++ 编译器是否像你上面描述的那样为这个表达式生成代码。所有现代 C++ 编译器都实现了一种称为“返回值优化”（http://en.wikipedia.org/wiki/Return_value_optimization）的优化。通过返回值优化，结果evolutionMatrix*stateMatrix直接存储在stateMatrix; 没有复制。

在这个话题上显然存在相当大的混乱，这也是我建议 Stackoverflow 可能是一个更好的论坛的原因之一。那里有很多 C++“语言律师”，而我们这里的大多数人宁愿把时间花在 CSE 上。;-)

我根据 Bangerth 教授的帖子创建了以下简单示例：

#ifndef NDEBUG
#include <iostream>

using namespace std;
#endif

class ExpensiveObject  {
public:
  ExpensiveObject () {
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor called." << endl;
#endif
    v = 0;
  }
  ExpensiveObject (int i) { 
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor(int) called." << endl;
#endif
    v = i; 
  }
  ExpensiveObject (const ExpensiveObject  &a) {
    v = a.v;
#ifndef NDEBUG
    cout << "ExpensiveObject  copy constructor called." << endl;
#endif
  }
  ~ExpensiveObject() {
#ifndef NDEBUG
    cout << "ExpensiveObject  destructor called." << endl;
#endif
  }
  ExpensiveObject  operator=(const ExpensiveObject  &a) {
#ifndef NDEBUG
    cout << "ExpensiveObject  assignment operator called." << endl;
#endif
    if (this != &a) {
      return ExpensiveObject (a);
    }
  }
  void print() const {
#ifndef NDEBUG
    cout << "v=" << v << endl;
#endif
  }
  int getV() const {
    return v;
  }
private:
  int v;
};

ExpensiveObject  operator+(const ExpensiveObject  &a1, const ExpensiveObject  &a2) {
#ifndef NDEBUG
  cout << "ExpensiveObject  operator+ called." << endl;
#endif
  return ExpensiveObject (a1.getV() + a2.getV());
}

int main()
{
  ExpensiveObject  a(2), b(3);
  ExpensiveObject  c = a + b;
#ifndef NDEBUG
  c.print();
#endif
}

它看起来比实际上更复杂，因为我想在优化模式下编译时完全删除所有用于打印输出的代码。当我运行使用调试选项编译的版本时，我得到以下输出：

ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
v=5
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.

首先要注意的是没有构建临时变量——只有 a、b 和 c。默认构造函数和赋值运算符永远不会被调用，因为在此示例中不需要它们。

Bangerth 教授提到了表达式模板。事实上，这种优化技术对于在矩阵类库中获得良好的性能非常重要。但仅当对象表达式比简单的 a + b 更复杂时才重要。例如，如果我的测试是：

  ExpensiveObject  a(2), b(3), c(9);
  ExpensiveObject  d = a + b + c;

我会得到以下输出：

ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  destructor called.
 v=14
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.

这种情况显示了临时构造的不良构造（对构造函数的 5 次调用和对运算符 + 的两次调用）。正确使用表达式模板（这个话题远远超出了本论坛的范围）可以暂时防止这种情况发生。（对于积极性高的人，可以在http://www.amazon.com/C-Templates-The-Complete-Guide/dp/0201734842的第 18 章中找到关于表达式模板的特别易读的讨论）。

最后，编译器实际在做什么的真正“证明”来自检查编译器输出的汇编代码。对于第一个示例，在优化模式下编译时，这段代码非常简单。所有函数调用都被优化掉了，汇编代码基本上将 2 加载到一个寄存器中，将 3 加载到第二个寄存器中，然后添加它们。

我想反驳你的前提，即这是“非常、非常低效”。性能方面，它实际上是无关紧要的，因为矩阵乘法是¹ $\mathcal{O}(n^{2.8})$，而复制这些矩阵之一只是$\mathcal{O}(n^2)$. 因此对于$n$与您一样大，乘法可能会完全超过复制时间。

也就是说，除非你在复制中引入了一个巨大的常数——这实际上并不是那么牵强，因为复制的版本在另一个方面要贵得多：它需要更多的内存。因此，如果您最终不得不在硬盘之间进行交换，那么复制可能确实会成为瓶颈。但是，即使您自己不复制任何内容，强并行化算法也可能会复制一些自己的副本。确实，确保在每个步骤中不会使用太多内存的唯一方法是将乘法拆分为 的列stateMatrix，因此一次只进行小的乘法运算。例如，您可以定义

class HChunkMatrix // optimised for destructive left-multiplication updates
{
  std::vector<arma::cx_mat> colChunks; // e.g. for an m×n matrix,
                                      //  use √n chunks, each an m×√n matrix
 public:
  ...

  HChunkMatrix& operator *= (const arma::cx_mat& lhMult) {
    for (&colChunk: colChunks) {
      colChunk = lhMult * colChunk;
    }
    return *this;
  }
}

您还应该首先考虑是否需要将其演变stateMatrix为一个。如果您基本上只是想要n状态 kets 的独立时间演化，那么您不妨一个一个地演化它们，这样内存成本要低得多。特别是 if evolutionMatrixis sparse，您绝对应该检查一下！因为这基本上只是一个哈密顿量，不是吗？哈密​​顿量通常是稀疏的或近似稀疏的。

¹_{或$\mathcal{O}(n^{2.38})$，但那个实际上并没有用。}