/*
 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 *  Copyright (C) 2008 - INRIA - Arnaud TORSET
 *
 *  This file must be used under the terms of the CeCILL.
 *  This source file is licensed as described in the file COPYING, which
 *  you should have received as part of this distribution.  The terms
 *  are also available at
 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 *
 */

#include <stdlib.h>
#ifndef WITHOUT_LAPACK
#include "lapack.h"
#else
#include "division.h"
#endif
#include "determ.h"
#include "multiplication.h"
#include "subtraction.h"
#include "addition.h"


doubleComplex zdeterma(doubleComplex * in, int size){
#ifndef WITHOUT_LAPACK
	doubleComplex tmp1,tmp2,out;
	int i=0,info=0;
	doubleComplex *inCopy;
	int* vectPivot;
	
	switch(size){
		case 2 : out = zdiffs(zmuls(in[0],in[3]),zmuls(in[1],in[2]));
			 break;
		case 3 : /*regle de Sarrus*/
			 out=DoubleComplex(0,0);
			 /*Addition Part*/
			 tmp1 = zmuls(in[0],in[4]);
			 tmp2 = zmuls(tmp1,in[8]);
			 out = zadds(out,tmp2);
			 tmp1 = zmuls(in[1],in[5]);
			 tmp2 = zmuls(tmp1,in[6]);
			 out = zadds(out,tmp2);
			 tmp1 = zmuls(in[2],in[3]);
			 tmp2 = zmuls(tmp1,in[7]);
			 out = zadds(out,tmp2);
			 /*Subtraction Part*/
			 tmp1 = zmuls(in[0],in[5]);
			 tmp2 = zmuls(tmp1,in[7]);
			 out = zdiffs(out,tmp2);
			 tmp1 = zmuls(in[1],in[3]);
			 tmp2 = zmuls(tmp1,in[8]);
			 out = zdiffs(out,tmp2);
			 tmp1 = zmuls(in[2],in[4]);
			 tmp2 = zmuls(tmp1,in[6]);
			 out = zdiffs(out,tmp2);
			 break;
			 
		default : 
			  /*Allocate inCopy and Copy in in inCopy*/
			  inCopy=(doubleComplex*)malloc((unsigned int)(size*size)*sizeof(doubleComplex));
			  for (i=0;i<size*size;i++) inCopy[i]=in[i];
			  
			  /*Calculation of determinant*/
			  vectPivot = (int*)malloc((unsigned int)size*sizeof(int));
			  zgetrf_ ( &size, &size, inCopy, &size, vectPivot, &info);
			  
			  out=DoubleComplex(1,0);
			  for (i=0;i<size;i++){
			  	if (vectPivot[i]!=i+1) out=DoubleComplex(-zreals(out),-zimags(out));
			  	out=zmuls(inCopy[i*(size+1)],out);
			  }
			  free(vectPivot);
			  free(inCopy);
			  break;
	}
	
#else
	int i=0,j=0,k=0;
	doubleComplex *inCopy;
	doubleComplex tmp1,tmp2;
	doubleComplex pivot;
	doubleComplex out;
	
		
	switch(size){
		case 2 : out = zdiffs(zmuls(in[0],in[3]),zmuls(in[1],in[2]));
			 break;
		case 3 : /*regle de Sarrus*/
			 out=DoubleComplex(0,0);
			 /*Addition Part*/
			 tmp1 = zmuls(in[0],in[4]);
			 tmp2 = zmuls(tmp1,in[8]);
			 out = zadds(out,tmp2);
			 tmp1 = zmuls(in[1],in[5]);
			 tmp2 = zmuls(tmp1,in[6]);
			 out = zadds(out,tmp2);
			 tmp1 = zmuls(in[2],in[3]);
			 tmp2 = zmuls(tmp1,in[7]);
			 out = zadds(out,tmp2);
			 /*Subtraction Part*/
			 tmp1 = zmuls(in[0],in[5]);
			 tmp2 = zmuls(tmp1,in[7]);
			 out = zdiffs(out,tmp2);
			 tmp1 = zmuls(in[1],in[3]);
			 tmp2 = zmuls(tmp1,in[8]);
			 out = zdiffs(out,tmp2);
			 tmp1 = zmuls(in[2],in[4]);
			 tmp2 = zmuls(tmp1,in[6]);
			 out = zdiffs(out,tmp2);
			 break;
			 
		default : 
				  /*Copy the input matrix*/
			  inCopy=malloc((unsigned int)(size*size)*sizeof(doubleComplex));
			  for (i=0;i<size*size;i++) inCopy[i]=in[i];
			  
    			  for (i=0;i<size;i++){
				for (j=i+1;j<size;j++){
					pivot = zrdivs(inCopy[i*size+j],inCopy[i*size+i]);
					for (k=0;k<size-i;k++){
						inCopy[i*size+j+k*size]=zdiffs(inCopy[i*size+j+k*size],zmuls(pivot,inCopy[i*size+i+k*size]));
					}
				}
			  }
			  out=DoubleComplex(1,0);
			  for (i=0;i<size;i++){
			  	out = zmuls(out,inCopy[i*size+i]);
			  }
			  free(inCopy);
			  break;
	}
#endif
	return out;
	
}