/* gausssj_main.c:
 *
 * Template for Gauss-Jordan linear equation solution using Numerical
 * Recipes routines. The files nrutil.h, nrutil.c, and gaussj.c (all
 * obtained from the NR source directory, /usr/local/Numerical_Recipes/C
 * on newton) should be present in your working directory.  Compile this
 * program with
 *
 *	gcc gaussj_main.c nrutil.c -o gaussj -lm
 *
 * Note the ideosyncratic definitions of vectors and matrices.  They
 * are at least consistent throughout the entire NR package...
 */

/*------------------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>

#include "nrutil.h"
#include "nr.h"

//extern "C" double** dmatrix(long int, long int, long int, long int);
//extern "C" void gaussj(double**, int, double**, int);

#define N 5	/* dimension of the system */

int main()
{
    /* Solve the linear system Ax = b by Gauss-Jordan elimination. */

    int i, j;

    /* NR declarations of matrices: */

    double **a  = dmatrix(1, N, 1, N);
    double **b  = dmatrix(1, N, 1, 1);	/* b is really a column vector,
					   but gaussj expects a matrix
					   since we can solve multiple
					   equations with a single call */

    /* Set up A and b.  Note the NR convention that indices start at 1. */

    for (i = 1; i <= N; i++)		/* i = row, j = column */
	for (j = 1; j <= N; j++)
	    a[i][j] = sqrt(i+j);	/* or whatever... */

    for (i = 1; i <= N; i++)
	b[i][1] = i;

    /*--------------------------------------------------------------------*/

    /* Print out the input vector. */

    printf("b:\n");
    for (i = 1; i <= N; i++)
	printf("    %10.2f\n", b[i][1]);

    /* Print out the inverse matrix. */

    printf("A:\n");
    for (i = 1; i <= N; i++) {
	for (j = 1; j <= N; j++)
	    printf("%10.2f ", a[i][j]);
	printf("\n");
    }

    /*--------------------------------------------------------------------*/

    /* Do the work!  Solve by Gauss-Jordan.  On return, A is replaced by
       its inverse and b is replaced by the solution matrix. */

    gaussj(a, N, b, 1);

    /*--------------------------------------------------------------------*/

    /* Print out the solution vector. */

    printf("x:\n");
    for (i = 1; i <= N; i++)
	printf("    %10.2f\n", b[i][1]);

    /* Print out the inverse matrix. */

    printf("Ainv:\n");
    for (i = 1; i <= N; i++) {
	for (j = 1; j <= N; j++)
	    printf("%10.2f ", a[i][j]);
	printf("\n");

    /*--------------------------------------------------------------------*/

    }
}