/* gausssj_main_nofrills.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_nofrills.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"

#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;

    /* Do the work... */

    gaussj(a, N, b, 1);
}