M4RI 1.0.1
Data Structures | Defines | Functions
packedmatrix.h File Reference

Dense matrices over GF(2) represented as a bit field. More...

#include <math.h>
#include <assert.h>
#include <stdio.h>

Go to the source code of this file.

Data Structures

struct  mzd_t
 Dense matrices over GF(2). More...

Defines

#define MZD_MUL_BLOCKSIZE   MIN(((int)sqrt((double)(4*CPU_L2_CACHE)))/2,2048)
 Matrix multiplication block-ing dimension.
#define mzd_free_window   mzd_free
 Free a matrix window created with mzd_init_window.
#define mzd_sub   mzd_add
 Same as mzd_add.
#define _mzd_sub   _mzd_add
 Same as mzd_sub but without any checks on the input.

Functions

mzd_tmzd_init (const size_t r, const size_t c)
 Create a new matrix of dimension r x c.
void mzd_free (mzd_t *A)
 Free a matrix created with mzd_init.
mzd_tmzd_init_window (const mzd_t *M, const size_t lowr, const size_t lowc, const size_t highr, const size_t highc)
 Create a window/view into the matrix M.
static void mzd_row_swap (mzd_t *M, const size_t rowa, const size_t rowb)
 Swap the two rows rowa and rowb.
void mzd_copy_row (mzd_t *B, size_t i, const mzd_t *A, size_t j)
 copy row j from A to row i from B.
void mzd_col_swap (mzd_t *M, const size_t cola, const size_t colb)
 Swap the two columns cola and colb.
static void mzd_col_swap_in_rows (mzd_t *M, const size_t cola, const size_t colb, const size_t start_row, const size_t stop_row)
 Swap the two columns cola and colb but only between start_row and stop_row.
static BIT mzd_read_bit (const mzd_t *M, const size_t row, const size_t col)
 Read the bit at position M[row,col].
static void mzd_write_bit (mzd_t *M, const size_t row, const size_t col, const BIT value)
 Write the bit value to position M[row,col].
void mzd_print (const mzd_t *M)
 Print a matrix to stdout.
void mzd_print_tight (const mzd_t *M)
 Print the matrix to stdout.
static void mzd_row_add_offset (mzd_t *M, size_t dstrow, size_t srcrow, size_t coloffset)
 Add the rows sourcerow and destrow and stores the total in the row destrow, but only begins at the column coloffset.
void mzd_row_add (mzd_t *M, const size_t sourcerow, const size_t destrow)
 Add the rows sourcerow and destrow and stores the total in the row destrow.
mzd_tmzd_transpose (mzd_t *DST, const mzd_t *A)
 Transpose a matrix.
mzd_tmzd_mul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Naive cubic matrix multiplication.
mzd_tmzd_addmul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Naive cubic matrix multiplication and addition.
mzd_t_mzd_mul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B, const int clear)
 Naive cubic matrix multiplication with the pre-transposed B.
mzd_t_mzd_mul_va (mzd_t *C, const mzd_t *v, const mzd_t *A, const int clear)
 Matrix multiplication optimized for v*A where v is a vector.
void mzd_randomize (mzd_t *M)
 Fill matrix M with uniformly distributed bits.
void mzd_set_ui (mzd_t *M, const unsigned value)
 Set the matrix M to the value equivalent to the integer value provided.
int mzd_gauss_delayed (mzd_t *M, const size_t startcol, const int full)
 Gaussian elimination.
int mzd_echelonize_naive (mzd_t *M, const int full)
 Gaussian elimination.
BIT mzd_equal (const mzd_t *A, const mzd_t *B)
 Return TRUE if A == B.
int mzd_cmp (const mzd_t *A, const mzd_t *B)
 Return -1,0,1 if if A < B, A == B or A > B respectively.
mzd_tmzd_copy (mzd_t *DST, const mzd_t *A)
 Copy matrix A to DST.
