/****************************************************************************
// Test which uses routine LU to solve a linear system using LU factorization.
//
// Usage: a.out < matrix.dat
//
// where matrix.dat is an ASCII file consisting of the
// matrix size (M,N) followed by its values.  For example,
// enter a 3x3 matrix as: 
//  8.1  1.2  4.3
//  1.3  4.3  2.9
//  0.4  1.3  6.1
//
     C.A. Bertulani 04/07/2000
****************************************************************************/

typedef double Number; 

#include <stdio.h>
#include<math.h>
#include<iostream.h>
#include<iomanip.h>
#include "butil.h"



int main()
{

void LU(Number **,int,int*,Number*);
void LUSOLVE(Number **,int,int*,Number b[]);

    Number **A, **O, **y, *b, *x, *col, d, ax;
	int N, *indx, i, j, k;

	cout << " Enter number of rows" << endl;
	cin >> N;

	A = matrix(N,N);
	O = matrix(N,N);
	y = matrix(N,N);
	x = vector(N);
	b = vector(N);
	col = vector(N);
	indx = ivector(N);

	for(i=1;i<=N;i++)b[i]=1;
    
	cout << " Enter Matrix A[i][j]" << endl;

	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++) cin >> A[i][j];
	}
	

    cout << "Original Matrix A: " <<  endl;
	
	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++){
			O[i][j] = A[i][j];			/* saving original matrix A */
			cout <<"  " <<  O[i][j];
		}
		cout  << endl;
	}

	LU(O,N,indx,&d);
	LUSOLVE(O,N,indx,b);
    

    cout << "Solution x for Ax=b, where b=[1,1,...] " <<endl;

	for(i=1;i<=N;i++)cout << " x["<<i<<"]: " << b[i] << endl;

    cout << "If your result is right"  <<endl;

    cout << "A*x should be the vector [1,1,...]. Let us check: "  <<endl;

//  original matrix A has been LU decomposed into O. 
	
	cout << "A.x = [ ";
	for(i=1;i<=N;i++){
		ax = 0.;
		for(j=1;j<=N;j++)ax += A[i][j]*b[j];
			cout << ax << " ";
	}
		cout  << "]" << endl;

// Now, to solve A.x=b with a different b, just enter new b and repeat the call
// LUSOLVE(O,N,indx,b), not with A, but with O and indx, as already set by LU;

//  Inverse of matrix A  

	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++){
			O[i][j] = A[i][j];			/* saving original matrix A */
		}
	}

	LU(O,N,indx,&d);		/* Decompose the matrix just once */
	for(j=1;j<=N;j++){		/* find inverse by columns */
		for(i=1;i<=N;i++)col[i]=0.0;
		col[j] = 1.0;
    	LUSOLVE(O,N,indx,col);
		for(i=1;i<=N;i++)y[i][j]=col[i];
	}

	cout << "Inverse of Matrix A: " <<  endl;
	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++) cout <<"  " <<  y[i][j];
		cout  << endl;
	}

    cout << "Is this right? Let us check. A.inv(A) ="  <<endl;

	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++){
			ax = 0.;
			for(k=1;k<=N;k++)ax += y[i][k]*A[k][j];
			cout << ax << " ";
		}
		cout  << endl;
	}

//  Determinant of matrix A  

/* In this case, call LU just ONE TIME: */

	for(i=1;i<=N;i++){
		for(j=1;j<=N;j++){
			O[i][j] = A[i][j];			/* saving original matrix A */
		}
	}

	LU(O,N,indx,&d);  /* this retuns d as +-1 */
	for(j=1;j<=N;j++) d *= O[j][j];
	cout << "... and you determinant is:  " << d;
	cout  << endl;


	return 0;
}

/******************************************************************** 
Solves the set of n linear equations A.X = B, where a[1..n][1..n] is
input, not as the matrix A but rather as its LU decomposition, determined
by the routine LU. indx[1..n] is input as the permutation vector 
returned by LU. b[1..n] is input as the right-hand side vector B, and
returns with the solution vector X. a,n, and indx are not modified by this
routine and can be left in place for successive calls with different 
rigt-hand sides b. This routine takes into account th possibility that
b will begin with many zero elements, so it is efficient for use in
matrix inversion. 
*********************************************************************/

void LUSOLVE(Number **a,int n,int *indx, Number b[])

{
	int i,ii=0,ip,j;
	Number sum;

	for (i=1;i<=n;i++) {  /* When ii is set to a positive value, it will */
		ip=indx[i];			/* become the index of the first novanishing element */
		sum=b[ip];			/* of b. We now do the forward substitution. The only */
		b[ip]=b[i];			/* new wrinkle is to unscramble the permutation as we go. */
		if (ii)
			for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
		else if (sum) ii=i;		/* A nonzero element was encountered, so from now on we */
		b[i]=sum;				/* will have to do the sums in the loop above */
	}
	for (i=n;i>=1;i--) {		/* Now we do the backsubstitution */
		sum=b[i];
		for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];       /* Store a component of the solution vector X */
	}
}

/******************************************************************** 
Given a matrix  a[1..n][1..n] this routine replaces it by the
LU decomposition of a rowwise permutation of itself. a and n are
input. a is also output, indx[1..n] is an output vector that records
the row permutation effected by the partial pivoting, d is output
as +-1 depending on whether the number of ro interchanges was even or
odd, respectively. 
*********************************************************************/

#define TINY 1.0e-20;

void LU(Number **a,int n,int *indx, Number *d)

{
	int i,imax,j,k;
	Number big,dum,sum,temp;
	Number *vv;

	vv=vector(n);			/* vv stores the implicit scaling of each row */
	*d=1.0;					/* No row interchanges yet */
	for (i=1;i<=n;i++) {		/* Loop over rows to get the implicit scaling */								
		big=0.0;				/* information */
		for (j=1;j<=n;j++)
			if ((temp=fabs(a[i][j])) > big) big=temp;
			if (big == 0.0) cout << "Singular matrix in routine LU" << endl;
			/* No nonzero largest element */
        		vv[i]=1.0/big;			/* Save the scaling */
	}        
			for (j=1;j<=n;j++) {		/* This is the loop over columns of */
				for (i=1;i<j;i++) {		/* Crout's method */
					sum=a[i][j];
					for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
					a[i][j]=sum;
				}
		big=0.0;			/* Initialize for the search for largest pivot element */
		for (i=j;i<=n;i++) {
			sum=a[i][j];
			for (k=1;k<j;k++)
				sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big) {
				/* Is the figure of merit for the pivot better than the best so far? */
				big=dum;
				imax=i;
			}
		}
		if (j != imax) {   /* Do we need to interchange rows? */
			for (k=1;k<=n;k++) {   /* Yes, do so... */
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			*d = -(*d);			/* ...and change the parity of d */
			vv[imax]=vv[j];		/* Also interchange the scale factor */
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != n) {			/* Now, finally, divide by the pivot element */
			dum=1.0/(a[j][j]);
			for (i=j+1;i<=n;i++) a[i][j] *= dum;
		}
	}				/* Go back for the next column in the reduction */
	delete [] vv;
}

#undef TINY