C中的二元向量和矩阵操作
我试图在C中实现一个数据结构，这将允许我有效地操作**二进制**矩阵（仅包含1或0）。 我将解释我必须对此矩阵应用哪些操作，并想知道使用哪种最佳数据结构？
操作在字段F_2中完成（这意味着1 + 1 = 0，其他操作保持不变）。 我有一个叫做H k * n矩阵（ k < n ）。 最多， k = 2325， n = 3009。
我必须对此矩阵执行的操作是：
我将仅使用行交换和行添加来部分对角化它。 一旦完成，我将不再使用行操作，并将在此矩阵上运行大量（！）列添加 （我的意思是“很多”是关于（（nk）/ 2）³列添加）
我正在考虑矩阵的数据结构：
对于矩阵系数，我考虑在一个单个unsigned int中一次存储多个位的序列。 例如，我可以将序列(11001011)存储到uint8_t 203 （从二进制转换为十进制）
- 这是个好主意吗 ？
如果我这样做，我有两个选择：
我可以使用uint16_t或uint64_t系数在许多4 * 4或8 * 8子矩阵中分割我的矩阵H.
- 这是一个很好的选择（在时间效率方面），如果是，是否更好地使用
uint16_t或uint64_t？
另外我在考虑将每一行存储在多个uint32_t或uint64_t ，然后操作我的部分对角化。 接下来切换到将矩阵编码为n列向量以处理剩余操作的结构。
- 你认为这更有效吗？
无论我使用什么方法，我都必须有效地访问unsigned int的第n位（ uint16或64 ）。 我怎么做 ？
为获得最佳性能，请使用行指针数组进行行交换和行添加。 使用uint_fast32_t ，除非您打算在16位或8位处理器上运行它。
完成所有行交换和行添加后，转置数组。 虽然此操作“慢”，但以下列操作将如此快以抵消转置成本。
考虑以下：
#include #include typedef uint_fast32_t word; #define WORD_BITS (CHAR_BIT * sizeof (word)) typedef struct { int rows; /* Number of row vectors */ int cols; /* Number of defined bits in each row */ int words; /* Number of words per row vector */ word **row; /* Array of pointers */ } row_matrix; typedef struct { int rows; /* Number of defined bits in each column */ int cols; /* Number of col vectors */ int words; /* Number of words per col vector */ word **col; } col_matrix;
虽然您可以使用单一类型来描述两种矩阵forms，但使用单独的类型可以使代码和函数更易于维护。 你最终会得到一些重复的代码，但与清晰，直观的类型相比，这是一个小问题。
在32位系统上， uint_fast32_t通常是32位类型。 在64位系统上，它通常是64位。 WORD_BITS宏扩展为word的位数 – 它不总是32！
对位进行编号的最简单方法是将矩阵中最左边的位指定为位0，并将这些位存储在每个字的最低有效位中。 如果你有row_matrix *rm ，那么行row ， col列的位是
!!(rm->row[row][col / WORD_BITS] & ((word)1U << (col % WORD_BITS)))
!! 是not-not运算符：如果参数非零，则产生1，否则它产生0.因为我们从字中屏蔽了一个位，所以“bit is set”值将是2的幂（1,2） ，4,8,16,32,64等）。
要设置该位，请使用
rm->row[row][col / WORD_BITS] |= (word)1U << (col % WORD_BITS);
要清除一点，你需要使用除目标位1之外的所有掩码进行二进制AND运算。使用not运算符很容易实现~ ：
rm->row[row][col / WORD_BITS] &= ~((word)1U << (col % WORD_BITS));
col_matrix *cm的相应操作是
!!(cm->col[col][row / WORD_BITS] & ((word)1U << (row % WORD_BITS))) cm->col[col][row / WORD_BITS] |= (word)1U << (row % WORD_BITS); cm->col[col][row / WORD_BITS] &= ~((word)1U << (row % WORD_BITS));
虽然除法/和模数（或余数） %通常是慢速的（与加法，减法甚至乘法相比），但这里WORD_BITS将是所有广泛使用的架构上的两个编译WORD_BITS量的幂。 我所知道的所有编译器都会将上述内容转换为快速位移和二进制AND运算符。
要将行srcrow添加到行dstrow ，您只需对所有单词执行二进制异或：
{ const word *const src = rm->row[srcrow]; word *const dst = rm->row[dstrow]; int w = rm->words; while (w-->0) dst[w] ^= src[w]; }
类似地，对于列矩阵，
{ const word *const src = cm->col[srccol]; word *const dst = cm->col[dstcol]; int w = cm->words; while (w-->0) dst[w] ^= src[w]; }
请注意，如果组合两行以上，则可以非常有效地执行此操作; 它会比连续添加更快。 Intel和AMD CPU非常擅长预测上述模式，因此您可以使用多个源行/列。 此外，目的地不必参与结果，但如果我猜对了你正在实施的算法，我猜你想要它。
如果你知道目标体系结构有SSE2或更好，甚至AVX，你可以分别使用emmintrin.h或immintrin.h头文件，用于编译器内置类型和运算符，允许你XOR 128位和256位，分别，立刻; 有时会给你很大的提升。
由于向量类型需要C标准称之为“过度对齐”，因此您还需要包含mm_malloc.h ，并使用_mm_malloc()和_mm_free()来分配word数据的行/列向量 - 和显然是圆的words ，所以你可以访问行/列作为一个合适的整数字类型（ __m128i为SSE ， __m256i为AVX）。
就个人而言，我总是首先实现非版本化版本，然后是unit testing的一些“讨厌”测试用例，然后才能看到它的矢量化。 这样做的好处是，您可以将非版本化版本作为初步版本提供给那些将使用它的人，并且您可以比较矢量化和非矢量化案例之间的测试用例结果，以查看其中一个或另一个是否有错误。
转置操作非常简单，虽然我建议使用三重循环：最内层循环一个字中的位。 此外，您可能想要检查哪个顺序 - 行或列主要 - 最适合外循环; 根据矩阵大小，您可能会看到巨大的差异。 （这是由于缓存行为：您希望CPU能够预测访问模式，而不必重新加载相同的缓存行。在最好的情况下，在最近几年的AMD和Intel x86- 64个处理器，如果两个矩阵都适合缓存，则可以获得接近缓存的速度。）
所有上述内容都可以在单个头文件中实现 - 如果目标体系结构支持SSE2 / AVX，甚至包括矢量化版本 - 因此实现起来应该不会太难。
有问题吗？
当你引用类型uint16_t ， uint64_t ，…我想这是矩阵系数。 因此，您应该比我们更了解您正在操纵哪个值：如果您可能生成大量数字，那么您需要一个大型类型来避免溢出。 关于效率我怀疑你会感觉到速度方面的差异，但你可以通过选择较小的类型来节省一些空间。
无论如何，这都是关于优化的：你不应该为强类型而烦恼。 首先，使用char （或uint8_t ）应该没问题，因为你只处理1和0 。
我没有看到从char matrix[][]切换到typedef struct matrix_columns兴趣，我认为你可以通过明智地使用行和列索引来执行操作。
最后，在unsigned int coef得到位置i的位：
unsigned int bit = (coef>>i) & 1;
[编辑]此外，这是一个关于二进制矩阵运算（乘法，加法，XOR）的相关问题 。
我想在多个uintX_t存储每一行是个好主意，但我会选择X来匹配处理器的字大小。 通过这种方式，您可以一次汇总多个矩阵条目：
uint8_t a = 3; // 00000011 uint8_t b = 2; // 00000010 b ^= a; //results 00000001 //for XOR is just addition in F_2
我认为你的问题的瓶颈是列操作，对吧？ 因此，在执行部分对角化之后，为什么不在其转置上操作列添加，就像你为部分对角化做的那样？ 你会花一些时间进行这种转换，因为它必须稍微完成，但在那之后，列的添加会简单得多。
假设您有一个64×64矩阵：
uint64_t matrix[64]; // This means that matrix[i] is the i-th row
你如何将第j列与第l列相加（mod 2），对于j和l在1和64之间？ 您将不得不迭代所有行（假设矩阵行中的较高位置处于较低有效位）：
int i; for (i=0; i<64; i++){ matrix[i] ^= (matrix[i] & (1<<(64-j))) << (64-l); }
但是如果你将matrix转换为matrixT ，则上面的for循环将等效于：
matrixT[l] ^= matrixT[j];
这将由处理器在一个步骤中完成（我猜，但不确定）。 因此，如果您在此之后花时间进行转置，您将受益于处理器在一步中对其字大小的数据进行算术运算的能力。