I have a simple 2D (row, column) matrix which I currently reorder according to the algorithm below, using another array as final container to swap items.
The problem is that I need to save memory (the code is running on a very low end device), and thus I need to figure a way to reorder the array in-place.
The algorithm is as follows:
for (int iRHS = 0; iRHS < NUM_VARS; iRHS )
for (int iRow = 0; iRow < _numEquations; iRow ) {
coef[iRHS][iRow] = _matrixCoef(iRow, iRHS);
}
Note: coef is a pointer to double accessed via subscript, _matrixCoef is a matrix helper class and uses a vector of double accessed by operator(row,col). Here I want to eliminate coef, so that all values are reordered in _matrixCoef in-place instead.
Edit: NUM_VARS is a define set to 2.
Is this possible in-place after all?
Edit 2:
Here is the matrix class which is accessed above via operator overload (row, col):
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;
};
CodePudding user response:
After you posted code, I will suggest another solution, that's rather simple and quick to implement.
In your current Matrix class:
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;
};
You may want to modify other member functions, depending on your application.
CodePudding user response:
Well, assuming your matrix is a square matrix, that is NUM_VARS == _numEquations, it is possible. Otherwise the size of the resulting matrix will not allow in-place transposition. In that case, a solution would be to modify the downstream computations to swap the row/column indices while doing the operations.
But if it's square, you can try this:
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]);
}
CodePudding user response:
If you want to explicitly re-order the data array in-place:
Consider this post that details in-place re-ordering of an array given a permutation. I have modified it only slightly. Algorithm to apply permutation in constant memory space
The expression for the permutation can be deduced as follows:
For index k = i j * rows
, j = k / rows
and i = k - j * rows
(integer division). The index of the transposed matrix is k_transpose = j i * cols
. Now for the code:
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]);
}
}
Example with a symbolic matrix:
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";
Output:
a,b,c,d,e,f,g,h,
a,e,b,f,c,g,d,h,
Of course, you will also have to update your struct with the correct rows and columns. The asymptotic computational complexity is not guaranteed to be linear, but there are no auxiliary data structures that consume your precious memory!