(I know I'm repeating myself here, but I cannot help it! 
)
Why not use a structure to describe the matrix?
Code:
#ifndef MATRIX_H
#define MATRIX_H
#define _ISO_C99_SOURCE
#include <stdlib.h>
#include <string.h>
#include <errno.h>
/* Matrix data owner structure.
* Many matrices can refer to the same data.
* The owner is only destroyed when the last matrix
* referring to the data is destroyed.
*/
typedef struct {
size_t size;
long refcount;
double data[]; /* C99. Use data[0] if your compiler complains. */
} owner_t;
/* Matrix data structure.
*/
typedef struct {
int rows;
int cols;
long rowstep;
long colstep;
double *origin;
owner_t *owner;
} matrix_t;
/* To declare a matrix variable, use either
* matrix_t m = MATRIX_INIT;
* or
* matrix_t m;
* matrix_init(&m);
*/
#define MATRIX_INIT { 0, 0, 0, 0, NULL, NULL }
static void matrix_init(matrix_t *const mptr)
{
if (mptr)
*mptr = (matrix_t)MATRIX_INIT; /* C99. Assign each member instead if your compiler complains. */
}
/* Macro referring to a specific element. Ealuates m three times!
*/
#define MATRIX_ELEMENT(m,row,col) ((m).origin[(row)*(m).rowstep + (col)*(m).colstep])
/* "outside" value returned by matrix_get() and matrix_set,
* if the specified element is outside the matrix.
*/
static double matrix_outside = 0.0;
/* Safe matrix element getter.
*/
static inline double matrix_get(const matrix_t *const mptr, const int row, const int col)
{
if (mptr && row >= 0 && col >= 0 && row < mptr->rows && col < mptr->cols)
return mptr->origin[row * mptr->rowstep + col * mptr->colstep];
else
return matrix_outside;
}
/* Safe matrix element setter.
*/
static inline double matrix_set(matrix_t *const mptr, const int row, const int col, const double value)
{
if (mptr && row >= 0 && col >= 0 && row < mptr->rows && col < mptr->cols)
return mptr->origin[row * mptr->rowstep + col * mptr->colstep] = value;
else
return matrix_outside;
}
/* Destroy a matrix.
* Every matrix you create, you must destroy.
* It is okay to destroy the same matrix twice,
* and to destroy a matrix that's only init'ed, not created.
* (Destroying a matrix does not affect related matrices,
* even if they refer to the same data.)
*/
static void matrix_destroy(matrix_t *const mptr)
{
if (mptr) {
if (mptr->owner && --mptr->owner->refcount < 1) {
mptr->owner->size = 0;
mptr->owner->refcount = -1073741824L;
free(mptr->owner);
}
*mptr = (matrix_t)MATRIX_INIT; /* C99. Assign each member instead if your compiler complains. */
}
}
/* Create a matrix. The elements are non-initialized.
* Note: If mptr points to an existing matrix, it'll be destroyed first.
* Returns 0 if successful, errno error code otherwise.
*/
static int matrix_create(matrix_t *const mptr, const int rows, const int cols)
{
const size_t total = (size_t)rows * (size_t)cols;
const size_t bytes = total * sizeof (double);
const size_t owned = bytes + sizeof (owner_t);
owner_t *owner;
/* Verify sane arguments. */
if (!mptr || rows < 1 || cols < 1)
return errno = EINVAL;
/* Destroy existing matrix. */
matrix_destroy(mptr);
/* Verify the size requested is not too large. */
if ((size_t)total / (size_t)cols != (size_t)rows ||
(size_t)bytes / sizeof (double) != total ||
owned <= bytes)
return errno = ENOMEM;
/* Allocate the data owner structure. */
owner = malloc(owned);
if (!owner)
return errno = ENOMEM;
/* Owner owns the matrix data, */
owner->size = total;
owner->refcount = 1L;
/* and the matrix just refers to it. */
mptr->owner = owner;
mptr->origin = owner->data;
mptr->rows = rows;
mptr->cols = cols;
/* By default, use row-major order */
mptr->rowstep = cols;
mptr->colstep = 1L;
return 0;
}
/* Transpose a matrix.
* Returns 0 if successful, errno error code otherwise.
*/
int matrix_transpose(matrix_t *const mptr)
{
if (mptr && mptr->rows > 0 && mptr->cols > 0) {
const int old_rows = mptr->rows;
const int old_cols = mptr->cols;
const long old_rowstep = mptr->rowstep;
const long old_colstep = mptr->colstep;
mptr->rows = old_cols;
mptr->cols = old_rows;
mptr->rowstep = old_colstep;
mptr->colstep = old_rowstep;
return 0;
} else
return errno = EINVAL;
}
/* Detach a matrix.
* If there are views to m, or m itself is a view,
* this duplicates the data, detaching all views.
* Returns 0 if successful, errno error code otherwise.
*/
int matrix_detach(matrix_t *mptr)
{
matrix_t temp = MATRIX_INIT;
int row, col;
/* No matrix specified? */
if (!mptr || mptr->rows < 1 || mptr->cols < 1)
return errno = EINVAL;
/* Create a new matrix large enough to hold m. */
if (matrix_create(&temp, mptr->rows, mptr->cols))
return errno;
/* Copy matrix contents. */
for (row = 0; row < temp.rows; row++)
for (col = 0; col < temp.cols; col++)
temp.origin[row * temp.rowstep + col * temp.colstep] =
mptr->origin[row * mptr->rowstep + col * mptr->colstep];
/* Destroy the old m, */
matrix_destroy(mptr);
/* replacing it with the new one. */
*mptr = temp;
return 0;
}
/* Set matrix d to reflect the contents of matrix m.
* Any changes in m are visible in d, and vice versa;
* destroying one does not affect the other.
* You can transpose one without affecting the other.
* Returns 0 if success, errno error code otherwise.
*/
int matrix_attach_all(matrix_t *const dptr, const matrix_t *const mptr)
{
if (dptr) {
matrix_destroy(dptr);
if (mptr && mptr->origin && mptr->owner->refcount++ > 0) {
*dptr = *mptr;
return 0;
} else
return errno = ENOENT; /* No matrix to view. */
} else
return errno = EINVAL; /* No d matrix specified. */
}
/* Set matrix d to reflect (partial) contents of matrix m.
* For example, you can set d to be a single row, column,
* or even diagonal.
* Returns 0 if success, errno error code otherwise.
*/
int matrix_attach(matrix_t *const dptr, const matrix_t *const mptr,
const int rows, const int cols,
const int origin_row, const int origin_col,
const long row_rowstep, const long row_colstep,
const long col_rowstep, const long col_colstep)
{
int lastrow_row, lastrow_col;
int lastcol_row, lastcol_col;
/* No d specified? Invalid parameters? */
if (!dptr || rows < 0 || cols < 0 || origin_row < 0 || origin_col < 0)
return errno = EINVAL;
matrix_destroy(dptr);
/* No m specified? */
if (!mptr || !mptr->origin)
return errno = ENOENT;
/* Origin of view outside m? */
if (origin_row >= mptr->rows || origin_col >= mptr->cols)
return errno = EINVAL;
lastrow_row = origin_row + (rows - 1) * row_rowstep;
lastrow_col = origin_col + (rows - 1) * row_colstep;
lastcol_row = origin_row + (cols - 1) * col_rowstep;
lastcol_col = origin_col + (cols - 1) * col_colstep;
/* Does the view extend outside m? */
if (lastrow_row < 0 || lastrow_row >= mptr->rows ||
lastrow_col < 0 || lastrow_col >= mptr->cols ||
lastcol_row < 0 || lastcol_row >= mptr->rows ||
lastcol_col < 0 || lastcol_col >= mptr->cols)
return errno = EINVAL;
/* d will refer to the same data. */
if (mptr->owner->refcount++ < 1)
return errno = ENOENT;
dptr->rows = rows;
dptr->cols = cols;
dptr->rowstep = row_rowstep * mptr->rowstep + row_colstep * mptr->colstep;
dptr->colstep = col_rowstep * mptr->rowstep + col_colstep * mptr->colstep;
dptr->origin = mptr->origin + origin_row * mptr->rowstep + origin_col * mptr->colstep;
dptr->owner = mptr->owner;
return 0;
}
#endif /* MATRIX_H */
I only use C99-capable compilers, but if you must work with C89/C90, then small changes are needed to make the above compile; they're commented in the code. (The data[0] is a hack, but it's a widely used one; that's why it was standardized later/better.)
Overall, the matrix representation it uses is much more powerful than the representation used by e.g. GNU Scientific Library.
While the matrix.h structures do use one extra multiplication per element access (compared to e.g. GSL), the versatility it allows in your algorithm implementation is usually worth it. This was not designed for maximum execution speed, but but maximum versatility and ease of use. (It does not mean it's slow, though; it shouldn't be measurably slower except on very, very large matrices.)
Typically, I leverage this for submatrices and diagonal vectors. However, any regular pattern from the original matrix that stays within the original matrix is possible -- even something like every third column starting from right, and every second element in each column.
If you do e.g. linear algebra, and you need matrix views that have extra elements, create a large enough "canvas" matrix first. Create your data matrix as a (centered) view to it. Then, you can create more views that extend past the data area, as long as they stay within the "canvas" matrix.
Here is a simple example:
Code:
#include <stdio.h>
#include <errno.h>
#include "matrix.h"
/* Read a matrix from a stream.
*
* The format of the input starts with two integers,
* ROWS COLS
* followed by the matrix element in natural order,
* upper left corner first, followed by the rest of the elements
* on the top row, followed by the second row, and so on.
*
* Returns 0 if success, errno error code otherwise. Special errors:
* ENOENT: No data, end of input.
* EDOM: Invalid number of rows or columns.
* EIO: Error reading matrix elements.
*/
int matrix_read(matrix_t *const mptr, FILE *const in)
{
int rows, cols, row, col;
/* No m specified? */
if (!mptr)
return errno = EINVAL;
matrix_destroy(mptr);
/* No input stream specified? */
if (!in)
return errno = EINVAL;
/* Skip leading whitespace. */
(void)fscanf(in, " ");
/* Read error? End of input? */
if (ferror(in))
return errno = EIO;
if (feof(in))
return errno = ENOENT;
/* Parse the number of rows and cols. */
if (fscanf(in, " %d %d", &rows, &cols) != 2)
return errno = EDOM;
if (rows < 1 || cols < 1)
return errno = EDOM;
/* Create a new matrix of desired size. */
if (matrix_create(mptr, rows, cols))
return errno;
/* Read loop. */
for (row = 0; row < rows; row++)
for (col = 0; col < cols; col++)
if (fscanf(in , " %lf", &(mptr->origin[row * mptr->rowstep + col * mptr->colstep])) != 1) {
matrix_destroy(mptr);
return errno = EIO;
}
/* Success. */
return 0;
}
/* Pretty-print a matrix using pattern (for each double).
*/
void matrix_fprintf(FILE *const out, const char *const pattern, const matrix_t *const mptr)
{
int r, c;
for (r = 0; r < mptr->rows; r++) {
fprintf(out, "[");
for (c = 0; c < mptr->cols; c++)
fprintf(out, pattern, mptr->origin[r * mptr->rowstep + c * mptr->colstep]);
fprintf(out, " ]\n");
}
}
#define matrix_printf(pattern, mptr) matrix_fprintf(stdout, pattern, mptr)
int main(void)
{
matrix_t m = MATRIX_INIT;
matrix_t d = MATRIX_INIT;
int n;
/* Read m from standard input. */
if (matrix_read(&m, stdin)) {
fprintf(stderr, "Error reading a matrix from standard input.\n");
return 1;
}
/* Print the matrix. */
printf("Read %d-row, %d-column matrix from standard input:\n", m.rows, m.cols);
matrix_printf(" %5g", &m);
/* Create a column vector of the diagonal elements. */
if (m.rows <= m.cols)
n = m.rows;
else
n = m.cols;
if (matrix_attach(&d, &m,
1, n, /* Column vector of n elements */
0, 0, /* Starts at the upper left corner of the matrix */
0L, 0L, /* Since there is just one row, the row step is empty */
1L, 1L /* Each vector element is one cell down and right in the matrix */
)) {
fprintf(stderr, "Cannot create a diagonal view.\n");
return 1;
}
printf("\nDiagonal elements in the matrix:\n");
matrix_printf(" %5g", &d);
/* Negate a middle element (of the first row) in the diagonal vector. */
MATRIX_ELEMENT(d, 0, n/2) *= -1.0;
printf("\nAfter negating the middle element in the diagonal vector:\n");
matrix_printf(" %5g", &d);
/* We don't need the diagonal vector anymore, so destroy it. */
matrix_destroy(&d);
/* Destroying one matrix does not affect the others,
* even if they refer to the same data. */
printf("\nThe original matrix reflects the change:\n");
matrix_printf(" %5g", &m);
matrix_destroy(&m);
return 0;
}
When running the example, supply it with matrix size (number of rows and columns) and the matrix elements in standard input; for example, if you supply it with
Code:
5 5
11 12 13 14 15
21 22 23 24 25
31 32 33 34 35
41 42 43 44 45
51 52 53 54 55
it will output
Code:
Read 5-row, 5-column matrix from standard input:
[ 11 12 13 14 15 ]
[ 21 22 23 24 25 ]
[ 31 32 33 34 35 ]
[ 41 42 43 44 45 ]
[ 51 52 53 54 55 ]
Diagonal elements in the matrix:
[ 11 22 33 44 55 ]
After negating the middle element in the diagonal vector:
[ 11 22 -33 44 55 ]
The original matrix reflects the change:
[ 11 12 13 14 15 ]
[ 21 22 23 24 25 ]
[ 31 32 -33 34 35 ]
[ 41 42 43 44 45 ]
[ 51 52 53 54 55 ]
