Originally Posted by
oogabooga
BTW, your comment on double *data says rows where it should say row
Good catch! Absolutely right. As you had fixed in the quote, the comment should be /* data[row*rowstep + col*colstep] */.
Originally Posted by
oogabooga
Also, you supposedly shouldn't end your own type names with _t.
It is true that the standard C library developers recommend that, but it's a "soft" recommendation, a warning of what kind of names they might need to reserve in the future. In other words, I read their recommendation as "if we define a new type, it will be named with a _t suffix; if we define a new error code, it will be defined with an E prefix", and so on.
For libraries, say for example "Nominal Animal's Library", or nal, the recommended type names are nal_matrix and nal_owner. The library abbreviation prefix helps avoid namespace pollution, although I have a nagging feeling nal_ is already used ..
For example code, I often use _t suffix for types to make the code easier to read.
If you have a better type name prefix or suffix I should use instead -- as some kind of suffix or prefix is needed to distinguish type names from variable names, making example code easier to read --, I'd be more than happy to hear suggestions!
Originally Posted by
oogabooga
Can you give an example of using the owner field?
Sure. For simplicity, let's keep this to single-threaded code.
That is, you can use these safely from multiple threads, as long as you manipulate a matrix (or shared matrices; that is, owners) from a single thread at a time.
First, we need the types. I'll also define safe accessors, in cases where you might have out-of-bounds column or row numbers:
Code:
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
typedef struct {
long refcount;
size_t size;
double data[];
} nal_owner;
typedef struct {
long rows;
long cols;
long rowstep;
long colstep;
double *data;
nal_owner *owner;
} nal_matrix;
/* Safe matrix accessor (getter).
* Returns 0.0 if row,col is outside the matrix m, otherwise the matrix element.
*/
double nal_matrix_get(const nal_matrix *const m, const long row, const long col)
{
if (row >= 0L && col >= 0L && m && row < m->rows && col < m->cols && m->data)
return m->data[row * m->rowstep + col * m->colstep];
else
return 0.0;
}
/* Safe matrix accessor (setter).
* Returns 0.0 if row,col is outside the matrix m, otherwise the value set.
*/
double nal_matrix_set(nal_matrix *const m, const long row, const long col, const double value)
{
if (row >= 0L && col >= 0L && m && row < m->rows && col < m->cols && m->data)
return m->data[row * m->rowstep + col * m->colstep] = value;
else
return 0.0;
}
The getter and setter are more for documentation, rather than actual use.
Next, we need a function that creates a data owner of requested size. This one tries to verify that the size does not overflow, but as I rewrote this one from scratch, it might be buggy.
Code:
/* Create a new owner structure, large enough to hold a matrix of the specified size.
*/
static nal_owner *nal_owner_create(const long row_count, const long col_count)
{
const size_t rows = (size_t)row_count,
cols = (size_t)col_count;
const size_t size = rows * cols;
const size_t bytes = size * sizeof (double);
const size_t total = bytes + sizeof (nal_owner);
nal_owner *owner;
/* Size too small? */
if (row_count < 1L || col_count < 1L) {
errno = EINVAL;
return NULL;
}
/* Size too large? */
if (size / rows != cols ||
bytes / sizeof (double) != size ||
bytes >= total) {
errno = ENOMEM;
return NULL;
}
owner = malloc(total);
if (!owner) {
errno = ENOMEM;
return NULL;
}
owner->refcount = 1L; /* Referenced by caller! */
owner->size = size; /* row_count*col_count doubles */
return owner;
}
Every time a matrix is no longer needed, it must be destroyed. Here is the destructor. Note that the owner is only destroyed when the refcount drops to zero.
Code:
/* Free/destroy a matrix.
* Does not affect other matrices sharing the data.
*/
void nal_matrix_destroy(nal_matrix *const m)
{
if (!m)
return;
if (m->owner)
if (m->owner->refcount-- < 1L) {
m->owner->size = 0;
free(m->owner);
}
m->rows = 0L;
m->cols = 0L;
m->rowstep = 0L;
m->colstep = 0L;
m->data = NULL;
m->owner = NULL;
}
When an independent matrix is created, the associated owner structure is also created:
Code:
/* Create a new independent matrix.
* Returns 0 if success, errno error code otherwise
* (EINVAL or ENOMEM).
*/
int nal_matrix_create(nal_matrix *const m, const long rows, const long cols)
{
nal_owner *owner;
if (!m)
return errno = EINVAL;
m->rows = 0L;
m->cols = 0L;
m->rowstep = 0L;
m->colstep = 0L;
m->data = NULL;
m->owner = NULL;
if (rows < 1L || cols < 1L)
return errno = EINVAL;
owner = nal_owner_create(rows, cols);
if (!owner)
return errno;
m->rows = rows;
m->cols = cols;
m->rowstep = cols;
m->colstep = 1;
m->data = owner->data;
m->owner = owner;
return 0;
}
Certain operations, like transpose, are trivial:
Code:
/* Transpose matrix m.
* Returns 0 if success, errno otherwise.
*/
int nal_matrix_transpose(nal_matrix *const m)
{
if (m) {
const long rows = m->rows,
cols = m->cols,
rowstep = m->rowstep,
colstep = m->colstep;
m->rows = cols;
m->cols = rows;
m->rowstep = colstep;
m->colstep = rowstep;
return 0;
} else
return errno = EINVAL;
}
We can create a linked duplicate of a matrix; that is, create a matrix that shares its data with another matrix:
Code:
/* Duplicate matrix, keeping the data shared between the two.
*/
int nal_matrix_dup(nal_matrix *const d, const nal_matrix *const m)
{
if (!d)
return errno = EINVAL;
if (!m || m->rows < 1L || m->cols < 1L || !m->owner || m->owner->refcount < 1L) {
d->rows = 0L;
d->cols = 0L;
d->rowstep = 0L;
d->colstep = 0L;
d->data = NULL;
d->owner = NULL;
return errno = ENOENT;
}
m->owner->refcount++;
d->rows = m->rows;
d->cols = m->cols;
d->rowstep = m->rowstep;
d->colstep = m->colstep;
d->data = m->data;
d->owner = m->owner;
return 0;
}
As you can see, duplicating a matrix, or referring to the data in a matrix from another matrix, involves just incrementing the refcount of the data owner.
In a nutshell: the owner owns the data, and the matrices just refer to it.
A plain duplicate/shared view is very dull; that's why we have rowstep and colstep. They allow very interesting shared matrices to be created:
Code:
/* Set d to be a row vector of matrix m.
* Returns 0 if success, errno otherwise.
*/
int nal_matrix_row_of(nal_matrix *const d, const nal_matrix *const m, const long row)
{
if (!d)
return errno = EINVAL;
d->rows = 0L;
d->cols = 0L;
d->rowstep = 0L;
d->colstep = 0L;
d->data = NULL;
d->owner = NULL;
if (!m || !m->owner || m->owner->size < 1)
return errno = EINVAL;
if (row < 0L || row >= m->rows)
return errno = EINVAL;
m->owner->refcount++;
d->rows = 1;
d->cols = m->cols;
d->rowstep = 0;
d->colstep = m->cols;
d->owner = m->owner;
d->data = m->data + row * m->rowstep;
return 0;
}
/* Set d to be a column vector of matrix m.
* Returns 0 if success, errno otherwise.
*/
int nal_matrix_column_of(nal_matrix *const d, const nal_matrix *const m, const long col)
{
if (!d)
return errno = EINVAL;
d->rows = 0L;
d->cols = 0L;
d->rowstep = 0L;
d->colstep = 0L;
d->data = NULL;
d->owner = NULL;
if (!m || !m->owner || m->owner->size < 1)
return errno = EINVAL;
if (col < 0L || col >= m->cols)
return errno = EINVAL;
m->owner->refcount++;
d->rows = m->rows;
d->cols = 1;
d->rowstep = m->rowstep;
d->colstep = 0;
d->owner = m->owner;
d->data = m->data + col * m->colstep;
return 0;
}
/* Set d to be a part of matrix m.
* Returns 0 if success, errno otherwise.
*/
int nal_matrix_part_of(nal_matrix *const d, const nal_matrix *const m,
const long row, const long col,
const long rows, const long cols,
const long down_rowstep, const long down_colstep,
const long right_rowstep, const long right_colstep)
{
const long lastrow = rows - 1L,
lastcol = cols - 1L;
const long upper_right_row = row + lastcol * right_rowstep,
upper_right_col = col + lastcol * right_colstep,
lower_left_row = row + lastrow * down_rowstep,
lower_left_col = col + lastcol * down_colstep,
lower_right_row = row + lastrow * down_rowstep + lastcol * right_rowstep,
lower_right_col = col + lastrow * down_colstep + lastcol * right_colstep;
if (!d)
return errno = EINVAL;
d->rows = 0L;
d->cols = 0L;
d->rowstep = 0L;
d->colstep = 0L;
d->data = NULL;
d->owner = NULL;
if (!m || !m->owner || m->owner->size < 1)
return errno = EINVAL;
/* Empty? */
if (rows < 1L || cols < 1L)
return errno = EINVAL;
/* Is the entire part within the matrix m? */
if (row < 0L || row >= m->rows ||
col < 0L || col >= m->cols ||
upper_right_row < 0L || upper_right_row >= m->rows ||
upper_right_col < 0L || upper_right_col >= m->cols ||
lower_left_row < 0L || lower_left_row >= m->rows ||
lower_left_col < 0L || lower_left_col >= m->cols ||
lower_right_row < 0L || lower_right_row >= m->rows ||
lower_right_col < 0L || lower_right_col >= m->cols)
return errno = EINVAL;
m->owner->refcount++;
d->rows = rows;
d->cols = cols;
d->rowstep = m->rowstep * down_rowstep + m->colstep * down_colstep;
d->colstep = m->rowstep * right_rowstep + m->colstep * right_colstep;
d->data = m->data + row * m->rowstep + col * m->colstep;
d->owner = m->owner;
return 0;
}
/* Set d to be the diagonal of matrix m.
* Returns 0 if success, errno otherwise.
*/
int nal_matrix_diagonal_of(nal_matrix *const d, const nal_matrix *const m)
{
if (!d)
return errno = EINVAL;
d->rows = 0L;
d->cols = 0L;
d->rowstep = 0L;
d->colstep = 0L;
d->data = NULL;
d->owner = NULL;
if (!m || !m->owner || m->owner->size < 1)
return errno = EINVAL;
m->owner->refcount++;
if (m->rows <= m->cols)
d->rows = m->rows;
else
d->rows = m->cols;
d->cols = 1;
d->rowstep = m->rowstep + m->colstep;
d->colstep = 0;
d->owner = m->owner;
d->data = m->data;
return 0;
}
While nal_matrix_row_of(), nal_matrix_column_of(), and nal_matrix_diagonal_of() could be easily written as a wrappers around nal_matrix_part_of(), I wanted to have a couple more examples of interesting shared matrix manipulations, in case it would help others to understand the logic.
To help debugging stuff, we'll need to print the matrices:
Code:
void nal_matrix_fprint(FILE *const out, const nal_matrix *const m,
const char *const prefix,
const char *const format,
const char *const suffix)
{
if (out && m && m->data) {
long r, c;
for (r = 0L; r < m->rows; r++) {
if (prefix && *prefix)
fputs(prefix, out);
for (c = 0L; c < m->cols; c++) {
if (c > 0L)
fputc(' ', out);
fprintf(out, format, m->data[r * m->rowstep + c * m->colstep]);
}
if (suffix && *suffix)
fputs(suffix, out);
fputc('\n', out);
}
}
}
A common format for matrices is to have the number of rows and number of columns precede the actual data. Here is a function that reads such matrices:
Code:
int nal_matrix_fread(FILE *const in, nal_matrix *const m)
{
long rows, cols, row, col;
if (!in || !m)
return errno = EINVAL;
if (fscanf(in, " %ld %ld", &rows, &cols) != 2 || rows < 1L || cols < 1L)
return errno = EDOM;
if (nal_matrix_create(m, rows, cols))
return errno;
for (row = 0L; row < rows; row++)
for (col = 0L; col < cols; col++)
if (fscanf(in, " %lf", &m->data[row * m->rowstep + col * m->colstep]) != 1)
return errno = EDOM;
return 0;
}
For illustration, here is a main() that reads a matrix (number of rows, number of columns, followed by each matrix element in row-major order, i.e. western written test order) from standard input, and manipulates it a bit.
Code:
int main(void)
{
nal_matrix m, m1, m2, m3;
long i;
/* Read a matrix from standard input.
* It is preceded by the number of rows, and the number of cols.
*/
if (nal_matrix_fread(stdin, &m)) {
if (errno == EDOM)
fprintf(stderr, "Invalid matrix data.\n");
else
fprintf(stderr, "%s.\n", strerror(errno));
return 1;
}
printf("Read a %ld-row, %ld-column matrix:\n", m.rows, m.cols);
nal_matrix_fprint(stdout, &m, "\t[ ", "%9.4f", " ]");
printf("\n");
if (nal_matrix_diagonal_of(&m1, &m)) {
fprintf(stderr, "nal_matrix_diagonal_of(): %s.\n", strerror(errno));
return 1;
}
printf("Column vector of diagonal elements:\n");
nal_matrix_transpose(&m1);
nal_matrix_fprint(stdout, &m1, "\t[ ", "%9.4f", " ]");
printf("\n");
if (nal_matrix_dup(&m2, &m)) {
fprintf(stderr, "nal_matrix_dup(): %s.\n", strerror(errno));
return 1;
}
if (!nal_matrix_part_of(&m3, &m, m.rows - 2, m.cols - 2,
m.rows - 2, m.cols - 2,
-1, 0,
0, -1)) {
printf("Center elements mirrored:\n");
nal_matrix_fprint(stdout, &m3, "\t[ ", "%9.4f", " ]");
printf("\n");
nal_matrix_destroy(&m3);
}
if (!nal_matrix_part_of(&m3, &m, 0, m.cols / 2, 1+m.rows/2, 1+m.cols/2, 1, -1, 1, 1)) {
printf("Checkerboard odd tiles, rotated 45 degrees:\n");
nal_matrix_fprint(stdout, &m3, "\t[ ", "%9.4f", " ]");
printf("\n");
nal_matrix_destroy(&m3);
}
nal_matrix_destroy(&m);
for (i = 0L; i < m1.cols; i++)
m1.data[i * (m1.rowstep + m1.colstep)] = -m1.data[i * (m1.rowstep + m1.colstep)];
printf("Original matrix transposed, diagonal elements negated:\n");
nal_matrix_transpose(&m2);
nal_matrix_fprint(stdout, &m2, "\t[ ", "%9.4f", " ]");
printf("\n");
nal_matrix_destroy(&m2);
printf("Negated diagonal elements as a column vector:\n");
nal_matrix_fprint(stdout, &m1, "\t[ ", "%9.4f", " ]");
printf("\n");
nal_matrix_destroy(&m1);
return 0;
}
As long as you nal_matrix_destroy() every matrix you create, you don't need to worry about which matrices are shared; there is no distinction between shared and unique matrices.
If thread safety is wanted, then the owner structures need to be chained, so that instead of immediately freed, they are only moved to a garbage chain. A couple of functions are used to manipulate an atomic counter (or semaphore, although semaphores may have a pretty low maximum count), say nal_lock() and nal_unlock(), which are used around accesses to matrix data via a temporary variable (and inside many of the library functions above); the garbage chain is then emptied whenever nal_unlock() is called without any other nal_lock()s pending. Some fields, like owner refcount, need either locking or atomic accesses, too.
More complex operations, like matrix multiplication, benefits from copying the matrix data to temporary storage that optimizes access patterns. The overhead of the above nal_matrix structure compared to those used by most libraries like GSL, that have only a stride (index = column + row * stride), seems to be neglible for all but smallest matrices.
If performance is important, smallest matrices -- especially 2x2, 3x3, 4x4 and 3x4 -- are best special-cased anyway.
The above example code has lots of defensive setting/clearing of structure fields that can safely be omitted in a real library. There are also lots of optimizations that can be done, and probably quite a few bugs too; I did just rewrite the above from scratch. (I'm terrible at finding copies of stuff I've already posted. I know I've posted the above already, but I just couldn't find it.)
Everyone is absolutely free to use the code shown in any way they like. I don't think it is creative enough to be copyrightable, but if it is, then it is free for any use, as long as you don't blame me if it breaks.
Questions?