M4RI 1.0.1
testsuite/test_pluq.c
#include <stdlib.h>
#include "m4ri/m4ri.h"

int test_pluq_full_rank (size_t m, size_t n){
  printf("pluq: testing full rank m: %5zu, n: %5zu",m,n);

  mzd_t* U = mzd_init (m,n);
  mzd_t* L = mzd_init (m,m);
  mzd_t* U2 = mzd_init (m,n);
  mzd_t* L2 = mzd_init (m,m);
  mzd_t* A = mzd_init (m,n);
  mzd_randomize (U);
  mzd_randomize (L);

  size_t i,j;
  for (i=0; i<m; ++i){
    for (j=0; j<i && j<n;++j)
      mzd_write_bit(U,i,j, 0);
    for (j=i+1; j<m;++j)
      mzd_write_bit(L,i,j, 0);
    if(i<n)
      mzd_write_bit(U,i,i, 1);
    mzd_write_bit(L,i,i, 1);
  }
  
  mzd_mul(A, L, U, 2048);

  mzd_t* Acopy = mzd_copy (NULL,A);

  mzp_t* P = mzp_init(m);
  mzp_t* Q = mzp_init(n);
  mzd_pluq(A, P, Q, 2048);

  for (i=0; i<m; ++i){
    for (j=0; j<i && j <n;++j)
      mzd_write_bit (L2, i, j, mzd_read_bit(A,i,j));
    for (j=i+1; j<n;++j)
      mzd_write_bit (U2, i, j, mzd_read_bit(A,i,j));
  }
  
  for (i=0; i<n && i<m; ++i){
    mzd_write_bit(L2,i,i, 1);
    mzd_write_bit(U2,i,i, 1);
  }
  mzd_addmul(Acopy,L2,U2,0);
  int status = 0;
  for ( i=0; i<m; ++i)
    for ( j=0; j<n; ++j){
      if (mzd_read_bit (Acopy,i,j)){
        status = 1;
      }
    }
  if (status){
    printf(" ... FAILED\n");
  }  else
    printf (" ... passed\n");
  mzd_free(U);
  mzd_free(L);
  mzd_free(U2);
  mzd_free(L2);
  mzd_free(A);
  mzd_free(Acopy);
  mzp_free(P);
  mzp_free(Q);
  return status;
}

int test_pluq_half_rank(size_t m, size_t n) {
  printf("pluq: testing half rank m: %5zd, n: %5zd",m,n);

  mzd_t* U = mzd_init(m, n);
  mzd_t* L = mzd_init(m, m);
  mzd_t* U2 = mzd_init(m, n);
  mzd_t* L2 = mzd_init(m, m);
  mzd_t* A = mzd_init(m, n);
  mzd_randomize (U);
  mzd_randomize (L);

  size_t i,j;
  for (i=0; i<m && i<n; ++i){
    mzd_write_bit(U,i,i, 1);
    for (j=0; j<i;++j)
      mzd_write_bit(U,i,j, 0);
    if (i%2)
      for (j=i; j<n;++j)
        mzd_write_bit(U,i,j, 0);
    for (j=i+1; j<m;++j)
      mzd_write_bit(L,i,j, 0);
    mzd_write_bit(L,i,i, 1);
  }
  
  mzd_mul(A, L, U, 0);

  mzd_t* Acopy = mzd_copy (NULL,A);



  mzp_t* Pt = mzp_init(m);
  mzp_t* Q = mzp_init(n);
  int r = mzd_pluq(A, Pt, Q, 0);

  for (i=0; i<r; ++i){
    for (j=0; j<i;++j)
      mzd_write_bit (L2, i, j, mzd_read_bit(A,i,j));
    for (j=i+1; j<n;++j)
      mzd_write_bit (U2, i, j, mzd_read_bit(A,i,j));
  }
  for (i=r; i<m; i++)
    for (j=0; j<r;++j)
      mzd_write_bit (L2, i, j, mzd_read_bit(A,i,j));
  for (i=0; i<r; ++i){
    mzd_write_bit(L2,i,i, 1);
    mzd_write_bit(U2,i,i, 1);
  }

  mzd_apply_p_left(Acopy, Pt);
  mzd_apply_p_right_trans(Acopy, Q);
  mzd_addmul(Acopy,L2,U2,0);

  int status = 0;
  for ( i=0; i<m; ++i) {
    for ( j=0; j<n; ++j){
      if (mzd_read_bit(Acopy,i,j)){
        status = 1;
      }
    }
    if(status)
      break;
  }
  if (status)
    printf(" ... FAILED\n");
  else
    printf (" ... passed\n");
  mzd_free(U);
  mzd_free(L);
  mzd_free(U2);
  mzd_free(L2);
  mzd_free(A);
  mzd_free(Acopy);
  mzp_free(Pt);
  mzp_free(Q);
  return status;
}

