96 lines
1.9 KiB
C
96 lines
1.9 KiB
C
|
#include <m4ri/config.h>
|
||
|
|
#include <stdlib.h>
|
||
|
|
#include <m4ri/m4ri.h>
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Check that inversion works.
|
||
|
|
*
|
||
|
|
* \param n Number of rows of A
|
||
|
|
* \param k Parameter k of M4RM algorithm, may be 0 for automatic choice.
|
||
|
|
*/
|
||
|
|
int invert_test(rci_t n, int k) {
|
||
|
|
int ret = 0;
|
||
|
|
printf("invert: n: %4d, k: %2d", n, k);
|
||
|
|
|
||
|
|
mzd_t *I2 = mzd_init(n, n);
|
||
|
|
mzd_set_ui(I2, 1);
|
||
|
|
|
||
|
|
mzd_t *U = mzd_init(n,n);
|
||
|
|
mzd_randomize(U);
|
||
|
|
for(rci_t i=0; i<n; i++) {
|
||
|
|
mzd_write_bit(U, i, i, 1);
|
||
|
|
for (rci_t j=0; j<i; j++)
|
||
|
|
mzd_write_bit(U, i, j, 0);
|
||
|
|
}
|
||
|
|
|
||
|
|
mzd_t *B = mzd_copy(NULL, U);
|
||
|
|
mzd_trtri_upper(B);
|
||
|
|
|
||
|
|
mzd_t *I1 = mzd_mul(NULL, U, B, 0);
|
||
|
|
|
||
|
|
if (mzd_equal(I1, I2) != TRUE) {
|
||
|
|
ret += 1;
|
||
|
|
printf(" U*~U != 1 ");
|
||
|
|
}
|
||
|
|
|
||
|
|
mzd_t *L = mzd_init(n, n);
|
||
|
|
mzd_randomize(L);
|
||
|
|
for (rci_t i = 0; i < n; ++i) {
|
||
|
|
for (rci_t j = i + 1; j < n; ++j)
|
||
|
|
mzd_write_bit(L,i,j, 0);
|
||
|
|
mzd_write_bit(L,i,i, 1);
|
||
|
|
}
|
||
|
|
mzd_t *A = mzd_mul(NULL, L, U, 0);
|
||
|
|
|
||
|
|
B = mzd_inv_m4ri(B, A, k);
|
||
|
|
|
||
|
|
I1 = mzd_mul(I1, A, B, 0);
|
||
|
|
|
||
|
|
if (mzd_equal(I1, I2) != TRUE) {
|
||
|
|
ret += 1;
|
||
|
|
printf(" A*~A != 1 ");
|
||
|
|
}
|
||
|
|
|
||
|
|
if(ret == 0) {
|
||
|
|
printf(" ... passed\n");
|
||
|
|
} else {
|
||
|
|
printf(" ... FAILED\n");
|
||
|
|
}
|
||
|
|
mzd_free(U);
|
||
|
|
mzd_free(L);
|
||
|
|
mzd_free(A);
|
||
|
|
mzd_free(B);
|
||
|
|
mzd_free(I1);
|
||
|
|
mzd_free(I2);
|
||
|
|
|
||
|
|
return ret;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
int main() {
|
||
|
|
int status = 0;
|
||
|
|
srandom(17);
|
||
|
|
|
||
|
|
for(int k=0; k<5; k++) {
|
||
|
|
status += invert_test( 1,k);
|
||
|
|
status += invert_test( 2,k);
|
||
|
|
status += invert_test( 3,k);
|
||
|
|
status += invert_test( 21,k);
|
||
|
|
status += invert_test( 64,k);
|
||
|
|
status += invert_test( 128,k);
|
||
|
|
status += invert_test( 193,k);
|
||
|
|
status += invert_test(1000,k);
|
||
|
|
status += invert_test(1024,k);
|
||
|
|
status += invert_test(1025,k);
|
||
|
|
status += invert_test(1290,k);
|
||
|
|
status += invert_test(1710,k);
|
||
|
|
status += invert_test(2048,k);
|
||
|
|
}
|
||
|
|
if (status == 0) {
|
||
|
|
printf("All tests passed.\n");
|
||
|
|
return 0;
|
||
|
|
} else {
|
||
|
|
return -1;
|
||
|
|
}
|
||
|
|
}
|