1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
//! Links to mzd.h
//!
//! FIXME implement missing functions
//!
use libc;
use std::mem::size_of;

use crate::misc::m4ri_one;
use crate::misc::m4ri_radix;
use crate::misc::Rci;
use crate::misc::Wi;
use crate::misc::Word;
use crate::misc::BIT;

/// Represents the blocks used by M4RI internally
#[repr(C)]
struct MzdBlock {
    size: libc::size_t,
    begin: *mut Word,
    end: *mut Word,
}

/// Represents the Mzd data type used by M4RI
#[repr(C)]
pub struct Mzd {
    /// Number of rows
    pub nrows: Rci,
    /// Number of columns
    pub ncols: Rci,
    /// Number of words with valid bits:
    /// width = ceil(ncols / m4ri_radix)
    pub width: Wi,

    /// Offset in words between rows
    /// ``
    /// rowstride = (width < mzd_paddingwidth || (width & 1) == 0) ? width : width + 1;
    /// ``
    /// where width is the width of the underlying non-windowed matrix
    rowstride: Wi,

    /// Offset in words from start of block to first word
    ///
    /// ``rows[0] = blocks[0].begin + offset_vector``
    offset_vector: Wi,

    /// Number of rows to the first row counting from the start of the
    /// first block
    row_offset: Wi,

    /// Booleans to speed up things
    ///
    /// The bits have the following meaning
    ///
    /// 1: Has non-zero excess
    /// 2. Is windowed, but has zero offset
    /// 3. Is windowed, but has zero excess
    /// 4. Is windowed, but owns the Blocks allocations
    /// 5. Spans more than 1 Block
    flags: u8,

    /// blockrows_log = log2(blockrows)
    /// where blockrows is the number of rows in one block,
    /// which is a power of 2.
    blockrows_log: u8,

    // Ensures sizeof(mzd_t) == 64
    padding: [u8; 62
        - 2 * size_of::<Rci>()
        - 4 * size_of::<Wi>()
        - size_of::<Word>()
        - 2 * size_of::<*const libc::c_void>()],

    /// Mask for valid bits in the word with the highest index (width - 1)
    high_bitmask: Word,
    /// Pointers to the actual blocks of memory containing the values packed into words
    blocks: *const MzdBlock,
    /// Address of first word in each row, so the first word of `row [i]` is in `m->rows[i]`
    pub rows: *const *mut Word,
}

/// Flag when `ncols%64 == 0`
pub static MZD_FLAG_NONZERO_EXCESS: u8 = 0x2;

/// Flag for windowed matrix
pub static MZD_FLAG_WINDOWED_ZEROOFFSET: u8 = 0x4;

/// Flag for windowed matrix where `ncols % 64 == 0`
pub static MZD_FLAG_WINDOWED_ZEROEXCESS: u8 = 0x8;

/// Flag for windowed matrix which owns its memory
pub static MZD_FLAG_WINDOWED_OWNSBLOCKS: u8 = 0x10;

/// Flag for multiple blocks
pub static MZD_FLAG_MULTIPLE_BLOCKS: u8 = 0x20;

