/* This is a sample program to show how to find the determinant and 
   the inverse of a given matrix. Here ludcmp() and lubksb() use index from 0 to n-1.
*/

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

#define TINY 1.0e-20;


main()
{
  float **a,**b,**c,**d,det,x,y;
  float **matrix(),*vector();
  int i,j,k,l,m,n;
  float invdet();

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

  a=matrix(0,n,0,n);
  b=matrix(0,n,0,n);
  c=matrix(0,n,0,n);
  d=matrix(0,n,0,n);

  for (i=0; i<n; i++)          /* creating an arbitrary matrix */
    for (j=0; j<n; j++)
      {	x=i; y=j;
	d[i][j]=a[i][j]=2*i+5*j+10.0/(2.0+abs(i-j)); 
      }

  if (n==3) {
    det=a[0][0]*a[1][1]*a[2][2]+a[1][0]*a[2][1]*a[0][2]+a[2][0]*a[0][1]*a[1][2]
      -a[0][2]*a[1][1]*a[2][0]-a[1][0]*a[0][1]*a[2][2]-a[0][0]*a[1][2]*a[2][1];
    printf("det=%f\n",det);
  }
  prtmtx(a,n,n); 

  det=invdet(a,n,b);

  printf("The determinant: %f\n",det);
  printf("\nThe inverse matrix:\n");

  prtmtx(b,n,n); 

  for (i=0; i<n; i++)
    for (j=0; j<n; j++){
      c[i][j]=0.0;
      for (k=0; k<n; k++)
	c[i][j]+=d[i][k]*b[k][j];
    }

  printf("\nVerification (should be identity matrix)\n");
  prtmtx(c,n,n); 
  
}

float invdet(a,n,inv)
/* This function computes both the determinant of matrix a and its inverse matrix */

float **a,**inv;
int n;
{
  int i,j,k,l,m,*indx,*ivector();
  float d,**x,**y,*col,*vector(),**matrix();
  void ludcmp();
  void lubksb();

  col=vector(0,n-1);
  indx=ivector(0,n-1);

  ludcmp(a,n,indx,&d);
  for (j=0; j<n; j++) {
    d *= a[j][j];
    for (i=0; i<n; i++) col[i]=0.0;
    col[j]=1.0;
    lubksb(a,n,indx,col);
    for (i=0; i<n; i++) inv[i][j]=col[i];
  }
  return d;

}


prtvtr(v,n)
/* print an n by 1 vector */
float v[];
int n;
{ int i;
  for (i=0; i<n; i++)
    printf("%7.2f ",v[i]);
}

prtmtx(a,m,n)
/* print an n by n matrix */
float **a;
int m,n;
{ int i,j;
  for (i=0; i<m; i++){
    for (j=0; j<n; j++){
      printf("%7.2f ",a[i][j]);
    }
    printf("\n");
  }
}


void ludcmp(a,n,indx,d)
int n,*indx;
float **a, *d;
{
  int i,imax,j,k;
  float big,dum,sum,temp;
  float *vv, *vector();

  vv=vector(0,n-1);
  *d=1.0;
  for (i=0; i<n;i++) {
    big=0.0;
    for (j=0; j<n; j++) 
      if ((temp=fabs(a[i][j])) > big) big=temp;
    if (big == 0.0) printf("Singular matrix\n");
    vv[i]=1.0/big;
  }
  for (j=0; j<n; j++) {
    for (i=0; i<j; i++) {
      sum=a[i][j];
      for (k=0; k<i; k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
    }
    big=0.0;
    for (i=j; i<n; i++) {
      sum=a[i][j];
      for (k=0; k<j; k++)
	sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
      if ( (dum=vv[i]*fabs(sum)) >= big) {
	big=dum;
	imax=i;
      }
    }
    if ( j != imax) {
      for (k=0; k<n; k++) {
	dum=a[imax][k];
	a[imax][k]=a[j][k];
	a[j][k]=dum;
      }
      *d = -(*d);
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j][j] == 0.0) a[j][j]=TINY;
    if (j != n-1) {
      dum=1.0/a[j][j];
      for (i=j+1; i<n; i++) a[i][j] *= dum;
    }
  }
}

void lubksb(a,n,indx,b)
float **a, b[];
int n,*indx;
{
  int i,ii=0,ip,j;
  float sum;

  for (i=0; i<n; i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii>=0) 
      for (j=ii; j<=i-1; j++) sum -= a[i][j]*b[j];
    else if (sum) ii=i;
    b[i]=sum;
  }
  for (i=n-1; i>=0; i--) {
    sum=b[i];
    for (j=i+1; j<n; j++) sum -= a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}




int *ivector(nl,nh)
int nl,nh;
{
  int *v;
  v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
  if (!v) printf("allocation failure in ivector()");
  return v-nl;
}


float *vector(nl,nh)
int nl,nh;
{
  float *v;
  v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
  if (!v) printf("allocation failure in vector()");
  return v-nl;
}


void free_vector(v,nl,nh)
float *v;
int nl,nh;
{
  free((char*) (v+nl));
}

float **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{ 
  int i;
  float **m;
  m=(float **)malloc((unsigned) (nch-ncl+1)*sizeof(float*));
  if (!m) printf("allocation failure 1 in matrix()");
  m -= nrl;

  for (i=nrl; i<=nrh; i++) {
    m[i]=(float *)malloc((unsigned) (nch-ncl+1)*sizeof(float));
    if (!m[i]) printf("allocation failure 2 in matrix()");
    m[i] -= ncl;
  }
  return m;
}