#include <stdio.h>          // for printf()
#include <stdlib.h>         // for malloc(), free(), rand(), srand()
#include <stddef.h>         // for size_t, ptrdiff_t
#include <stdbool.h>        // for typedef bool

#include <ulmaux.h>
#include <ulmblas.h>


//-- simple test program -------------------------------------------------------

#ifndef COLMAJOR
#define COLMAJOR 1
#endif

#ifndef SEED_RAND
#define SEED_RAND 0
#endif

int
main()
{
    printf("COLMAJOR = %d\n", COLMAJOR);

    srand(SEED_RAND);

    size_t      n       = 6;
    ptrdiff_t   incRowA = COLMAJOR ? 1 : n;
    ptrdiff_t   incColA = COLMAJOR ? n : 1;
    double      *A      = malloc(n*n*sizeof(double));

    ptrdiff_t   incX    = 1;
    double      *x      = malloc(n*sizeof(double));

    ptrdiff_t   incB    = 1;
    double      *b      = malloc(n*sizeof(double));

    ptrdiff_t   incP    = 1;
    size_t      *p      = malloc(n*sizeof(size_t));

    if (!A || !x || !b || !p) {
        abort();
    }

    // init A, x and b
    randDGeMatrix(n, n, false, A, incRowA, incColA);

    for (size_t i=0; i<n; ++i) {
        x[i*incX] = i+1;
    }

    dgemv(n, n, 1, A, incRowA, incColA, x, incX, 0, b, incB);

    printf("A=\n");
    printDGeMatrix(n, n, A, incRowA, incColA);

    printf("b =\n");
    printDGeMatrix(1, n, b, 0, incB);

    // compute LU factorization and overwrite b with solution
    dgetrf(n, n, A, incRowA, incColA, p, incP);
    dlaswp(n, 0, n-1, p, incP, A, incRowA, incColA);
    dtrsv(n, true, true, A, incRowA, incColA, b, incB);
    dtrsv(n, false, false, A, incRowA, incColA, b, incB);

    // compare computed solution (in b) with exact solution (in x)
    printf("computed solution: x =\n");
    printDGeMatrix(1, n, b, 0, incB);

    printf("exact solution: x =\n");
    printDGeMatrix(1, n, x, 0, incX);

    free(A);
    free(x);
    free(b);
    free(p);
}