extern "C" {
    /// Create a new rows x columns matrix
    pub fn mzd_init(rows: Rci, columns: Rci) -> *mut Mzd;

    /// Free a matrix created with mzd_init.
    /// Automatically done by the Deref trait on Mzd
    pub fn mzd_free(matrix: *mut Mzd);

    /// \brief Create a window/view into the matrix M.
    ///
    /// A matrix window for M is a meta structure on the matrix M. It i
    /// setup to point into the matrix so M \em must \em 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.
    ///
    /// \param M Matrix
    /// \param lowr Starting row (inclusive)
    /// \param lowc Starting column (inclusive, **must be multiple of m4ri_radix**)
    /// \param highr End row (exclusive)
    /// \param highc End column (exclusive)
    pub fn mzd_init_window(
        matrix: *mut Mzd,
        lowr: Rci,
        lowc: Rci,
        highr: Rci,
        highc: Rci,
    ) -> *mut Mzd;

    /// Swap the two rows rowa and rowb
    pub fn mzd_row_swap(matrix: *mut Mzd, rowa: Rci, rowb: Rci);

    /// \brief 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.
    ///
    /// \param B Target matrix.
    /// \param i Target row index.
    /// \param A Source matrix.
    /// \param j Source row index.
    pub fn mzd_copy_row(b: *mut Mzd, i: Rci, a: *const Mzd, j: Rci);

    /// Swap the two columns cola and colb
    pub fn mzd_col_swap(matrix: *mut Mzd, cola: Rci, colb: Rci);

    /// Transpose a matrix
    /// Dest may be null for automatic creation
    pub fn mzd_transpose(dest: *mut Mzd, source: *const Mzd) -> *mut Mzd;

    /// naive cubic matrix multiplication
    /// the first argument may be null for automatic creation
    pub fn mzd_mul_naive(dest: *mut Mzd, a: *const Mzd, b: *const Mzd) -> *mut Mzd;

    ///
    /// brief Naive cubic matrix multiplication with the pre-transposed B.
    ///
    /// That is, compute C such that C == AB^t.
    ///
    /// param C Preallocated product matrix.
    /// param A Input matrix A.
    /// param B Pre-transposed input matrix B.
    /// param clear Whether to clear C before accumulating AB
    pub fn _mzd_mul_naive(
        dest: *mut Mzd,
        a: *const Mzd,
        b: *const Mzd,
        clear: libc::c_int,
    ) -> *mut Mzd;

    /// naive cubic matrix multiplication and addition
    ///
    /// C == C + AB
    pub fn mzd_addmul_naive(c: *mut Mzd, a: *const Mzd, b: *const Mzd) -> *mut Mzd;

    /// Matrix multiplication optimized for v*A where v is a vector
    ///
    /// param C: preallocated product matrix
    /// param v: input matrix v
    /// param A: input matrix A
    /// param clear: if set clear C first, otherwise add result to C
    pub fn _mzd_mul_va(c: *mut Mzd, v: *const Mzd, a: *const Mzd, clear: libc::c_int) -> *mut Mzd;

    /// Fill the matrix m with uniformly distributed bits.
    pub fn mzd_randomize(m: *mut Mzd);

    /// Return true if A == B
    pub fn mzd_equal(a: *const Mzd, b: *const Mzd) -> libc::c_int;

    /// Copy a matrix to dest
    ///
    /// Dest may be null for automatic creation
    pub fn mzd_copy(dest: *mut Mzd, src: *const Mzd) -> *mut Mzd;

    /// Concatenate B to A and write the result to C
    pub fn mzd_concat(c: *mut Mzd, a: *const Mzd, b: *const Mzd) -> *mut Mzd;

    /// Set to identity matrix if the second argument is 1
    pub fn mzd_set_ui(a: *mut Mzd, n: libc::c_uint);

    /// Stack A on top of B into C
    pub fn mzd_stack(c: *mut Mzd, a: *mut Mzd, b: *const Mzd) -> *mut Mzd;

    /// Copy a submatrix
    /// first argument may be preallocated space or null
    pub fn mzd_submatrix(
        s: *mut Mzd,
        m: *const Mzd,
        lowr: Rci,
        lowc: Rci,
        highr: Rci,
        highc: Rci,
    ) -> *mut Mzd;

    /// Invert the target matrix using gaussian elimination
    /// To avoid recomputing the identity matrix over and over again,
    /// I may be passed in as identity parameter
    /// The first parameter may be null to have the space automatically allocated
    pub fn mzd_invert_naive(inv: *mut Mzd, a: *const Mzd, identity: *const Mzd) -> *mut Mzd;

    /// Set C = A + B
    /// If C is passed in, the result is written there
    /// otherwise a new matrix is created
    pub fn mzd_add(c: *mut Mzd, a: *const Mzd, b: *const Mzd) -> *mut Mzd;

    /// Set C = A - B
    /// If C is passed in, the result is written there
    /// otherwise a new matrix is created
    ///
    /// Secretly an alias for mzd_add
    pub fn mzd_sub(c: *mut Mzd, a: *const Mzd, b: *const Mzd) -> *mut Mzd;

    /// Zero test for matrix
    pub fn mzd_is_zero(a: *const Mzd);

    /// Clear the given row, but only begins at the column coloffset.
    ///
    /// param M Matrix
    /// param row Index of row
    /// param coloffset Column offset
    pub fn mzd_row_clear_offset(m: *mut Mzd, row: Rci, coloffset: Rci);

}

