#include "genocop.h"

/********************************************************************************/
/*                                                                              */
/*           FUNCTION NAME     :   inverse()                                    */
/*                                                                              */
/*           SYNOPSIS          :   void inverse(a,a_inverse,n)                  */
/*                                                                              */
/*           DESCRIPTION       :   This function finds the inverse matrix of a  */
/*                                  given float matrix                          */
/*                                                                              */
/*           FUNCTIONS CALLED  :   det(),                                       */
/*                                 matrix().                                    */
/*                                                                              */
/*           CALLING FUNCITONS :   main()                                       */
/*                                                                              */
/*           AUTHOR            :   Swarnalatha Swaminathan                      */
/*                                                                              */
/*           DATE              :   1/17/92                                      */
/*                                                                              */
/*                                                                              */
/*           REV            DATE            BY           DESCRIPTION            */
/*           ---            ----            --           -----------            */
/*                                                                              */
/*                                                                              */
/********************************************************************************/




void inverse(a,a_inverse,n)
MATRIX a,a_inverse;    /*the original and the inverse matrices*/
int n;                     /*the dimension of the matrix*/
{
 int i,j,k,k1,l,l1,n1;     /*counter variables*/
 float d1,                 /*determinat of the matrix*/
       v,w;
 MATRIX b;                /*a temporary float matrix*/

 b = matrix(1, n, 1, n);
 for(i=1; i<=n; i++)
   for(j=1; j<=n; j++)
    b[i][j] = 0.0;

 d1 = det(a, n);

 if(d1==0.0) printf("No inverse exists");
 else
 {
   v=-1.0;
   for(i=1;i<=n;i++)
    {
      v=-v;
      w=-1.0;
        for(j=1;j<=n;j++)
         {
           w=-w;
           n1=n-1;
             for(k=1;k<=n1;k++)
              {
                k1=k;
                if(k>=i) k1=k+1;
                  for(l=1;l<=n1;l++)
                   {
                     l1=l;
                     if(l>=j) l1=l+1;
                     b[k][l]=a[k1][l1];
                   }
              }
               a_inverse[j][i]=det(b,n1)/d1*v*w;
         }
    }
 }
 free_matrix(b,1,n,1);
}

/********************************************************************************/
/*                                                                              */
/*           FUNCTION NAME     :   det()                                        */
/*                                                                              */
/*           SYNOPSIS          :   float det(input_matrix,nl)                   */
/*                                                                              */
/*           DESCRIPTION       :   This function returns the determinant of a   */
/*                                  given float matrix                          */
/*                                                                              */
/*           FUNCTIONS CALLED  :   matrix().                                    */
/*                                                                              */
/*           CALLING FUNCITONS :   inverse(),                                   */
/*                                 p_equalities().                              */
/*                                                                              */
/*           AUTHOR            :   Swarnalatha Swaminathan                      */
/*                                                                              */
/*           DATE              :   1/17/92                                      */
/*                                                                              */
/*                                                                              */
/*           REV            DATE            BY           DESCRIPTION            */
/*           ---            ----            --           -----------            */
/*                                                                              */
/*                                                                              */
/********************************************************************************/




float det(input_matrix, nl)
MATRIX input_matrix;                /*input float matrix*/
int nl;                             /*order of the matrix*/
{
 int ia,ib,i1,j1,n2,nz,n3,l,l1,m,m1;
 float d,x;
 MATRIX temp_input_matrix,b1;

 temp_input_matrix = matrix(1,nl,1,nl);
 b1 = matrix(1,nl,1,nl);

 for(i1=1;i1<=nl;i1++)
  for(j1=1;j1<=nl;j1++)
   temp_input_matrix[i1][j1]=input_matrix[i1][j1];

 d=1.0;
 for(ia=1;ia<=nl;ia++)
  {
   x=-1.0;
   n2=nl-ia+1;
     n3=n2-1;
   nz=0;
   for(ib=1;ib<=n2;ib++)
    if(temp_input_matrix[ib][1]!=0.0) nz=ib;
   if(nz==0.0) d=0.0;
   else
    {
     for(i1=1;i1<=nz;i1++)
      x=-x;
     d*=temp_input_matrix[nz][1]*x;

     for(l=1;l<=n3;l++)
      {
       l1=l;
       if(l>=nz) l1=l+1;
       for(m=1;m<=n3;m++)
        {
         m1=m+1;
         b1[l][m] = temp_input_matrix[l1][m1] - temp_input_matrix[l1][1]
                    /temp_input_matrix[nz][1]*temp_input_matrix[nz][m1];
       }
     }
   }
     for(i1=1;i1<=n3;i1++)
       for(j1=1;j1<=n3;j1++)
         temp_input_matrix[i1][j1]=b1[i1][j1];

 }
 free_matrix(b1,1,nl,1);
 free_matrix(temp_input_matrix,1,nl,1);

 return(d);
 }