就地重新排序简单的二维矩阵

Question

我有一个简单的 2D（行、列）矩阵，我目前根据下面的算法重新排序，使用另一个数组作为最终容器来交换项目。

问题是我需要节省内存（代码在非常低端的设备上运行），因此我需要想办法就地重新排序数组。

算法如下：

for (int iRHS = 0; iRHS < NUM_VARS; iRHS++)
    for (int iRow = 0; iRow < _numEquations; iRow++) {
        coef[iRHS][iRow] = _matrixCoef(iRow, iRHS); 
    }

注意：coef 是指向通过下标进行双重访问的指针，_matrixCoef 是一个矩阵辅助类，并使用由 operator(row,col) 访问的双重向量。在这里我想消除coef，以便所有值都在_matrixCoef中重新排序。

编辑：NUM_VARS 是定义设置为 2。

毕竟这可能就地吗？

编辑 2：

这是上面通过运算符重载（行，列）访问的矩阵类：

struct Matrix
{
    /// Creates a matrix with zero rows and columns.
    Matrix() = default;
    /// Creates a matrix with \a rows rows and \a col columns
    /// Its elements are initialized to 0.
    Matrix(int rows, int cols) : n_rows(rows), n_cols(cols), v(rows * cols, 0.) {}
    /// Returns the number or rows of the matrix
    inline int getNumRows() const { return n_rows; }
    /// Returns the number or columns of the matrix.
    inline int getNumCols() const { return n_cols; }
    /// Returns the reference to the element at the position \a row, \a col.
    inline double & operator()(int row, int col) { return v[row + col * n_rows]; }
    /// Returns the element at the position \a row, \a col by value.
    inline double operator()(int row, int col) const { return  v[row + col * n_rows]; }
    /// Returns the values of the matrix in column major order.
    double const * data() const { return v.data(); }
    /// Returns the values of the matrix in column major order.
    double * data() { return v.data(); }
    /// Initialize the matrix with given size. All values are set to zero.
    void initialize(int iRows, int iCols)
    {
        n_rows = iRows;
        n_cols = iCols;
        v.clear();
        v.resize(iRows * iCols);
    }
    
    void resize(int iRows, int iCols)
    {
        n_rows = iRows;
        n_cols = iCols;
        v.resize(iRows * iCols);
    }
private:
    int n_rows = 0;
    int n_cols = 0;
    std::vector<double> v;
};

Answer 1

在您发布代码后，我会建议另一种解决方案，该解决方案相当简单且易于实施。

在您当前的 Matrix 类中：

struct Matrix
{
    // ...

    // add this:
       void transpose()
       {
           is_transposed = !is_transposed;
       }
    // ...

    // modify these:

    /// Returns the number or rows of the matrix
    inline int getNumRows() const { return (is_transposed) ? n_cols : n_rows; }
    /// Returns the number or columns of the matrix.
    inline int getNumCols() const { return (is_transposed) ? n_rows : n_cols; }
    /// Returns the reference to the element at the position \a row, \a col.
    inline double & operator()(int row, int col) 
    {
        if (is_transposed) 
            return v[col + row * n_rows]; 
        return v[row + col * n_rows]; 
    }
    /// Returns the element at the position \a row, \a col by value.
    inline double operator()(int row, int col) const 
    { 
        if (is_transposed) 
            return  v[col + row * n_rows]; 
        return  v[row + col * n_rows]; 
    }

private:
    // ...

    // add this:
    bool is_transposed = false;
};

您可能想要修改其他成员函数，具体取决于您的应用程序。

Answer 2

好吧，假设您的矩阵是方阵，即 NUM_VARS == _numEquations，这是可能的。否则，生成的矩阵的大小将不允许就地转置。在这种情况下，解决方案是修改下游计算以在执行操作时交换行/列索引。

但如果是方形的，你可以试试这个：

for (size_t r = 0; r < NUM_VARS; ++r)
    for (size_t c = 0; c < NUM_VARS; ++c)
    {
         if (&_matrixCoef[r][c] < &_matrixCoef[c][r])
             std::swap(_matrixCoef[r][c], _matrixCoef[c][r]);
    }

Answer 3

如果你想明确地重新排序数据数组：

考虑一下这篇文章，它详细介绍了给定排列的数组的就地重新排序。我只是稍微修改了一下。 Algorithm to apply permutation in constant memory space

置换表达式可推导出如下：

对于索引k = i + j * rows、j = k / rows 和i = k - j * rows（整数除法）。转置矩阵的索引是k_transpose = j + i * cols。现在看代码：

int perm(int k, int rows, int cols)
{
    int j = k / rows;
    int i = k - j * rows;
    return j + i * cols;
}
template<typename Scalar>
void transpose_colmajor(Scalar* A, int rows, int cols)
{
//tiny optimization: for (int i = 1; i < rows * cols - 1; i++)
    for (int i = 0; i < rows * cols; i++) {
        int ind = perm(i, rows, cols);
        while (ind > i)
            ind = perm(ind, rows, cols);
        std::swap(A[i], A[ind]);
    }
}

符号矩阵示例：

    int rows = 4;
    int cols = 2;
    char entry = 'a';
    std::vector<char> symbolic_matrix;
    for (int col = 0; col < cols; col++)
    {
        for (int row = 0; row < rows; row++)
        {
            symbolic_matrix.push_back(entry);
            entry++;
        }
    }
    for (char coefficient : symbolic_matrix)
        std::cout << coefficient << ",";
    std::cout << "\n\n";

    transpose_colmajor(symbolic_matrix.data(), rows, cols);

    for (char coefficient : symbolic_matrix)
        std::cout << coefficient << ",";
    std::cout << "\n\n";

输出：

a,b,c,d,e,f,g,h,

a,e,b,f,c,g,d,h,

当然，您还必须使用正确的行和列更新您的结构。渐近计算复杂度不保证是线性的，但是没有辅助数据结构会消耗你宝贵的内存！

编辑： 我不得不承认，我从这个问题中学到了很多东西。我想知道为什么转置方阵如此直观，而矩形则不然。它与排列的循环有关。这个算法https://www.geeksforgeeks.org/minimum-number-swaps-required-sort-array/ 计算给定排列所需的最小交换。如果您插入上面的置换表达式，您可以看到在rows != cols 时，最小交换次数急剧增加。例如，4 x 4 矩阵需要6 交换，而4 x 3 矩阵需要8，这对我来说非常违反直觉！我发布的算法需要rows * cols - 2 交换，这不是最优的。然而，基于矩形矩阵的随机实验，最小互换之间的差距似乎非常小，并且对于较大的矩阵不太重要。绝对值得专门研究方阵的转置例程，因为最小交换总是n * (n - 1)/2。