我必须生成一个复数多维数组 N (j, k, m, l)
,大约为 10*100*100*1000。
我想执行以下操作来计算这个 N 并返回。
for j<10 ...
{
for k<100 ...
{
......
some matrix multiplication to generate a 2D complex valued matrix n(100*1000)
......
N(j,k,:,:)= n
}
}
我的问题:
如何高效地实现
N(j,k,:,:)= n
。对于目前的问题规模,我应该从头开始编码还是使用一些现有的库?
最佳答案
你说的是 10*100*100*1000 = 100,000,000 个复数,如果 2 个 float
可能每个 8 个字节,或者 2 个 double
可能每个 16 个字节,所以大约 800 兆字节或 1.6 GB。在普通台式 PC 的容量范围内,这是一个好的开始。
有效分配的主要内容是确保内存布局使得分配处理连续内存。您可以编写几个类来提供一个漂亮的界面 - 比如说 Matrix_2D
然后是 Matrix_4D
,例如:
template <typename T>
class Matrix_4D
{
public:
Matrix_4D(size_t j, size_t k, size_t l, size_t m)
: j_(j), k_(k), l_(l), m_(m), data_(new T[j * k * l * m]),
klm_(k * l * m), lm_(l * m),
{ /* optionally, initialise elements */ }
~Matrix_4D() { delete data_; }
T& operator()(size_t j, size_t k, size_t l, size_t m)
{
return data_[j * klm_ + k * lm_ + l * m_ + m];
}
const T& operator()(size_t j, size_t k, size_t l, size_t m) const
{
return data_[j * klm_ + k * lm_ + l * m_ + m];
}
void set(size_t l, size_t m, const Matrix_2D& m2)
{
if (m2.j_ != l_ || m2.k_ != m_)
throw std::runtime_error("mismatched dimensions");
std::copy(m2.data_[0], m2.data_[lm_], (*this)(l, m, 0, 0));
}
private:
size_t j_, k_, l_, m_;
size_t klm_, lm_; // needed so often -> save
T* data_;
};
矩阵类应该是 friend
,这样它们就可以相互提取数据。如果你想更漂亮,你实际上可以提供一个代理对象——将以下内容添加到 Matrix_4D
struct Proxy_2D
{
Proxy_2D(Matrix_4D& m4, size_t l, size_t m) : m4_(m4), l_(l), m_(m) { }
Proxy_2D& operator=(const Matrix2D& m2)
{
m4_.set(l_, m_, m2);
return *this;
}
Matrix_4D& m4_;
size_t l_, m_;
};
Proxy_2D operator()(size_t l, size_t m) { return Proxy_2D(*this, l, m); }
然后你可以这样做:
Matrix_4D m4(10, 20, 30, 40);
Matrix_2D m2(30, 40);
... set stuff in m2 ...
m4(2, 4) = m2;
编辑:对于您评论中的代码 - m2= m2 * transpose(m2)
- 如果您想追求这种自己动手的实现来学习 C++ 而不是获取现有的使用模板表达式等高性能技术的高效库(太复杂了,无法进入此处),然后在 Matrix_2D
中:
Matrix_2D transpose() const
{
Matrix_2D result(m_, l_);
for (size_t l = 0; l < l_; ++l)
for (size_t m = 0; m < m_; ++m)
result(m, l)= (*this)(l, m);
return result;
}
Matrix_2D& operator+=(const Matrix_2D& rhs)
{
for (size_t l = 0; l < l_; ++l)
for (size_t m = 0; m < m_; ++m)
(*this)(l, m) += rhs(l, m);
return *this;
}
Matrix_2D operator+(const Matrix_2D& rhs) const
{
Matrix_2D result(*this); // copy *this
return result += rhs;
}
有趣的是,您还可以在不复制数据的情况下将转置作为矩阵的一种动态透视,但是您需要确保底层矩阵对象的生命周期跨越转置对象的使用:
template <typename T>
class Transpose_2D
{
public:
Transpose_2D(Matrix_2D<T>& m) : m_(m) { }
T& operator()(size_t l, size_t m) { return m_(m, l); }
const T& operator()(size_t l, size_t m) const { return m_(m, l); }
private:
Matrix_2D<T>& m_;
};
相应地更改 Matrix_2D 加法函数签名允许使用它,例如:
template <typename U>
Matrix_2D& operator+=(const U& rhs)
...
然后你可以这样做:
m2 += Transpose_2D(m2);
而且它会相当高效。
关于c++ - C++ 中的多维复数值数组,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23076925/