source: CIVL/examples/omp/m4ri/tests/test_invert.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;
  }
}