Some old matrix profiling code

matperf.c

#include "flint/fmpz_mat.h"
#include "flint/profiler.h"

#define TIMEIT_PRINT1(__var, __timer, __reps) \
    __var = __timer->cpu*0.001/__reps;

#define TIMEIT_REPEAT1(__timer, __reps) \
    do \
    { \
        slong __timeit_k; \
        __reps = 1; \
        while (1) \
        { \
            timeit_start(__timer); \
            for (__timeit_k = 0; __timeit_k < __reps; __timeit_k++) \
            {

#define TIMEIT_END_REPEAT1(__timer, __reps) \
            } \
            timeit_stop(__timer); \
            if (__timer->cpu >= 30) \
                break; \
            __reps *= 10; \
        } \
    } while (0);

#define TIMEIT_START1 \
    do { \
        timeit_t __timer; slong __reps; \
        TIMEIT_REPEAT1(__timer, __reps)

#define TIMEIT_STOP1(__var) \
        TIMEIT_END_REPEAT1(__timer, __reps) \
        TIMEIT_PRINT1(__var, __timer, __reps) \
    } while (0);

#define TIMEIT_REP_START \
    slong __reps = 0; \
    double __best = 1e300, __time = 0; \
    do { \
        TIMEIT_START1

#define TIMEIT_REP_STOP(__result) \
        TIMEIT_STOP1(__time); \
        __best = FLINT_MIN(__best, __time); \
        __reps++; \
    } while (__reps < 10 && __time < 2.0); \
   __result = __best;

void
fmpz_mat_pascal(fmpz_mat_t mat, int triangular)
{
    slong R, C, i, j;

    R = fmpz_mat_nrows(mat);
    C = fmpz_mat_ncols(mat);

    if (R == 0 || C == 0)
        return;

    if (triangular == 1)
    {
        for (i = 1; i < R; i++)
            for (j = 0; j < i && j < C; j++)
                fmpz_zero(fmpz_mat_entry(mat, i, j));

        for (j = 0; j < C; j++)
            fmpz_one(fmpz_mat_entry(mat, 0, j));

        for (i = 1; i < R && i < C; i++)
            fmpz_one(fmpz_mat_entry(mat, i, i));

        for (i = 1; i < R; i++)
            for (j = i + 1; j < C; j++)
                fmpz_add(fmpz_mat_entry(mat, i, j),
                    fmpz_mat_entry(mat, i, j - 1),
                    fmpz_mat_entry(mat, i - 1, j - 1));
    }
    else if (triangular == -1)
    {
        for (i = 0; i < R; i++)
            for (j = i + 1; j < C; j++)
                fmpz_zero(fmpz_mat_entry(mat, i, j));

        for (i = 0; i < R; i++)
            fmpz_one(fmpz_mat_entry(mat, i, 0));

        for (i = 1; i < R && i < C; i++)
            fmpz_one(fmpz_mat_entry(mat, i, i));

        for (i = 2; i < R; i++)
            for (j = 1; j < i && j < C; j++)
                fmpz_add(fmpz_mat_entry(mat, i, j),
                    fmpz_mat_entry(mat, i - 1, j - 1),
                    fmpz_mat_entry(mat, i - 1, j));
    }
    else
    {
        for (j = 0; j < C; j++)
            fmpz_one(fmpz_mat_entry(mat, 0, j));

        for (i = 1; i < R; i++)
            fmpz_one(fmpz_mat_entry(mat, i, 0));

        for (i = 1; i < R; i++)
            for (j = 1; j < C; j++)
                fmpz_add(fmpz_mat_entry(mat, i, j),
                    fmpz_mat_entry(mat, i, j - 1),
                    fmpz_mat_entry(mat, i - 1, j));
    }
}

int main()
{
    fmpz_mat_t A, B, C;
    slong n, i, bits;
    double t1, t2, t3;

    flint_rand_t state;
    flint_randinit(state);

    /* Comparison 1: randtest matrices; show relative timings to look for outliers  */
    if (0)
    {
        for (i = 0; i < 100; i++)
        {
            slong m, n, k, bits2;

            do {
                m = n_randint(state, 300);
                n = n_randint(state, 300);
                k = n_randint(state, 300);
                bits = n_randint(state, 1000);
                bits2 = n_randint(state, 1000);
            } while ((double) m * n * k * (bits + bits2) > 400.0 * 400 * 400 * 400);

            fmpz_mat_init(A, m, n);
            fmpz_mat_init(B, n, k);
            fmpz_mat_init(C, m, k);

            fmpz_mat_randtest(A, state, bits);
            fmpz_mat_randtest(B, state, bits2);

            TIMEIT_START1
            fmpz_mat_mul(C, A, B);
            TIMEIT_STOP1(t1)

            TIMEIT_START1
            fmpz_mat_mul_multi_mod(C, A, B);
            TIMEIT_STOP1(t2)

            printf("%5ld %5ld %5ld %5ld %5ld    %4ld\n", m, n, k, fmpz_mat_max_bits(A), fmpz_mat_max_bits(B), FLINT_MIN(999, (slong) (100 * t2 / t1)));

            fmpz_mat_clear(A);
            fmpz_mat_clear(B);
            fmpz_mat_clear(C);
        }
    }

    /* Comparison 2: Pascal matrices */
    if (0)
    {
        for (n = 2; n <= 1000; n = FLINT_MAX(n+3,n*1.1))
        {
            fmpz_mat_init(A, n, n);
            fmpz_mat_init(B, n, n);
            fmpz_mat_init(C, n, n);

            fmpz_mat_pascal(A, 1);
            fmpz_mat_pascal(B, 1);

            TIMEIT_START1
            fmpz_mat_mul(C, A, B);
            TIMEIT_STOP1(t1)

            TIMEIT_START1
            fmpz_mat_mul_multi_mod(C, A, B);
            TIMEIT_STOP1(t2)

            printf("%5ld  %4ld\n", n, (slong) (100 * t2 / t1));

            fmpz_mat_clear(A);
            fmpz_mat_clear(B);
            fmpz_mat_clear(C);
        }
    }

    /* Comparison 3: big table for uniform matrices */
    if (0)
    {
        /* print the sizes */
        printf("       ");
        for (n = 2; n <= 2000; n = FLINT_MAX(n+3,n*1.3))
        {
            if (n < 1000)
                printf("%4ld", n);
            else
                printf("%5ld", n);
        }
        printf("\n\n");

        for (bits = 1; bits <= 25000; bits = FLINT_MAX(bits+3,bits*1.3))
        {
            printf("%5ld  ", bits);

            for (n = 2; n <= 2000; n = FLINT_MAX(n+3,n*1.3))
            {
                fmpz_mat_init(A, n, n);
                fmpz_mat_init(B, n, n);
                fmpz_mat_init(C, n, n);

                fmpz_mat_randbits(A, state, bits);
                fmpz_mat_randbits(B, state, bits);

                TIMEIT_START1
                fmpz_mat_mul(C, A, B);
                TIMEIT_STOP1(t1)

                TIMEIT_START1
                fmpz_mat_mul_multi_mod(C, A, B);
                TIMEIT_STOP1(t2)

                if (n < 1000)
                    printf("%4ld", (slong) (100.0 * t2 / t1));
                else
                    printf("%5ld", (slong) (100.0 * t2 / t1));

                fflush(stdout);

                if (t1 > 10.0)
                    break;

                fmpz_mat_clear(A);
                fmpz_mat_clear(B);
                fmpz_mat_clear(C);
            }

            printf("\n");
        }
    }
}