How to work on a sub-matrix in a matrix by pointer?
I have a matrix of size n. Take an example:
My recursive function does the processing on the elements that lie in the border of the matrix. Now I want to call it (the recursive call) on the inner square matrix:
This is the prototype of my recursive function:
void rotate(int** mat, size_t n);
I know that a 2D array is an array within an array. I know that *(mat+1) + 1) will give the memory address that should be the base address of my new matrix. This is what I tried:
rotate((int **)(*(mat+1) + 1), n-2)
But it does not work, and I get a segfault when I try to access it with [][].
Solution 1:
You cannot dereference mat+1 and reinterpret that as a pointer to a whole matrix. Instead provide the offsets as arguments to your function (I assume n-by-n square matrices):
void rotate(int** mat, size_t i, size_t j, size_t n) {
// assuming row-wise storage
int *row0 = mat[j]; // assumes j < n
int *row1 = mat[j + 1]; // assumes j + 1 < n
// access row0[i..] and row1[i..]
}
If you had continuous storage for your matrix, you could do the following instead:
rotate(int* mat, size_t i, size_t j, size_t n) {
int atIJ = mat[j * n + i]; // assuming row-wise storage
// ...
}
Solution 2:
I am not sure of your application, but I wonder if using #define for your matrix size would help....
#define X_SIZE 4
#define Y_SIZE 4
or even
#define N_SIZE 4
... because then you can use X_SIZE and Y_SIZE (OR N_SIZE) in your function without having to pass them explicitly.
in main you might put
int matrix[X_SIZE * Y_SIZE];
or
int matrix2[N_SIZE * N_SIZE];
then you can call the ith row and jth column element with
*(pmatrix + X_SIZE*j + i)
or
matrix[X_SIZE*j + i]
or
*(pmatrix2 + N_SIZE*j + i)
or
matrix2[N_SIZE*j + i]
where pmatrix and pmatrix2 are pointers to matrix and matrix2.
I am pretty sure there is no clever trick to be able to easily pass the inner square 2x2 matrix to a function, unless you were to copy the elements from the centre of your matrix into a new matrix and then copy back the result afterwards.
Solution 3:
This is not an answer to the stated question, but it is an answer to the underlying problem: management of matrices and views to matrices with minimal effort.
This will garner downvotes, but it has been so useful in solving the underlying problems whenever the type of question the OP poses has been asked, I find it is worth showing this alternative approach here. It is not interesting for small, fixed-size matrices, as the features only show their benefits when the sizes are larger or vary.
I use the following two structures to describe matrices. I shall leave out memory pool support (which allows one to manage a set of matrices as a pool, releasing them all at once, without having to manage each matrix separately) and everything related to multithreaded operation and thread-safety, for simplicity.
The code might contain typos; if you notice any, please leave a comment, and I'll fix them.
typedef int data_t; /* Matrix element data type */
struct owner {
long refcount; /* Number of referenced to this data */
size_t size; /* Number of elements in data[] */
data_t data[]; /* C99 flexible array member */
};
typedef struct {
int rows; /* Number of rows in this matrix */
int cols; /* Number of columns in this matrix */
long rowstride;
long colstride;
data_t *origin; /* Pointer to element at row 0, col 0 */
struct owner *owner; /* Owner structure origin points to */
} matrix;
#define MATRIX_INIT { 0, 0, 0L, 0L, NULL, NULL }
Matrix m element at row r, column c, is m.origin[r * m.rowstride + c * m.colstride], assuming 0 <= r && r < m.rows and 0 <= c < m.cols.
Matrices are typically declared as local variables, not as pointers. You do need to remember to free each individual matrix after you no longer need it. (The pool mechanism I omitted lets you avoid this, as all matrices in a pool are freed at once.)
Each matrix refers to exactly one owner structure. Owner structures record the number of references (the number of matrices referring to data in that structure), and are released when the reference count drops to zero:
void matrix_free(matrix *const m)
{
if (m != NULL) {
if (m->owner != NULL && --(m->owner.refcount) < 1L) {
m->owner.size = 0;
free(m->owner);
}
m->rows = 0;
m->cols = 0;
m->rowstride = 0L;
m->colstride = 0L;
m->origin = NULL;
m->owner = NULL;
}
}
Whenever a new matrix is created, the corresponding owner structure is created:
int matrix_new(matrix *const m, const int rows, const int cols)
{
const size_t size = (size_t)rows * (size_t)cols;
struct owner *o;
if (m == NULL)
return errno = EINVAL;
m->rows = 0;
m->cols = 0;
m->rowstride = 0L;
m->colstride = 0L;
m->origin = NULL;
m->owner = NULL;
if (rows < 1 || cols < 1)
return errno = EINVAL;
o = malloc(sizeof (struct owner) + size * sizeof (data_t));
if (o == NULL) {
return errno = ENOMEM;
o->refcount = 1L;
o->size = size;
m->rows = rows;
m->cols = cols;
m->origin = o->data;
m->owner = o;
#if DEFAULT_COLUMN_MAJOR > 0
/* Default to column-major element order */
m->rowstride = 1L;
m->colstride = (long)rows;
#else
/* Default to row-major element order */
m->rowstride = (long)cols;
m->colstride = 1L;
#endif
return m;
}
Note that the above does not initialize the matrix elements to any value, so they initially contain garbage.
Matrix transpose is a trivial, fast operation:
void matrix_transpose(matrix *const m)
{
if (m->rows > 0 && m->cols > 0) {
const int rows = m->rows;
const int cols = m->cols;
const long rowstride = m->rowstride;
const long colstride = m->colstride;
m->rows = cols;
m->cols = rows;
m->rowstride = colstride;
m->colstride = rowstride;
}
}
Similarly, you can rotate and mirror matrices, just remember to modify the origin member in those cases, too.
The interesting and useful cases are being able to create "views" into other matrices. The data referenced is exactly the same -- modifying one is immediately visible in the other(s); this is true aliasing --, no memory copying is needed. Unlike in most libraries (such as GSL, GNU Scientific Library), these "views" are perfectly ordinary matrices themselves. Here are some examples:
int matrix_submatrix_from(matrix *const m, const matrix *const src,
const int firstrow, const int firstcol,
const int rows, const int cols)
{
if (m == NULL || m == src)
return errno = EINVAL;
m->rows = 0;
m->cols = 0;
m->rowstride = 0L;
m->colstride = 0L;
m->origin = NULL;
m->owner = NULL;
if (firstrow + rows > src->rows ||
firstcol + cols > src->cols)
return errno = EINVAL;
if (src == NULL || src->owner == NULL)
return errno = EINVAL;
if (src->owner.refcount < 1L || src->owner.size == 0)
return errno = EINVAL;
else {
++(src->owner.refcount);
m->owner = src->owner;
}
m->origin = src->origin + src->rowstride * firstrow
+ src->colstride * firstcol;
m->rows = rows;
m->cols = cols;
m->rowstride = src->rowstride;
m->colstride = src->colstride;
return 0;
}
int matrix_transposed_from(matrix *const m, const matrix *const src)
{
if (m == NULL || m == src)
return errno = EINVAL;
m->rows = 0;
m->cols = 0;
m->rowstride = 0L;
m->colstride = 0L;
m->origin = NULL;
m->owner = NULL;
if (src == NULL || src->owner == NULL)
return errno = EINVAL;
if (src->owner.refcount < 1L || src->owner.size == 0)
return errno = EINVAL;
else {
++(src->owner.refcount);
m->owner = src->owner;
}
m->origin = src->origin;
m->rows = src->cols;
m->cols = src->rows;
m->rowstride = src->colstride;
m->colstride = src->rowstride;
return 0;
}
Using code similar to above, you can create one-row or one-column matrix views describing any row, column, or diagonal. (The diagonals are especially useful in certain situations.) Submatrices can be mirrored or rotated, and so on. You can safely free a matrix you only need a submatrix or other view from, as the owner structure reference count keeps track of when the data can be safely discarded.
Matrix multiplication, and other similar complex operations for larger matrices, are very sensitive to cache locality issues. This means that you are better off copying the source matrix data into compact arrays (with the arrays properly aligned and in elements in correct order for that operand). The overhead caused by both row and column having a separate stride (instead of only one, as is typical) is actually minimal; in my own tests, neglible.
The best feature of this approach, however, is that it lets you write efficient code without worrying about what is a "real" matrix, what is a "view", and how the actual underlying data is stored in an array, unless you care. Finally, it is simple enough for anyone who grasps basic dynamic memory management in C to completely understand.