For best performance, use an array of row pointers for the row swap and row additions. Use <stdint.h>, and a fast unsigned integer type of minimum supported word size -- I recommend either uint_fast32_t, unless you intend to run this on 16- or 8-bit processors.
When all the row swap and row additions are done, transpose the array. Although this operation is "slow", the following column operations will be so fast to offset the transpose cost.
Consider the following:
#include <stdint.h>
#include <limits.h>
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;
Although you could use a single type to describe the two matrix forms, using separate types makes the code and functions easier to maintain. You'll end up having some duplicated code, but that's a minor issue compared to having clear, intuitive types.
On 32-bit systems, uint_fast32_t is typically a 32-bit type. On 64-bit systems, it is typically 64-bit. The WORD_BITS macro expands to the number of bits in a word -- it is NOT always 32!
The easiest way to number the bits is to designate the leftmost bit in the matrix as bit 0, and store the bits in the least significant bits in each word. If you have row_matrix *rm, then the bit at row row, column col is
!!(rm->row[row][col / WORD_BITS] & ((word)1U << (col % WORD_BITS)))
The !! is the not-not operator: if the argument is nonzero, it yields 1, otherwise it yields 0. Because we mask a single bit from the word, the "bit is set" value would otherwise be some power of two (1, 2, 4, 8, 16, 32, 64, and so on).
To set the bit, use
rm->row[row][col / WORD_BITS] |= (word)1U << (col % WORD_BITS);
To clear a bit, you need to binary-AND with a mask having all except the target bit 1. This is easy to achieve using the not operator ~:
rm->row[row][col / WORD_BITS] &= ~((word)1U << (col % WORD_BITS));
The corresponding operations for the col_matrix *cm are
!!(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));
Although the division / and modulus (or remainder) % are generally slow-ish (compared to addition, subtraction, and even multiplication), here WORD_BITS will be a power of two compile-time constant on all widely used architectures. All compilers I know of will turn the above into fast bit shifts and binary-AND operators.
To add row srcrow to row dstrow, you simply do a binary exclusive-or on all the words:
{
const word *const src = rm->row[srcrow];
word *const dst = rm->row[dstrow];
int w = rm->words;
while (w-->0)
dst[w] ^= src[w];
}
and analogously, for the column matrix,
{
const word *const src = cm->col[srccol];
word *const dst = cm->col[dstcol];
int w = cm->words;
while (w-->0)
dst[w] ^= src[w];
}
Note that if you combine more than two rows, you can do so very efficiently; it'll be MUCH faster than doing the additions consecutively. Intel and AMD CPUs are very good at predicting the above pattern, so you can just use more than one source row/column. Also, the destination does not have to participate in the result, although if I guess correctly what algorithm you're implementing, I guess you want it to.
If you know the target architectures have SSE2 or better, or even AVX, you can use the emmintrin.h or immintrin.h header files, respectively, for compiler built-in types and operators that allow you to XOR 128 bits and 256 bits, respectively, at once; sometimes giving you quite a bit of a boost.
Since the vector types require what the C standards call "excess alignment", you'll then also need to include mm_malloc.h, and use _mm_malloc() and _mm_free() to allocate the row/column vectors for the word data -- and obviously round words up so you can access the row/column as a suitable integer word type (__m128i for SSE*, __m256i for AVX).
Personally, I always implement the unvectorized version first, then some "nasty" test cases for unit testing, and only then see about vectorizing it. This has the benefit that you can give the unvectorized version as a preliminary version for those who will be using it, and you can compare the test case results between the vectorized and unvectorized cases, to see if one or the other has a bug in it.
The transpose operation is quite straightforward, although I do recommend using a triple loop: innermost looping over the bits in a single word. Also, you'll probably want to check which order -- row or column major -- works best for the outer loop; depending on the matrix size, you may see a huge difference. (This is due to cache behaviour: you want the CPU to be able to predict the access pattern, and not have to reload the same cachelines. In the best case, on at-most-few-years-old AMD and Intel x86-64 processors, you can get near cache speeds if both matrixes fit in the cache.)
All of the above can be implemented in a single header file -- even including the vectorized versions if the target architecture supports SSE2/AVX --, so it should not be too difficult to implement.
Questions?
