/****************************************************
* Calculate the determinant of a real square matrix *
* by the LU decomposition method                    *
* ------------------------------------------------- *
* SAMPLE RUN:                                       *
* (Calculate the determinant of real square matrix: *
*        | 1  0.5  0  0    0  0    0  0    0 |      *
*        |-1  0.5  1  0    0  0    0  0    0 |      *
*        | 0 -1    0  2.5  0  0    0  0    0 |      *
*        | 0  0   -1 -1.5  4  0    0  0    0 |      *
*        | 0  0    0 -1   -3  5.5  0  0    0 |      *
*        | 0  0    0  0   -1 -4.5  7  0    0 |      *
*        | 0  0    0  0    0 -1   -6  8.5  0 |      *
*        | 0  0    0  0    0  0   -1 -7.5 10 |      *
*        | 0  0    0  0    0  0    0 -1   -9 |  )   *
*                                                   *
*                                                   *
*  Determinant =   1.000000                         *
*                                                   *
*                                                   *
*                 C++ version by J-P Moreau, Paris. *
*                         (www.jpmoreau.fr)         *
*****************************************************
 Exact value is: 1
---------------------------------------------------*/
#include <basis.h>
#include <vmblock.h>

#include "lu.h"
  
int main()  {

  REAL **A;       // input matrix of size n x n
  int  *INDX;     // dummy integer vector 
  void *vmblock = NULL; 

// NOTA: index zero not used  here. 

  int i, id, j, n, rc;
  REAL det;

  n=9;

  vmblock = vminit();  
  A       = (REAL **)vmalloc(vmblock, MATRIX,  n+1, n+1);  
  INDX    = (int *)  vmalloc(vmblock, VEKTOR,  n+1, 0);

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;  //not enough memory
  }

  for (i=0; i<=n; i++)
    for (j=0; j<=n; j++)
      A[i][j] = 0.0;

// define matrix A column by column
  A[1][1]=1.0; A[2][1]=-1.0;
  A[1][2]=0.5; A[2][2]=0.5; A[3][2]=-1.0;
  A[2][3]=1.0; A[4][3]=-1.0;
  A[3][4]=2.5; A[4][4]=-1.5; A[5][4]=-1.0;
  A[4][5]=4.0; A[5][5]=-3.0; A[6][5]=-1.0;
  A[5][6]=5.5; A[6][6]=-4.5; A[7][6]=-1.0;
  A[6][7]=7.0; A[7][7]=-6.0; A[8][7]=-1.0;
  A[7][8]=8.5; A[8][8]=-7.5; A[9][8]=-1.0;
  A[8][9]=10.0; A[9][9]=-9.0;

  //call LU decomposition
  rc = LUDCMP(A,n,INDX,&id);

  //calculate determinant and display
  det = id;
  if (rc==0)  {
    for (i=1; i<=n; i++)  det *= A[i][i];
    printf("\n Determinant = %f\n\n", det);
	return 0;
  }
  else {
    printf("\n Error in LU decomposition.\n\n");
    return 2;
  }
} // main

// end of file deter1.cpp