#include <stdio.h>

#define NMAX 100
#define tiny 1.0E-20


void ludcmp(double A[][NMAX+1], int N, int *INDX, int D);
void lubksb(double A[][NMAX+1], int N, int *INDX, double *B);

main()
{
      int i,j;
      int N,D,dummy=NMAX;
      double A[NMAX+1][NMAX+1];
      double B[NMAX+1];
      int INDX[NMAX+1];


      do
      {
         printf("Input N (<=%d) - dimension of the array : ",dummy);
         scanf("%d",&N);
      } while(N>NMAX);
	 
      for(i=1;i<=N;i++)
      {
         for(j=1;j<=N;j++)
         {
            printf("Element (%d, %d) of the LHS : ",i,j);
            scanf("%lf",&A[i][j]);
         }
      printf("Element %d of the RHS : ",i);
      scanf("%lf",&B[i]);
      }

      ludcmp(A,N,INDX,D);
      lubksb(A,N,INDX,B);

      for(i=1;i<=N;i++)
      {
         printf("Element %d of the result : %lf\n",i,B[i]);
      }

      return 0;

}


/*
* Given an NxN matrix A, with physical dimension NP, this routine replaces it
* by the LU decomposition of a rowwise permutation of itself.
* A and N are input. A is output; INDX is output vector which records the
* row permutation effected by the partial pivoting; D id output as +1 or -1
* depending on whether the number of row interchanges was odd or even,
* respectively. This routine is used in combination with lubksb to solve
* linear equations or invert a matrix.
*/

void ludcmp(double A[][NMAX+1],int N,int *INDX,int D)
{
      int i,j,k,imax;
      double VV[NMAX];
      double aamax,sum,dummy;

      D=1;
      for(i=1;i<=N;i++)
      {
         aamax=0;
         for(j=1;j<=N;j++)
         {
            if (abs(A[i][j]>aamax))
               aamax=abs(A[i][j]);
         }
         if (aamax==0.0)
         {
            printf("Singular Matrix. \n");
            exit(1);
         }
         VV[i]=1.0/aamax;
      }

      for(j=1;j<=N;j++)
      {
         for(i=1;i<=j-1;i++)
         {
            sum=A[i][j];
            for(k=1;k<=i-1;k++)
            {
               sum -= A[i][k]*A[k][j];
            }
            A[i][j]=sum;
         }
         aamax=0;
         for(i=j;i<=N;i++)
         {
            sum=A[i][j];
            for(k=1;k<=j-1;k++)
            {
               sum -= A[i][k]*A[k][j];
            }
            A[i][j]=sum;
            dummy=VV[i]*abs(sum);
            if (dummy>=aamax)
            {
               imax=i;
               aamax=dummy;
            }
         }
         if (j!=imax)
         {
            for(k=1;k<=N;k++)
            {
               dummy=A[imax][k];
               A[imax][k]=A[j][k];
               A[j][k]=dummy;
            }
            D=-D;
            VV[imax]=VV[j];
         }
         INDX[j]=imax;
         if (A[j][j]==0.0)
            A[j][j]=tiny;
         if (j!=N)
         {
            dummy=1.0/A[j][j];
            for(i=j+1;i<=N;i++)
            {
               A[i][j]=A[i][j]*dummy;
            }
         }
      }
}

/*
* Solves the set of N linear equations A.X=D. Here A is input, not as a matrix
* A but rather as its LU decomposition. INDX is the input as the permutation
* vector returned bu ludcmp. B is input as the right-hand side vector B, and
* returns with the solution vector X.
*/

void lubksb(double A[][NMAX+1],int N,int *INDX,double *B)
{
      int ii,i,j,k,ll;
      double sum;

      ii=0;
      for(i=1;i<=N;i++)
      {
         ll=INDX[i];
         sum=B[ll];
         B[ll]=B[i];
         if (ii!=0)
         {
            for(j=ii;j<=i-1;j++)
            {
               sum -= A[i][j]*B[j];
            }
	 }
         else if (sum!=0.0)
            ii=i;

         B[i]=sum;
      }

      for(i=N;i>=1;i--)
      {
         sum=B[i];
         if (i<N)
         {
            for(j=i+1;j<=N;j++)
            {
               sum -= A[i][j]*B[j];
            }
         }
         B[i]=sum/A[i][i];
      }
}