#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <stdbool.h>
#include <ulmaux.h>
#include <ulmblas.h>
#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();
}
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);
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);
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);
}