/// Write the bit to position M[row, col]
#[inline]
pub unsafe fn mzd_write_bit(matrix: *mut Mzd, row: Rci, col: Rci, value: BIT) {
    let therow: *const *mut Word = (*matrix).rows.offset(row as isize);
    let column: *mut Word = (*therow).offset((col / m4ri_radix) as isize);
    let pos = col % m4ri_radix;
    let column_bitmasked: Word = *column & !(m4ri_one << pos);
    let column_newbit: Word = (value as Word & m4ri_one) << pos;
    debug_assert_eq!(column_newbit.count_ones(), value as u32);
    *column = column_bitmasked | column_newbit;
}

/// Read the bit at position M[row, col]
///
/// # Unsafe behaviour
/// No bounds checking
///
/// Reimplemented in Rust as the C library declares it as inline
#[inline]
pub unsafe fn mzd_read_bit(matrix: *const Mzd, row: Rci, col: Rci) -> BIT {
    let therow: *const *mut Word = (*matrix).rows.offset(row as isize);
    let column: Word = *(*therow).offset((col / m4ri_radix) as isize);

    let thebit = (column >> (col % m4ri_radix)) & m4ri_one;
    thebit as BIT
}

/// Get a pointer to the first word
///
/// param: M: matrix
///
/// Return a pointer to the first word of the first row
#[inline]
pub unsafe fn mzd_first_row(matrix: *const Mzd) -> *mut Word {
    let result: *mut Word = (*(*matrix).blocks)
        .begin
        .offset((*matrix).offset_vector as isize);
    debug_assert!(
        (*matrix).nrows == 0 || result == *(*matrix).rows,
        "Result is not the expected ptr"
    );
    result
}

/// Get pointer to first word of row
///
/// Param M Matrix
/// Param row the row index
#[inline]
pub unsafe fn mzd_row(matrix: *const Mzd, row: Rci) -> *mut Word {
    debug_assert!(row >= 0);
    let big_vector: Wi = (*matrix).offset_vector + row * (*matrix).rowstride;
    // FIXME __M4RI_UNLIKELY -> _builtin_expect
    let result: *mut Word = if (*matrix).flags & MZD_FLAG_MULTIPLE_BLOCKS != 0 {
        let n = ((*matrix).row_offset + row) >> (*matrix).blockrows_log;
        (*(*matrix).blocks.add(n as usize)).begin.offset(
            (big_vector - n * ((*(*matrix).blocks).size / ::std::mem::size_of::<Word>()) as i32)
                as isize,
        )
    } else {
        (*(*matrix).blocks).begin.add(big_vector as usize)
    };

    debug_assert_eq!(
        result,
        *(*matrix).rows.add(row as usize),
        "Result is not the expected ptr"
    );
    result
}

/// Create a const window/view into a const matrix
///
/// Note that this function still allocates a new Mzd struct that needs to be dropped.
/// You **must** call ``mzd_free_window_const``.
///
/// Also, this function copies a bunch of references
/// As a result, you could get multiple mut references into the same memory.
/// This is very unsafe if not used properly.
pub unsafe fn mzd_init_window_const(
    matrix: *const Mzd,
    lowr: Rci,
    lowc: Rci,
    highr: Rci,
    highc: Rci,
) -> *const Mzd {
    mzd_init_window(std::mem::transmute(matrix), lowr, lowc, highr, highc)
}

