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

#define Z 8

float** request(rows,cols)
{
        float **matrix;
        int i;
        
        matrix=malloc(rows*sizeof(float*));
        for (i=0;i<rows;i++)
                matrix[i]=malloc(cols*sizeof(float));
        
        return matrix;
}

void libera(float **pointer, int sizei)
{
        int i;

        for (i=0;i<sizei;i++)
                free(pointer[i]);
        free(pointer);
}

void pinta(float **A, int rows, int cols, int decimal)
{
        int i,j;

        for (i=0;i<rows;i++)
        {
                for (j=0;j<cols;j++)
                {
                        if (decimal)
                                printf("%3.2f ",A[i][j]);
                        else
                                printf("%3.0f ",A[i][j]);
                }
                printf("\n");
        }
        printf("---------------------------------------------------\n");        
}

/* See pg 106 Kernighan & Ritchie (Prentice Hall) */
void swap(float *px, float *py)
{
        float tmp;        
        tmp=*px;        
        *px=*py;        
        *py=tmp;        
}

float** eye(int order)
{
        int l,m;
        float **I;

        I=request(order,order);

        for (l=0;l<order;l++)
                for (m=0;m<order;m++)
                        I[l][m]=0;

        for (l=0;l<order;l++)
                I[l][l]=1;

        return I;        
}

/* Gauss-Jordan */
/* http://www.cs.caltech.edu/courses/cs171/c2-1.pdf */
/* Original method solves Ax=B and B is not necessarily the Identity as here is        */
/* float** inverse(float** A, int n, float** B, int m)                                 */
float** inverse(float** input, int n)
{
        float **b,**A;
        int i,icol,irow,j,k,l,ll,*indxc,*indxr,*ipiv,m=n;
        float big,dum,pivinv;

        A=request(n,n);
        for (i=0;i<n;i++)
                for (j=0;j<n;j++)
                        A[i][j]=input[i][j];

        b=eye(m);
        indxc=malloc(n*sizeof(int));
        indxr=malloc(n*sizeof(int));
        ipiv=malloc(n*sizeof(int));

        for (j=0;j<n;j++)
                ipiv[j]=0;

        for (i=0;i<n;i++)
        {
                big=0;
                for (j=0;j<n;j++)
                {
                        if (ipiv[j]!=1)
                        {
                                for (k=0;k<n;k++)
                                {
                                        if (ipiv[k]==0)
                                        {
                                                if (fabs(A[j][k])>=big)
                                                {
                                                        big=fabs(A[j][k]);        
                                                        irow=j;
                                                        icol=k;
                                                }
                                        }
                                }
                        }
                }
                ++(ipiv[icol]);
                if (irow!=icol)
                {
                        for (l=0;l<n;l++)
                                swap(&A[irow][l],&A[icol][l]);
                        for (l=0;l<m;l++)
                                swap(&b[irow][1],&b[icol][l]);
                }
                indxr[i]=irow;
                indxc[i]=icol;
                if (A[icol][icol]==0)
                        printf("Error: singular matrix (result will be wrong) \n");
                else
                        pivinv=1/A[icol][icol];
                A[icol][icol]=1;                
                for (l=0;l<n;l++)
                {
                        A[icol][l]=A[icol][l]*pivinv;
                        b[icol][l]=b[icol][l]*pivinv;
                }                
                for (ll=0;ll<n;ll++)
                {
                        if (ll!=icol)
                        {
                                dum=A[ll][icol];
                                A[ll][icol]=0;
                                for (l=0;l<n;l++)
                                {
                                        A[ll][l]=A[ll][l]-(A[icol][l]*dum);
                                        b[ll][l]=b[ll][l]-(b[icol][l]*dum);
                                }
                        }
                }        
        }
        for (l=n-1;l>=0;l--)
        {
                if (indxr[l]!=indxc[l])
                {
                        for (k=0;k<n;k++)
                                swap(&A[k][indxr[l]],&A[k][indxc[l]]);
                }
        }

        free(indxc);
        free(indxr);
        free(ipiv);
        libera(b,n);

        return A;
}

float** scalar_division(float **matrix, int rows, int cols, float scalar)
{
        float **result;
        int i,j;

        result=request(rows,cols);

        for (i=0;i<rows;i++)
                for (j=0;j<cols;j++)
                        result[i][j]=matrix[i][j]/scalar;

        return result;
}

float** scalar_product(float **matrix, int rows, int cols, float scalar)
{
        float **result;
        int i,j;

        result=request(rows,cols);

        for (i=0;i<rows;i++)
                for (j=0;j<cols;j++)
                        result[i][j]=matrix[i][j]*scalar;

        return result;
}

float** dot_product(float **A, float **B,int rows, int cols)
{
        int i,j;
        float **dot_product;

        dot_product=request(rows,cols);        
        for (i=0;i<rows;i++)
                for (j=0;j<cols;j++)
                        dot_product[i][j]=A[i][j]*B[i][j];
        
        return dot_product;
}

float** matrix_product(float **A, float **B,int Arows, int Acols, int Brows, int Bcols)
{
        float **product,sumatorio;
        int i,j,l;

        if ((Arows==Bcols)&&(Acols==Brows))
        {
                product=request(Arows,Arows);                
                for (i=0;i<Arows;i++)
                        for (j=0;j<Arows;j++)
                        {
                                sumatorio=0;
                                for (l=0;l<Acols;l++)
                                        sumatorio=sumatorio+(A[i][l]*B[l][j]);
                                product[i][j]=sumatorio;
                        }                        
                return product;
        }        
        else
                return NULL;

}

int main()
{
	float source[Z][Z]={
	      {157,157,159,154,158,154,155,158},{158,158,156,155,157,158,155,159},
	      {160,157,156,153,156,171,156,155},{154,158,156,155,154,160,159,155},
	      {151,160,157,156,156,155,155,159},{158,159,155,156,154,171,157,155},
	      {153,155,158,157,155,157,153,156},{159,158,154,157,155,156,155,153}};
        float **A,**B,**C,**D,**I,**iA;
        int i,j,z=Z;

        A=request(z,z);
        for (i=0;i<z;i++)
          for (j=0;j<z;j++)
            A[i][j]=source[i][j];
        printf("A:\n");
        pinta(A,z,z,0);

        B=scalar_division(A,z,z,10);
        printf("B=A/10\n");
        pinta(B,z,z,1);

        C=dot_product(B,B,z,z);
        printf("C = B .x B :\n");
        pinta(C,z,z,1);

        D=matrix_product(A,B,z,z,z,z);
        printf("D = A x B :\n");
        pinta(D,z,z,1);

        printf("Inv(A): \n");        

        iA=inverse(A,z);        
        pinta(iA,z,z,1);

        printf("I = Inv(A) x A \n");        
        I=matrix_product(iA,A,z,z,z,z);
        pinta(I,z,z,0);

        libera(A,z);
        libera(B,z);
        libera(C,z);
        libera(D,z);
        libera(I,z);
}