mzd_tmzd_concat (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Concatenate B to A and write the result to C.
mzd_tmzd_stack (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Stack A on top of B and write the result to C.
mzd_tmzd_submatrix (mzd_t *S, const mzd_t *M, const size_t lowr, const size_t lowc, const size_t highr, const size_t highc)
 Copy a submatrix.
mzd_tmzd_invert_naive (mzd_t *INV, mzd_t *A, const mzd_t *I)
 Invert the matrix target using Gaussian elimination.
mzd_tmzd_add (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Set C = A+B.
mzd_t_mzd_add (mzd_t *C, const mzd_t *A, const mzd_t *B)
 Same as mzd_add but without any checks on the input.
void mzd_combine (mzd_t *DST, const size_t row3, const size_t startblock3, const mzd_t *SC1, const size_t row1, const size_t startblock1, const mzd_t *SC2, const size_t row2, const size_t startblock2)
 row3[col3:] = row1[col1:] + row2[col2:]
static word mzd_read_bits (const mzd_t *M, const size_t x, const size_t y, const int n)
static void mzd_xor_bits (const mzd_t *M, const size_t x, const size_t y, const int n, word values)
 XOR n bits from values to M starting a position (x,y).
static void mzd_and_bits (const mzd_t *M, const size_t x, const size_t y, const int n, word values)
 AND n bits from values to M starting a position (x,y).
static void mzd_clear_bits (const mzd_t *M, const size_t x, const size_t y, const int n)
 Clear n bits in M starting a position (x,y).
int mzd_is_zero (mzd_t *A)
 Zero test for matrix.
void mzd_row_clear_offset (mzd_t *M, const size_t row, const size_t coloffset)
 Clear the given row, but only begins at the column coloffset.
int mzd_find_pivot (mzd_t *M, size_t start_row, size_t start_col, size_t *r, size_t *c)
 Find the next nonzero entry in M starting at start_row and start_col.
double mzd_density (mzd_t *A, int res)
 Return the number of nonzero entries divided by nrows * ncols.
double _mzd_density (mzd_t *A, int res, size_t r, size_t c)
 Return the number of nonzero entries divided by nrows * ncols considering only the submatrix starting at (r,c).
size_t mzd_first_zero_row (mzd_t *A)
 Return the first row with all zero entries.

Detailed Description

Dense matrices over GF(2) represented as a bit field.

Author:
Gregory Bard <bard@fordham.edu>
Martin Albrecht <M.R.Albrecht@rhul.ac.uk>

Define Documentation

#define _mzd_sub   _mzd_add

Same as mzd_sub but without any checks on the input.

Parameters:
CPreallocated difference matrix, may be NULL for automatic creation.
AMatrix
BMatrix
Ignores offset atrtribute of packedmatrix.
#define mzd_free_window   mzd_free

Free a matrix window created with mzd_init_window.

Parameters:
AMatrix
#define MZD_MUL_BLOCKSIZE   MIN(((int)sqrt((double)(4*CPU_L2_CACHE)))/2,2048)

Matrix multiplication block-ing dimension.

Defines the number of rows of the matrix A that are processed as one block during the execution of a multiplication algorithm.

#define mzd_sub   mzd_add

Same as mzd_add.

Parameters:
CPreallocated difference matrix, may be NULL for automatic creation.
AMatrix
BMatrix
Ignores offset atrtribute of packedmatrix.

Function Documentation

mzd_t* _mzd_add ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Same as mzd_add but without any checks on the input.

Parameters:
CPreallocated sum matrix, may be NULL for automatic creation.
AMatrix
BMatrix
Ignores offset atrtribute of packedmatrix.
double _mzd_density ( mzd_t A,
int  res,
size_t  r,
size_t  c 
)

Return the number of nonzero entries divided by nrows * ncols considering only the submatrix starting at (r,c).

If res = 0 then 100 samples per row are made, if res > 0 the function takes res sized steps within each row (res = 1 uses every word).

Parameters:
AMatrix
resResolution of sampling
rRow to start counting
cColumn to start counting
mzd_t* _mzd_mul_naive ( mzd_t C,
const mzd_t A,
const mzd_t B,
const int  clear 
)

Naive cubic matrix multiplication with the pre-transposed B.

That is, compute C such that C == AB^t.

Parameters:
CPreallocated product matrix.
AInput matrix A.
BPre-transposed input matrix B.
clearWhether to clear C before accumulating AB
mzd_t* _mzd_mul_va ( mzd_t C,
const mzd_t v,
const mzd_t A,
const int  clear 
)

Matrix multiplication optimized for v*A where v is a vector.

Parameters:
CPreallocated product matrix.
vInput matrix v.
AInput matrix A.
clearIf set clear C first, otherwise add result to C.
mzd_t* mzd_add ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Set C = A+B.

C is also returned. If C is NULL then a new matrix is created which must be freed by mzd_free.

Parameters:
CPreallocated sum matrix, may be NULL for automatic creation.
AMatrix
BMatrix
Examples:
testsuite/test_multiplication.c.
mzd_t* mzd_addmul_naive ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Naive cubic matrix multiplication and addition.

That is, compute C such that C == C + AB.

Parameters:
CPreallocated product matrix.
AInput matrix A.
BInput matrix B.
Note:
Normally, if you will multiply several times by b, it is smarter to calculate bT yourself, and keep it, and then use the function called _mzd_mul_naive
static void mzd_and_bits ( const mzd_t M,
const size_t  x,
const size_t  y,
const int  n,
word  values 
) [inline, static]

AND n bits from values to M starting a position (x,y).

Parameters:
MSource matrix.
xStarting row.
yStarting column.
nNumber of bits (<= RADIX);
valuesWord with values;
static void mzd_clear_bits ( const mzd_t M,
const size_t  x,
const size_t  y,
const int  n 
) [inline, static]

Clear n bits in M starting a position (x,y).

Parameters:
MSource matrix.
xStarting row.
yStarting column.
nNumber of bits (<= RADIX);
int mzd_cmp ( const mzd_t A,
const mzd_t B 
)

Return -1,0,1 if if A < B, A == B or A > B respectively.

Parameters:
AMatrix.
BMatrix.
Note:
This comparison is not well defined mathematically and relatively arbitrary since elements of GF(2) don't have an ordering.
Ignores offset atrtribute of packedmatrix.
void mzd_col_swap ( mzd_t M,
const size_t  cola,
const size_t  colb 
)

Swap the two columns cola and colb.

Parameters:
MMatrix.
colaColumn index.
colbColumn index.
static void mzd_col_swap_in_rows ( mzd_t M,
const size_t  cola,
const size_t  colb,
const size_t  start_row,
const size_t  stop_row 
) [inline, static]

Swap the two columns cola and colb but only between start_row and stop_row.

Parameters:
MMatrix.
colaColumn index.
colbColumn index.
start_rowRow index.
stop_rowRow index (exclusive).
void mzd_combine ( mzd_t DST,
const size_t  row3,
const size_t  startblock3,
const mzd_t SC1,
const size_t  row1,
const size_t  startblock1,
const mzd_t SC2,
const size_t  row2,
const size_t  startblock2 
)

row3[col3:] = row1[col1:] + row2[col2:]

Adds row1 of SC1, starting with startblock1 to the end, to row2 of SC2, starting with startblock2 to the end. This gets stored in DST, in row3, starting with startblock3.

Parameters:
DSTdestination matrix
row3destination row for matrix dst
startblock3starting block to work on in matrix dst
SC1source matrix
row1source row for matrix sc1
startblock1starting block to work on in matrix sc1
SC2source matrix
startblock2starting block to work on in matrix sc2
row2source row for matrix sc2
Ignores offset atrtribute of packedmatrix.

Todo:
this code is slow if offset!=0

mzd_t* mzd_concat ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Concatenate B to A and write the result to C.

That is,

 [ A ], [ B ] -> [ A  B ] = C
 

The inputs are not modified but a new matrix is created.

Parameters:
CMatrix, may be NULL for automatic creation
AMatrix
BMatrix
Note:
This is sometimes called augment.
Ignores offset atrtribute of packedmatrix.
mzd_t* mzd_copy ( mzd_t DST,
const mzd_t A 
)

Copy matrix A to DST.

Parameters:
DSTMay be NULL for automatic creation.
ASource matrix.
Ignores offset atrtribute of packedmatrix.
Examples:
testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.
void mzd_copy_row ( mzd_t B,
size_t  i,
const mzd_t A,
size_t  j 
)

copy row j from A to row i from B.

The offsets of A and B must match and the number of columns of A must be less than or equal to the number of columns of B.

Parameters:
BTarget matrix.
iTarget row index.
ASource matrix.
jSource row index.
double mzd_density ( mzd_t A,
int  res 
)

Return the number of nonzero entries divided by nrows * ncols.

If res = 0 then 100 samples per row are made, if res > 0 the function takes res sized steps within each row (res = 1 uses every word).

Parameters:
AMatrix
resResolution of sampling
int mzd_echelonize_naive ( mzd_t M,
const int  full 
)

Gaussian elimination.

This will do Gaussian elimination on the matrix m. If full=FALSE, then it will do triangular style elimination, and if full=TRUE, it will do Gauss-Jordan style, or full elimination.

Parameters:
MMatrix
fullGauss-Jordan style or upper triangular form only.
Ignores offset atrtribute of packedmatrix.
See also:
mzd_echelonize_m4ri(), mzd_echelonize_pluq()
Examples:
testsuite/bench_elimination.c, and testsuite/test_elimination.c.
BIT mzd_equal ( const mzd_t A,
const mzd_t B 
)

Return TRUE if A == B.

Parameters:
AMatrix
BMatrix
Ignores offset atrtribute of packedmatrix.
Examples:
testsuite/test_elimination.c, and testsuite/test_multiplication.c.
int mzd_find_pivot ( mzd_t M,
size_t  start_row,
size_t  start_col,
size_t *  r,
size_t *  c 
)

Find the next nonzero entry in M starting at start_row and start_col.

This function walks down rows in the inner loop and columns in the outer loop. If a nonzero entry is found this function returns 1 and zero otherwise.

If and only if a nonzero entry is found r and c are updated.

Parameters:
MMatrix
start_rowIndex of row where to start search
start_colIndex of column where to start search
rRow index updated if pivot is found
cColumn index updated if pivot is found
size_t mzd_first_zero_row ( mzd_t A)

Return the first row with all zero entries.

If no such row can be found returns nrows.

Parameters:
AMatrix
void mzd_free ( mzd_t A)

Free a matrix created with mzd_init.

Parameters:
AMatrix
Examples:
testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.
int mzd_gauss_delayed ( mzd_t M,
const size_t  startcol,
const int  full 
)

Gaussian elimination.

This will do Gaussian elimination on the matrix m but will start not at column 0 necc but at column startcol. If full=FALSE, then it will do triangular style elimination, and if full=TRUE, it will do Gauss-Jordan style, or full elimination.

Parameters:
MMatrix
startcolFirst column to consider for reduction.
fullGauss-Jordan style or upper triangular form only.
Ignores offset atrtribute of packedmatrix.
mzd_t* mzd_init ( const size_t  r,
const size_t  c 
)

Create a new matrix of dimension r x c.

Use mzd_free to kill it.

Parameters:
rNumber of rows
cNumber of columns
Examples:
testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.
mzd_t* mzd_init_window ( const mzd_t M,
const size_t  lowr,
const size_t  lowc,
const size_t  highr,
const size_t  highc 
)

Create a window/view into the matrix M.

A matrix window for M is a meta structure on the matrix M. It is setup to point into the matrix so M must not be freed while the matrix window is used.

This function puts the restriction on the provided parameters that all parameters must be within range for M which is not enforced currently .

Use mzd_free_window to free the window.

Parameters:
MMatrix
lowrStarting row (inclusive)
lowcStarting column (inclusive)
highrEnd row (exclusive)
highcEnd column (exclusive)
mzd_t* mzd_invert_naive ( mzd_t INV,
mzd_t A,
const mzd_t I 
)

Invert the matrix target using Gaussian elimination.

To avoid recomputing the identity matrix over and over again, I may be passed in as identity parameter.

Parameters:
INVPreallocated space for inversion matrix, may be NULL for automatic creation.
AMatrix to be reduced.
IIdentity matrix.
Ignores offset atrtribute of packedmatrix.
int mzd_is_zero ( mzd_t A)

Zero test for matrix.

Parameters:
AInput matrix.
mzd_t* mzd_mul_naive ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Naive cubic matrix multiplication.

That is, compute C such that C == AB.

Parameters:
CPreallocated product matrix, may be NULL for automatic creation.
AInput matrix A.
BInput matrix B.
Note:
Normally, if you will multiply several times by b, it is smarter to calculate bT yourself, and keep it, and then use the function called _mzd_mul_naive
Examples:
testsuite/test_multiplication.c.
void mzd_print ( const mzd_t M)

Print a matrix to stdout.

The output will contain colons between every 4-th column.

Parameters:
MMatrix
void mzd_print_tight ( const mzd_t M)

Print the matrix to stdout.

Parameters:
MMatrix
void mzd_randomize ( mzd_t M)

Fill matrix M with uniformly distributed bits.

Parameters:
MMatrix
Todo:
Allow the user to provide a RNG callback.
Ignores offset atrtribute of packedmatrix.
Examples:
testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.
static BIT mzd_read_bit ( const mzd_t M,
const size_t  row,
const size_t  col 
) [inline, static]

Read the bit at position M[row,col].

Parameters:
MMatrix
rowRow index
colColumn index
Note:
No bounds checks whatsoever are performed.
Examples:
testsuite/test_pluq.c.
static word mzd_read_bits ( const mzd_t M,
const size_t  x,
const size_t  y,
const int  n 
) [inline, static]

Get n bits starting a position (x,y) from the matrix M.

Parameters:
MSource matrix.
xStarting row.
yStarting column.
nNumber of bits (<= RADIX);
void mzd_row_add ( mzd_t M,
const size_t  sourcerow,
const size_t  destrow 
)

Add the rows sourcerow and destrow and stores the total in the row destrow.

Parameters:
MMatrix
sourcerowIndex of source row
destrowIndex of target row
Note:
this can be done much faster with mzd_combine.
static void mzd_row_add_offset ( mzd_t M,
size_t  dstrow,
size_t  srcrow,
size_t  coloffset 
) [inline, static]

Add the rows sourcerow and destrow and stores the total in the row destrow, but only begins at the column coloffset.

Parameters:
MMatrix
dstrowIndex of target row
srcrowIndex of source row
coloffsetColumn offset
void mzd_row_clear_offset ( mzd_t M,
const size_t  row,
const size_t  coloffset 
)

Clear the given row, but only begins at the column coloffset.

Parameters:
MMatrix
rowIndex of row
coloffsetColumn offset
static void mzd_row_swap ( mzd_t M,
const size_t  rowa,
const size_t  rowb 
) [inline, static]

Swap the two rows rowa and rowb.

Parameters:
MMatrix
rowaRow index.
rowbRow index.
void mzd_set_ui ( mzd_t M,
const unsigned  value 
)

Set the matrix M to the value equivalent to the integer value provided.

Specifically, this function does nothing if value%2 == 0 and returns the identity matrix if value%2 == 1.

If the matrix is not square then the largest possible square submatrix is set to the identity matrix.

Parameters:
MMatrix
valueEither 0 or 1
mzd_t* mzd_stack ( mzd_t C,
const mzd_t A,
const mzd_t B 
)

Stack A on top of B and write the result to C.

That is,

 [ A ], [ B ] -> [ A ] = C
                 [ B ]
 

The inputs are not modified but a new matrix is created.

Parameters:
CMatrix, may be NULL for automatic creation
AMatrix
BMatrix
Ignores offset atrtribute of packedmatrix.
mzd_t* mzd_submatrix ( mzd_t S,
const mzd_t M,
const size_t  lowr,
const size_t  lowc,
const size_t  highr,
const size_t  highc 
)

Copy a submatrix.

Note that the upper bounds are not included.

Parameters:
SPreallocated space for submatrix, may be NULL for automatic creation.
MMatrix
lowrstart rows
lowcstart column
highrstop row (this row is not included)
highcstop column (this column is not included)
mzd_t* mzd_transpose ( mzd_t DST,
const mzd_t A 
)

Transpose a matrix.

This function uses the fact that:

   [ A B ]T    [AT CT]
   [ C D ]  =  [BT DT] 
 

and thus rearranges the blocks recursively.

Parameters:
DSTPreallocated return matrix, may be NULL for automatic creation.
AMatrix
static void mzd_write_bit ( mzd_t M,
const size_t  row,
const size_t  col,
const BIT  value 
) [inline, static]

Write the bit value to position M[row,col].

Parameters:
MMatrix
rowRow index
colColumn index
valueEither 0 or 1
Note:
No bounds checks whatsoever are performed.
Examples:
testsuite/bench_elimination.c, and testsuite/test_pluq.c.
static void mzd_xor_bits ( const mzd_t M,
const size_t  x,
const size_t  y,
const int  n,
word  values 
) [inline, static]

XOR n bits from values to M starting a position (x,y).

Parameters:
MSource matrix.
xStarting row.
yStarting column.
nNumber of bits (<= RADIX);
valuesWord with values;