/************************************************************
* This program calculates the determinant of a real square  *
* matrix using a recursive function Deter based on Kramer's *
* formula.                                                  *
* --------------------------------------------------------- *
* Method:                                                   *
*                                                           *
* Let us take the example:                                  *
* n=5                                                       *
* Input matrix is:   ( 3 2 1 6 3 )                          *
*                    ( 6 4 4 0 3 )                          *
*                    ( 4 0 1 0 2 )                          *
*                    ( 3 4 2 4 2 )                          *
*                    ( 4 6 1 5 2 )                          *
*                                                           *
* According to Kramer's rule, the determinant det(A) is     *
* calculated by the following:                              *
*                                                           *
*              (4 4 0 3)         (6 4 0 3)       (6 4 0 3)  *
* det = 3 * det(0 1 0 2) -2 * det(4 1 0 2) + 1 * (4 0 0 2)  *
*              (4 2 4 2)         (3 2 4 2)       (3 4 4 2)  *
*              (6 1 5 2)         (4 1 5 2)       (4 6 5 2)  *
*                                                           *
*              (6 4 4 3)         (6 4 4 0)                  *
*      -6 * det(4 0 1 2) +3 * det(4 0 1 0)                  *
*              (3 4 2 2)         (3 4 2 4)                  *
*              (4 6 1 2)         (4 6 1 5)                  *
*                                                           *
* As can be seen, 3, -2, 1, -6 and 3 are the elements of    *
* the 1st row alternatively multiplied by 1 or -1 and the   *
* the submatrices are obtained by successively eliminating  *
* row=1 and column=1, then row=1 and col=2, etc. Then you   *
* can recursively calculate the detrminants of these sub-   *
* matrices by using the same Function Deter with size n-1   *
* instead of n, and so forth...                             *
* When n=2, the result is immediate by the formula:         *
* A(1,1)*A(2,2)-A(2,1)*A(1,2).                              *
*                                                           *
* Program organization:                                     *
* After input of given square matrix A, here of size n=5,   *
* we call the function Deter(n,A) which returns the value   *
* of determinant. This function calls itself recursively    *
* with n decreasing by one at each step, until n=2.         *
* Auxiliary function is pow(n) returning -(-1)^n and auxi-  *
* liary procedure is Submatrix(n,A,B) that extracts sub-    *
* matrix B from A by eliminating row 1 and column col. Note *
* that B is dynamically allocated and disposed of at each   *
* step to spare memory.                                     *
*                                                           *
* SAMPLE RUN:                                               *
*                                                           *
* Size of matrix: 5                                         *
* Input matrix  :                                           *
*   3.0000  2.0000  1.0000  6.0000  3.0000                  *
*   6.0000  4.0000  4.0000  0.0000  3.0000                  *
*   4.0000  0.0000  1.0000  0.0000  2.0000                  *
*   3.0000  4.0000  2.0000  4.0000  2.0000                  *
*   4.0000  6.0000  1.0000  5.0000  2.0000                  *
*                                                           *
* Determinant = -7.00000000000000E+0001                     *
*                                                           *
* --------------------------------------------------------- *
*                                                           *
*                C++ release by Jean-Pierre Moreau, Paris.  *
*                           (www.jpmoreau.fr)               *
*************************************************************
Note: do not use for n>10                                  */
#include <stdio.h>
#include <math.h>

#include <basis.h>
#include <vmblock.h>

#define  SIZE  11

  int  n;      // size of matrix
  REAL **A;    // pointer to matrix SIZExSIZE
  REAL det;    // value of determinant
  FILE *fp;    // input.file
  
  void *vmblock = NULL;
  int  i, j;
  REAL temp;

// return -(-1)^n
int pow(int n) {
  if (n % 2 == 0) return 1; else return -1;
}

//return submatrix without 1st line and colth column of A
void Submatrix(int n, int col, REAL **A, REAL **B) {
  int i,j;
  boolean flag;

  if (n<3) return;
  for (i=1; i<n; i++) {
    flag=FALSE;
    for (j=1; j<n; j++)
	  if (j!=col) {
        if (!flag)  B[i][j]=A[i+1][j];
        else B[i][j]=A[i+1][j+1];
      }
      else {
        flag=TRUE;
        B[i][j]=A[i+1][j+1];
      }
  }
}

//print square matrix A(n,n)
void Matprint(int n, REAL **A) {
  int i,j;
  for (i=1; i<=n; i++) {
    for (j=1; j<=n; j++)
      printf(" %3.0f", A[i][j]);
    printf("\n");
  }
}

//return determinant of matrix A
REAL Deter(int n, REAL **A) {
  int i;
  REAL d;
  REAL **B;
  void *vmblock1 = NULL;

  d=ZERO;
  if (n==2)
	  return (A[1][1]*A[2][2]-A[1][2]*A[2][1]);
  else if (n==3) {
    d = A[1][1]*(A[2][2]*A[3][3]-A[3][2]*A[2][3]);
    d = d-A[1][2]*(A[2][1]*A[3][3]-A[3][1]*A[2][3]);
    d = d+A[1][3]*(A[2][1]*A[3][2]-A[3][1]*A[2][2]);
    return d;
  }
  else {
	for (i=1; i<=n; i++) {
// allocate memory to store submatrix of size n-1 x n-1
      vmblock1 = vminit();  
      B = (REAL **) vmalloc(vmblock1, MATRIX,  n, n);  
      if (! vmcomplete(vmblock1))
        LogError ("No Memory", 0, __FILE__, __LINE__);
      if (n>2) {
        Submatrix(n,i,A,B);
        //recursive call of Deter()
        d += pow(i)*A[1][i]*Deter(n-1,B);
      }
      vmfree(vmblock1);
    }
    return d;
  }
}

void main()  {

// read size n of matrix A from input file
  fp = fopen("mat10.dat","r");
  fscanf(fp, "%d", &n);

// allocate memory for matrix A of size n x n
  vmblock = vminit();  
  A = (REAL **) vmalloc(vmblock, MATRIX,  n+1, n+1);  
  if (! vmcomplete(vmblock))
    LogError ("No Memory", 0, __FILE__, __LINE__);

// read matrix A from input file
  for (i=1; i<=n; i++)
	for (j=1; j<=n; j++) {
	  fscanf(fp, "%lf", &temp);
	  A[i][j] = temp;
	}

  fclose(fp);

  printf("\n Size of matrix: %d\n", n);
  printf(" Input matrix:\n");
  Matprint(n, A);
  printf(" Computing determinant...\n");

  det = Deter(n,A);

  printf("\n Determinant = %f\n\n", det);

  vmfree(vmblock);

}

// end of file deter2a.cpp