/// Test if a matrix is windowed
///
/// return a non-zero value if the matrix is windowed, otherwise return zero
#[inline]
pub unsafe fn mzd_is_windowed(m: *const Mzd) -> u8 {
    (*m).flags & MZD_FLAG_WINDOWED_ZEROOFFSET
}

/// Test if this mzd_t should free blocks
#[inline]
pub unsafe fn mzd_owns_blocks(m: *const Mzd) -> bool {
    !(*m).blocks.is_null()
        && (mzd_is_windowed(m) == 0 || ((*m).flags & MZD_FLAG_WINDOWED_OWNSBLOCKS != 0))
}

impl Drop for Mzd {
    #[inline]
    fn drop(&mut self) {
        unsafe {
            mzd_free(self);
        }
    }
}

/// Free a matrix window created with mzd_init_window
///
/// This is actually just `mzd_free` so call `ptr::drop_in_place` instead
#[inline]
pub unsafe fn mzd_free_window(matrix: *mut Mzd) {
    std::ptr::drop_in_place(matrix)
}

/// Free a "const" window created with `mzd_init_window_const`.
///
/// This function *MUST* be called for const windows.
#[inline]
pub unsafe fn mzd_free_window_const(matrix: *const Mzd) {
    let matrix: *mut Mzd = std::mem::transmute(matrix);
    std::ptr::drop_in_place(matrix)
}

#[cfg(test)]
mod test {
    use super::*;
    use std::mem;
    use std::ptr;

    #[test]
    fn init() {
        for _ in 0..100 {
            let result: libc::c_int;
            unsafe {
                assert_eq!(mem::size_of::<Mzd>(), 64);
                let matrix = mzd_init(10, 10);
                assert!(!(*matrix).blocks.is_null());
                assert!(!(*matrix).rows.is_null());
                mzd_randomize(matrix);
                result = mzd_equal(matrix, matrix);
                assert_eq!(result, 1);
                let m2 = mzd_copy(ptr::null_mut(), matrix);
                mzd_randomize(m2);
                assert_eq!(mzd_equal(m2, matrix), 0);
                ptr::drop_in_place(matrix);
                ptr::drop_in_place(m2);
            }
        }
    }
    #[test]
    fn test_mzd_first_row() {
        for _ in 0..100 {
            unsafe {
                let matrix = mzd_init(10, 10);
                mzd_set_ui(matrix, 0);
                mzd_first_row(matrix);
                ptr::drop_in_place(matrix);
            }
        }
    }

    #[test]
    fn test_mzd_row() {
        for _ in 0..100 {
            unsafe {
                let matrix = mzd_init(10, 10);
                mzd_set_ui(matrix, 0);
                mzd_row(matrix, 5);
                ptr::drop_in_place(matrix);
            }
        }
    }

    #[test]
    fn test_mzd_read_bit() {
        for _ in 0..10 {
            unsafe {
                let matrix = mzd_init(1000, 1000);
                mzd_set_ui(matrix, 1);
                for i in 0..1000 {
                    for j in 0..1000 {
                        let bit = mzd_read_bit(matrix, i as Rci, j as Rci);
                        assert_eq!(bit == 1, i == j, "Should be unit matrix");
                    }
                }
                ptr::drop_in_place(matrix);
            }
        }
    }

    #[test]
    fn test_mzd_write_bit() {
        for _ in 0..10 {
            unsafe {
                let matrix = mzd_init(1000, 1000);
                for i in 0..1000 {
                    for j in 0..1000 {
                        mzd_write_bit(matrix, i as Rci, j as Rci, if i == j { 1 } else { 0 });
                    }
                }
                for i in 0..1000 {
                    for j in 0..1000 {
                        let bit = mzd_read_bit(matrix, i as Rci, j as Rci);
                        assert_eq!(bit == 1, i == j, "Should be unit matrix");
                    }
                }
                ptr::drop_in_place(matrix);
            }
        }
    }
}