#include "sym_mat.h"
#include <stdio.h>

// solve Ax=b, A=[2 -1 0; -1 2 -1; 0 -1 2], b = [1; 3; 2]
// compute c=A*b

int main() {
    sym_mat A;
    int i;
    double berr, ferr;
    double x[3], b[3] = {1.0, 3.0, 2.0}, c[3];

    sym_mat_init(&A, 3);

    sym_mat_add_entry(&A, 0, 0, 2.0); // ignores upper triangle
    sym_mat_add_entry(&A, 1, 1, 2.0);
    sym_mat_add_entry(&A, 2, 2, 2.0);
    sym_mat_add_entry(&A, 1, 0, -1.0);
    sym_mat_add_entry(&A, 2, 1, -1.0);

    sym_mat_factor(&A, "Equil", "ReOrder");

    printf("flops to factor: %g\n", A.flops);
    printf("           nnzA: %d\n", A.nnzA);
    printf("           nnzL: %d\n", A.nnzL);
    printf("reciprocal of condition number: %g\n", sym_mat_rcond(&A));

    sym_mat_solve(&A, b, x);

    // optional iterative refinement:
    sym_mat_refine(&A, b, x, &ferr, &berr);
    printf("berr = %g, ferr = %g\n", berr, ferr);

    sym_mat_apply(&A, b, c);

    for (i=0; i<3; i++)
	printf("b[%d] = %g\n", i, b[i]);
    for (i=0; i<3; i++)
	printf("x[%d] = %g\n", i, x[i]);
    for (i=0; i<3; i++)
	printf("c[%d] = %g\n", i, c[i]);

    sym_mat_free(&A);

    return 0;
}