/* Homework 2, problem 2.  Compile with
 *
 *	make hw2_2
 *
 * Run with
 *
 *	hw2_2  < hw2.2.dat
 */

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

#include <stdio.h>
#include <math.h>
#include "nrutil.h"
#include "nr.h"

int main()
{
    int i, j, k, n;
    int *index;
    double d;
    double **a, **asave, **ainv, **prod, *col;
    double max_diag, max_offd;

    /* Read in the matrix dimension and initialize. */

    scanf("%d", &n);
    printf("n = %d\n", n);

    a = dmatrix(1, n, 1, n);
    asave = dmatrix(1, n, 1, n);
    index = ivector(1, n);

    /* Read in the matrix. */

    for (i = 1; i <= n; i++)
	for (j = 1; j <= n; j++) {
	    scanf("%lf", &a[i][j]);
	    asave[i][j] = a[i][j];	/* (a will be overwritten below) */
	}

    /* LU decompose a.  On return, L and U are stored in a. */

    ludcmp(a, n, index, &d);

    /* Compute the determinant (product of d and the diagonal elements of U). */

    for (i = 1; i <= n; i++) {
	printf("a[%d][%d] = %f\n", i, i, a[i][i]);
	d *= a[i][i];
    }
    printf("determinant = %g\n", d);

    /* Invert the original matrix column by column (NR, p. 48). */

    ainv = dmatrix(1, n, 1, n);
    col = dvector(1, n);

    for (j = 1; j <= n; j++) {
	for (i = 1; i <= n; i++) col[i] = 0;
	col[j] = 1;
	lubksb(a, n, index, col);
	for (i = 1; i <= n; i++) ainv[i][j] = col[i];
    }

    /* Print out some elements of the inverse matrix. */

    printf("\ninverse first row:\n");
    k = n;
    if (n > 20) k = 20;
    for (j = 1; j <= k; j++)
	printf("%f ", ainv[1][j]);
    printf("\n");

    printf("\ninverse last row:\n");
    k = n - 19;
    if (k < 1) k = 1;
    for (j = k; j <= n; j++)
	printf("%f ", ainv[n][j]);
    printf("\n");

    /* Check that a*ainv = 1. */

    prod = dmatrix(1, n, 1, n);

    max_diag = 0;
    max_offd = 0;

    for (i = 1; i <= n; i++) {
	for (j = 1; j <= n; j++) {
	    prod[i][j] = 0;
	    for (k = 1; k <= n; k++)
		prod[i][j] += asave[i][k]*ainv[k][j];
	    if (i == j) {
		if (fabs(prod[i][j] - 1) > max_diag)
		    max_diag = fabs(prod[i][j] - 1);
	    } else {
		if (fabs(prod[i][j]) > max_offd)
		    max_offd = fabs(prod[i][j]);
	    }
	}
    }

    printf("\nmax diagonal error = %e\nmax off-diagonal error = %e\n",
	   max_diag, max_offd);
}