int test_pluq_structured(size_t m, size_t n) {

  printf("pluq: testing structured m: %5zd, n: %5zd", m, n);

  size_t i,j;
  mzd_t* A = mzd_init(m, n);
  mzd_t* L = mzd_init(m, m);
  mzd_t* U = mzd_init(m, n);

  for(i=0; i<m; i+=2)
    for (j=i; j<n; j++)
      mzd_write_bit(A, i, j, 1);

  mzd_t* Acopy = mzd_copy (NULL,A);

  mzp_t* P = mzp_init(m);
  mzp_t* Q = mzp_init(n);
  int r;
  r=mzd_pluq(A, P, Q, 0);
  printf(", rank: %5d ",r);

  for (i=0; i<r; ++i){
    for (j=0; j<i;++j)
      mzd_write_bit(L, i, j, mzd_read_bit(A,i,j));
    for (j=i+1; j<n;++j)
      mzd_write_bit(U, i, j, mzd_read_bit(A,i,j));
  }
  for (i=r; i<m; i++)
    for (j=0; j<r;++j)
      mzd_write_bit(L, i, j, mzd_read_bit(A,i,j));
  for (i=0; i<r; ++i){
    mzd_write_bit(L,i,i, 1);
    mzd_write_bit(U,i,i, 1);
  }

  mzd_apply_p_left(Acopy, P);
  mzd_apply_p_right_trans(Acopy, Q);

  mzd_addmul(Acopy, L, U, 0);
  int status = 0;
  for ( i=0; i<m; ++i)
    for ( j=0; j<n; ++j){
      if (mzd_read_bit (Acopy,i,j)){
        status = 1;
        break;
      }
    }

  if (status) {
    printf("\n");
    printf(" ... FAILED\n");
  }  else
    printf (" ... passed\n");
  mzd_free(U);
  mzd_free(L);
  mzd_free(A);
  mzd_free(Acopy);
  mzp_free(P);
  mzp_free(Q);
  return status;
}

int test_pluq_random(size_t m, size_t n) {
  printf("pluq: testing random m: %5zd, n: %5zd",m,n);

  size_t i,j;
  mzd_t* U = mzd_init(m, n);
  mzd_t* L = mzd_init(m, m);
  mzd_t* A = mzd_init(m, n);
  mzd_randomize(A);

  mzd_t* Acopy = mzd_copy (NULL,A);

  mzp_t* P = mzp_init(m);
  mzp_t* Q = mzp_init(n);
  int r;
  r=mzd_pluq(A, P, Q, 0);
  printf(", rank: %5d ",r);

  for (i=0; i<r; ++i){
    for (j=0; j<i;++j)
      mzd_write_bit(L, i, j, mzd_read_bit(A,i,j));
    for (j=i+1; j<n;++j)
      mzd_write_bit(U, i, j, mzd_read_bit(A,i,j));
  }
  for (i=r; i<m; i++)
    for (j=0; j<r;++j)
      mzd_write_bit(L, i, j, mzd_read_bit(A,i,j));
  for (i=0; i<r; ++i){
    mzd_write_bit(L,i,i, 1);
    mzd_write_bit(U,i,i, 1);
  }

  mzd_apply_p_left(Acopy, P);
  mzd_apply_p_right_trans(Acopy, Q);

  mzd_addmul(Acopy, L, U, 0);

  int status = 0;
  for ( i=0; i<m; ++i)
    for ( j=0; j<n; ++j){
      if (mzd_read_bit (Acopy,i,j)){
        status = 1;
        break;
      }
    }
  if (status) {
    printf(" ... FAILED\n");
  }  else
    printf (" ... passed\n");
  mzd_free(U);
  mzd_free(L);
  mzd_free(A);
  mzd_free(Acopy);
  mzp_free(P);
  mzp_free(Q);
  return status;
}


int main(int argc, char **argv) {
  int status = 0;

  status += test_pluq_structured(37, 37);
  status += test_pluq_structured(63, 63);
  status += test_pluq_structured(64, 64);
  status += test_pluq_structured(65, 65);
  status += test_pluq_structured(128, 128);

  status += test_pluq_structured(37, 137);
  status += test_pluq_structured(65, 5);
  status += test_pluq_structured(128, 18);

  status += test_pluq_full_rank(13,13);
  status += test_pluq_full_rank(37,37);
  status += test_pluq_full_rank(63,63);
  status += test_pluq_full_rank(64,64);
  status += test_pluq_full_rank(65,65);
  status += test_pluq_full_rank(97,97); 
  status += test_pluq_full_rank(128,128);
  status += test_pluq_full_rank(150,150);
  status += test_pluq_full_rank(256,256);
  status += test_pluq_full_rank(1024,1024);

  status += test_pluq_full_rank(13,11);
  status += test_pluq_full_rank(37,39);
  status += test_pluq_full_rank(64,164);
  status += test_pluq_full_rank(97,92);
  status += test_pluq_full_rank(128,121);
  status += test_pluq_full_rank(150,153);
  status += test_pluq_full_rank(256,258);
  status += test_pluq_full_rank(1024,1023);

  status += test_pluq_half_rank(64,64);
  status += test_pluq_half_rank(65,65);
  status += test_pluq_half_rank(66,66);
  status += test_pluq_half_rank(127,127);
  status += test_pluq_half_rank(129,129);
  status += test_pluq_half_rank(148,148);
  status += test_pluq_half_rank(132,132);
  status += test_pluq_half_rank(256,256);
  status += test_pluq_half_rank(1024,1024);

  status += test_pluq_half_rank(129,127);
  status += test_pluq_half_rank(132,136);
  status += test_pluq_half_rank(256,251);
  status += test_pluq_half_rank(1024,2100);

  status += test_pluq_random(63,63);
  status += test_pluq_random(64,64);
  status += test_pluq_random(65,65);

  status += test_pluq_random(128,128);
  status += test_pluq_random(128, 131);
  status += test_pluq_random(132, 731);
  status += test_pluq_random(150,150);
  status += test_pluq_random(252, 24);
  status += test_pluq_random(256,256);
  status += test_pluq_random(1024,1022);
  status += test_pluq_random(1024,1024);

  status += test_pluq_random(128,1280);
  status += test_pluq_random(128, 130);
  status += test_pluq_random(132, 132);
  status += test_pluq_random(150,151);
  status += test_pluq_random(252, 2);
  status += test_pluq_random(256,251);
  status += test_pluq_random(1024,1025);
  status += test_pluq_random(1024,1021);

  if (!status) {
    printf("All tests passed.\n");
    return 0;
  } else {
    return -1;
  }
}