/* matrices.h 
   written by Marco on august 1994 
   revised in july 2021
   fixed a bug in function allocm on 06 jan 2022
   fixed a bug in function freem on 06 jan 2022
   
   requires stdio.h and stdlib.h */


/* This header file is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This header file is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details. */




typedef struct matrix{
        double **el;
        int ordi;
        int ordj;}MAT;

const MAT MATNULL={NULL,0,0};  /* MATNULL = null matrix */

void submat(MAT *pmin, MAT mat, int h, int k);
void allocm(MAT *pmat, int ordi, int ordj);
void freem(MAT *pmat);
void writem(MAT mat);
void fwritem(FILE *fp, MAT mat);
void readm(MAT *pmat);
void freadm(FILE *fp, MAT *pmat);
void mcpy(MAT *pmat, MAT mat);
double sgn(int i, int j);
double det(MAT mat);
void transp(MAT *pdest, MAT sour);
void inv(MAT *pdest, MAT sour);
void mmult(MAT *dest, MAT mat1, MAT mat2);
void msum(MAT *dest, MAT mat1, MAT mat2);
void mmulsc(MAT *pdest, MAT source, double lambda);
void linsys(MAT *psol, MAT mat, MAT vcterms);



/* extract a submatrix from a matrix eliminating row h and col k */
void submat(MAT *pmin, MAT mat, int h, int k)
{
 int i,j,x,y;

 /* allocate the memory for the submatrix */
 allocm(pmin,mat.ordi-1,mat.ordj-1);

 /* fill the submatrix with the matrix elements, excluding row h and col k */
 for(i=1,x=1;i<=mat.ordi;i++,x++)
  {
   if(i==h)
      x--;
   else
    {
      for(j=1,y=1;j<=mat.ordj;j++,y++)
       {
        if(j==k)
          y--;
        else
         pmin->el[x][y]=mat.el[i][j];
       }
    }
  }
}




/* allocate memory for a matrix */
void allocm(MAT *pmat, int ordi, int ordj)
{
 int i,j;

 /* rows and columns must be >=1 */
 if(ordi<=0||ordj<=0)
           {fprintf(stderr,"Illegal order.\n");
            exit(EXIT_FAILURE);}

 /* allocate an array of pointers to the rows */
 if(!(pmat->el=(double **)(malloc((ordi+1)*sizeof(double *)))))
                           {fprintf(stderr,"Can't allocate memory.\n");
                            exit(EXIT_FAILURE);}
 /* allocate the rows */
 for(i=1;i<=ordi;i++)
   if(!(pmat->el[i]=(double *)(malloc((ordj+1)*sizeof(double)))))
                           {fprintf(stderr,"Can't allocate memory.\n");
			    exit(EXIT_FAILURE);}

 /* set the number of rows and columns into the matrix struct */
 pmat->ordi=ordi;
 pmat->ordj=ordj;
}




/* free the memory associated to a matrix */
void freem(MAT *pmat)
{
 int i;

 /* free the memory of the rows */
 for(i=1;i<=pmat->ordi;i++)
              free(pmat->el[i]);

 /* free the memory of the pointers to the rows */
 free(pmat->el);

 /* set to NULL the pointer of the matrix struct and to 0,0 rows and cols */
 *pmat=MATNULL;
}




/* write a matrix to stdout */
void writem(MAT mat)
{
 fwritem(stdout,mat);
}




/* write a matrix to file */
void fwritem(FILE *fp, MAT mat)
{
 int i,j;

 /* write the number of rows and cols */
 fprintf(fp,"%d\t%d\n",mat.ordi,mat.ordj);

 /* write the elements */
 for(i=1;i<=mat.ordi;i++)
  {for(j=1;j<=mat.ordj;j++)
        {fprintf(fp,"%9.6f",mat.el[i][j]);
         if(j<mat.ordj)
              fprintf(fp,"\t");
        }
   fprintf(fp,"\n");}
}



/* read a matrix from stdin */
void readm(MAT *pmat)
{
 int i,j;
 double c;

 /* read the number of rows and cols */
 printf("Rows? ");
 scanf("%d",&i);
 printf("Cols? ");
 scanf("%d",&j);

 /* allocate the memory for the matrix */
 allocm(pmat,i,j);

 /* read the elements */
 for(i=1;i<=pmat->ordi;i++)
  {for(j=1;j<=pmat->ordj;j++)
       {printf("a[%d][%d] = ",i,j);
        scanf("%lf",&c);
        pmat->el[i][j]=c;}
   printf("\n");}
}



/* read a matrix from file */
void freadm(FILE *fp, MAT *pmat)
{
 int i,j;
 double c;

 /* read the number of rows and cols */
 fscanf(fp,"%d",&i);
 fscanf(fp,"%d",&j);

 /* allocate the memory for the matrix */
 allocm(pmat,i,j);

 /* read the elements */
 for(i=1;i<=pmat->ordi;i++)
   for(j=1;j<=pmat->ordj;j++)
              {fscanf(fp,"%lf",&c);
               pmat->el[i][j]=c;}
}




