这不是对所述问题的回答,而是对根本问题的回答:以最少的努力管理矩阵和矩阵视图。
这将获得反对票,但每当提出 OP 提出的问题类型时,它在解决潜在问题方面非常有用,我发现值得在这里展示这种替代方法。
对于固定大小的小矩阵来说,这并不有趣,因为这些特征仅在大小较大或变化时才显示出它们的好处。
我使用以下两种结构来描述矩阵。为了简单起见,我将省略内存池支持(它允许将一组矩阵作为一个池来管理,一次释放它们,而不必单独管理每个矩阵)以及与多线程操作和线程安全相关的所有内容。
代码可能包含拼写错误;如果您发现任何问题,请发表评论,我会修复它们。
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 }
矩阵m元素在行r,列c,是m.origin[r * m.rowstride + c * m.colstride],假设0 <= r && r < m.rows和0 <= c < m.cols。
矩阵通常被声明为局部变量,而不是指针。您确实需要记住在不再需要后释放每个单独的矩阵。 (我省略的池机制可以避免这种情况,因为池中的所有矩阵都会立即释放。)
每个矩阵仅指一个所有者结构。所有者结构记录引用的数量(引用该结构中数据的矩阵的数量),并在引用计数降至零时释放:
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;
}
}
每当创建新矩阵时,都会创建相应的所有者结构:
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;
}
请注意,上面没有将矩阵元素初始化为任何值,因此它们最初包含垃圾。
矩阵转置是一种简单而快速的操作:
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;
}
}
同样,您可以旋转和镜像矩阵,只需记住在这些情况下也修改 origin 成员。
有趣且有用的案例是能够在其他矩阵中创建“视图”。引用的数据完全相同——修改一个在另一个中立即可见;这是真正的别名——不需要内存复制。与大多数库(例如 GSL、GNU 科学库)不同,这些“视图”本身就是完全普通的矩阵。下面是一些例子:
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;
}
使用与上述类似的代码,您可以创建描述任何行、列或对角线的单行或单列矩阵视图。 (对角线在某些情况下特别有用。)
子矩阵可以被镜像或旋转,等等。
您可以安全地释放只需要子矩阵或其他视图的矩阵,因为所有者结构引用计数会跟踪何时可以安全地丢弃数据。
矩阵乘法和其他类似的大型矩阵复杂运算对缓存局部性问题非常敏感。这意味着您最好将源矩阵数据复制到紧凑的数组中(数组正确对齐,并且该操作数的元素顺序正确)。由具有单独步幅的行和列(而不是典型的只有一个)引起的开销实际上是最小的;在我自己的测试中,可以忽略不计。
然而,这种方法的最大特点是它可以让您编写高效的代码,而不必担心什么是“真实”矩阵、什么是“视图”以及实际的底层数据如何存储在数组中, 除非你在乎。
最后,对于掌握 C 中基本动态内存管理的任何人来说,它都足够简单,完全可以理解。