test.c

#define _XOPEN_SOURCE 600
#include <stdio.h>
#include <gmp.h>
#include "bw.h"
#include "ulong_extras.h"
#include <time.h>

static inline double myrand() {
    /* ******************************************************************
       Generates a random floating-point number between -1.0 and
       1.0. Not sure if this is the best precision possible, but is
       based on bitwise operations, and I believe it should be faster
       than using (float)rand()/RAND_MAX.

       We pick up the number generated by rand, and leave just the bits
       over the mantissa region. The signal and exponent part are
       replaced to generate a number between 1.0 and 1.9999999. We
       then subtract 1.0 to bring the range down to 0.0 and 0.9999999.

       Coded in 2010-09-24 by Nicolau Werneck <nwerneck@gmail.com>.
     *************************************************************** */
    static union {
        unsigned  int i;
        float f;
    } myrand;
    myrand.i = (rand() & 0x007fffff) | 0x3f800000;
    return myrand.f-1.0;
}

int
main(int argc, char * argv[])
{
    nmod_sparse_mat_t M;
    n_primes_t primes;
    slong N = atoi(argv[1]), i, j;
    mp_limb_t p;
    float prob;
    mzd_t *K;
    struct timespec start, finish;
    double elapsed;

    printf("using %ldx%ld matrix\n", N, N);

    /* create M randomly */
    n_primes_init(primes);
    nmod_sparse_mat_init(M, N, N, 2);
    /*
    for (j = 0; j < N; j++)
        nmod_sparse_mat_ensure_row_alloc(M, j, 30);
    */
    if (1)
    {
        for (p = n_primes_next(primes), i = 0; i < N; p = n_primes_next(primes), i++)
        {
            prob = 1.0f - pow(1.0f - 6.0f/(float)p, 8);
            for (j = 0; j < M->c; j++)
            {
                if (myrand() < prob)
                {
                    _nmod_sparse_mat_set_entry(M, i, j, M->row_supports[i], 1);
                }
            }
        }
    }
    else
    {
        for (j = 1; j < N; j++)
        {
            _nmod_sparse_mat_set_entry(M, j, 0, 0, 1);
            _nmod_sparse_mat_set_entry(M, j, j, 1, 1);
        }
    }


    /*nmod_sparse_mat_print_pretty(M);*/

    K = mzd_init(N, 1);

    clock_gettime(CLOCK_MONOTONIC, &start);

    bw(K, M);

    clock_gettime(CLOCK_MONOTONIC, &finish);

    elapsed = (finish.tv_sec - start.tv_sec);
    elapsed += (finish.tv_nsec - start.tv_nsec) / 1000000000.0;

    printf("%lf secs taken\n", elapsed);

    mzd_free(K);

    nmod_sparse_mat_clear(M);

    n_primes_clear(primes);

    return 0;
}