/* duplicate a matrix */
void mcpy(MAT *pmat, MAT mat)
{
 int i,j;

 /* allocate the memory for the matrix */
 allocm(pmat,mat.ordi,mat.ordj);

 /* copy the source to the destination element by element */
 for(i=1;i<=mat.ordi;i++)
    for(j=1;j<=mat.ordj;j++)
       pmat->el[i][j]=mat.el[i][j];
}




/* return the sign of the cofactor */
double sgn(int i, int j)
{
 double s;
 int k,r;

 /* check if (i+j) is odd */
 k=i+j;
 r=k & 1;

 if(r==0) s=1.0;   /* if even, the sign of the cofactor is "+" */
 else s=-1.0;      /* if  odd, the sign of the cofactor is "-" */

 return(s);
}




/* Calculates the determinant of a square matrix */
double det(MAT mat)
{
 MAT temp=MATNULL;
 int j;
 double d=0.0;

 /* exit if the matrix is not square */
 if(mat.ordi!=mat.ordj)
    {fprintf(stderr,"Illegal operation on rectangular matrix.\n");
     exit(EXIT_FAILURE);}

 /* calculates the determinant recursively with Laplace expansion */
 if(mat.ordi==1)
   d=mat.el[1][1];
 else
   for(j=1;j<=mat.ordi;j++)
      {submat(&temp,mat,1,j);
       d+=sgn(1,j)*mat.el[1][j]*det(temp);
       freem(&temp);}

 /* return the determinant previously calculated */
 return(d);
}



/* transpose a matrix */
void transp(MAT *pdest, MAT sour)
{
 int i,j;
 
 /* allocate the memory for the transpose */
 allocm(pdest,sour.ordj,sour.ordi);

 /* fill the transpose */
 for(i=1;i<=sour.ordi;i++)
   for(j=1;j<=sour.ordj;j++)
     pdest->el[j][i]=sour.el[i][j];
}



/* invert a square non singular matrix */
void inv(MAT *pdest, MAT sour)
{
 MAT temp=MATNULL;
 int i,j;
 double delta;

 /* calculate the determinant ot the matrix */
 delta=det(sour);
 
 /* if the matrix is singular print an error message and exit */
 if(delta==0.0)
      {fprintf(stderr,"Illegal operation on singular matrix.\n");
       exit(EXIT_FAILURE);}
 
 /* allocate the memory for the inverse */
 allocm(pdest,sour.ordi,sour.ordj);
 
 /* calculate the elements ot the inverse with the cofactors */
 for(j=1;j<=sour.ordj;j++)
   for(i=1;i<=sour.ordi;i++)
       {
        submat(&temp,sour,j,i);
        pdest->el[i][j]=(sgn(i,j)*det(temp))/delta;
        freem(&temp);
       }
}




/* multiply two matrices mat1*mat2, rows by cols */
void mmult(MAT *dest, MAT mat1, MAT mat2)
{
 int n,m,p;

 /* if rows(mat1) is != cols(mat2) you can't multiply mat1*mat2 */
 if(mat2.ordi!=mat1.ordj)
         {fprintf(stderr,"Orders do not match.\n");
          exit(EXIT_FAILURE);}
 
 /* allocate the memory for the product */
 allocm(dest,mat1.ordi,mat2.ordj);
 
 for(n=1;n<=mat1.ordi;n++)
    for(m=1;m<=mat2.ordj;m++)
      {
       dest->el[n][m]=0.0;
       for(p=1;p<=mat1.ordj;p++)
              dest->el[n][m]+=mat1.el[n][p]*mat2.el[p][m];
      }
}


/* sum two matrices */
void msum(MAT *dest, MAT mat1, MAT mat2)
{
 int i,j;

 /* if the two matrices do not have the same number of rows and cols
    print an error message and exit */ 
 if(mat1.ordi!=mat2.ordi || mat1.ordj!=mat2.ordj)
                 {fprintf(stderr,"Orders do not match.\n");
                  exit(EXIT_FAILURE);}

 /* allocate the memory for the sum */
 allocm(dest,mat1.ordi,mat1.ordj);

 /* sum the two matrices element by element */
 for(i=1;i<=mat1.ordi;i++)
     for(j=1;j<=mat1.ordj;j++)
        dest->el[i][j]=mat1.el[i][j]+mat2.el[i][j];
}



/* multiply a matrix by the scalar lambda */
void mmulsc(MAT *pdest, MAT sour, double lambda)
{
 int i,j;
 
 /* allocate the memory for the product */
 allocm(pdest,sour.ordi,sour.ordj);

 /* multiply each element of the matrix by the scalar lambda */
 for(i=1;i<=sour.ordi;i++)
    for(j=1;j<=sour.ordj;j++)
       pdest->el[i][j]=lambda*sour.el[i][j];
}



/* solve a linear system with Cramer's rule */
void linsys(MAT *psol, MAT mat, MAT vcterms)
{
 MAT temp;

 /* calculate the inverse matrix of the coefficients */
 inv(&temp,mat);

 /* multiply the inverse matrix by the vector of the constant terms;
    this give the vector of the solutions */
 mmult(psol,temp,vcterms);

 freem(&temp);
}