diff options
Diffstat (limited to '2.3-1/src/c/matrixOperations/determ/cdeterma.c')
| -rw-r--r-- | 2.3-1/src/c/matrixOperations/determ/cdeterma.c | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/2.3-1/src/c/matrixOperations/determ/cdeterma.c b/2.3-1/src/c/matrixOperations/determ/cdeterma.c new file mode 100644 index 00000000..6f2d1e8f --- /dev/null +++ b/2.3-1/src/c/matrixOperations/determ/cdeterma.c @@ -0,0 +1,150 @@ +/* + * 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" + + + +floatComplex cdeterma(floatComplex *in, int size){ +#ifndef WITHOUT_LAPACK + int i=0,info=0; + doubleComplex *inCopy,out_tmp; + int *vectPivot; + floatComplex out, tmp1,tmp2; + + + /*Calculation of the determinant*/ + switch(size){ + case 2 : out = cdiffs(cmuls(in[0],in[3]),cmuls(in[1],in[2])); + break; + case 3 : /*regle de Sarrus*/ + out=FloatComplex(0,0); + /*Addition part*/ + tmp1 = cmuls(in[0],in[4]); + tmp2 = cmuls(tmp1,in[8]); + out = cadds(out,tmp2); + tmp1 = cmuls(in[1],in[5]); + tmp2 = cmuls(tmp1,in[6]); + out = cadds(out,tmp2); + tmp1 = cmuls(in[2],in[3]); + tmp2 = cmuls(tmp1,in[7]); + out = cadds(out,tmp2); + /*Subtraction part*/ + tmp1 = cmuls(in[0],in[5]); + tmp2 = cmuls(tmp1,in[7]); + out = cdiffs(out,tmp2); + tmp1 = cmuls(in[1],in[3]); + tmp2 = cmuls(tmp1,in[8]); + out = cdiffs(out,tmp2); + tmp1 = cmuls(in[2],in[4]); + tmp2 = cmuls(tmp1,in[6]); + out = cdiffs(out,tmp2); + break; + + default : + + + /*We must have a doubleComplex matrix cause + zgetrf accept only double parameters*/ + + /*Copy the input matrix*/ + inCopy=(doubleComplex*)malloc((unsigned int)(size*size)*sizeof(doubleComplex)); + for (i=0;i<size*size;i++){ + inCopy[i]=DoubleComplex((double)creals(in[i]),(double)cimags(in[i])); + } + + + + vectPivot=(int*)malloc((unsigned int)size*sizeof(int)); + zgetrf_ ( &size, &size, inCopy, &size, vectPivot, &info); + + out_tmp=DoubleComplex(1,0); + for (i=0;i<size;i++){ + if (vectPivot[i]!=i+1) out_tmp=DoubleComplex(-zreals(out_tmp),-zimags(out_tmp)); + out_tmp=zmuls(inCopy[i*(size+1)],out_tmp); + } + out=FloatComplex((float)zreals(out_tmp),(float)zimags(out_tmp)); + + free(vectPivot); + free(inCopy); + } +#else + int i=0,j=0,k=0; + floatComplex *inCopy; + floatComplex tmp1,tmp2; + floatComplex pivot; + floatComplex out; + + + switch(size){ + case 2 : out = cdiffs(cmuls(in[0],in[3]),cmuls(in[1],in[2])); + break; + case 3 : /*regle de Sarrus*/ + out=FloatComplex(0,0); + /*Addition part*/ + tmp1 = cmuls(in[0],in[4]); + tmp2 = cmuls(tmp1,in[8]); + out = cadds(out,tmp2); + tmp1 = cmuls(in[1],in[5]); + tmp2 = cmuls(tmp1,in[6]); + out = cadds(out,tmp2); + tmp1 = cmuls(in[2],in[3]); + tmp2 = cmuls(tmp1,in[7]); + out = cadds(out,tmp2); + /*Subtraction part*/ + tmp1 = cmuls(in[0],in[5]); + tmp2 = cmuls(tmp1,in[7]); + out = cdiffs(out,tmp2); + tmp1 = cmuls(in[1],in[3]); + tmp2 = cmuls(tmp1,in[8]); + out = cdiffs(out,tmp2); + tmp1 = cmuls(in[2],in[4]); + tmp2 = cmuls(tmp1,in[6]); + out = cdiffs(out,tmp2); + break; + + default : + /*Copy the input matrix*/ + inCopy=malloc((unsigned int)(size*size)*sizeof(floatComplex)); + 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 = crdivs(inCopy[i*size+j],inCopy[i*size+i]); + for (k=0;k<size-i;k++){ + inCopy[i*size+j+k*size]=cdiffs(inCopy[i*size+j+k*size],cmuls(pivot,inCopy[i*size+i+k*size])); + } + } + } + out=FloatComplex(1,0); + for (i=0;i<size;i++){ + out = cmuls(out,inCopy[i*size+i]); + } + free(inCopy); + break; + } +#endif + + return out; +} + + |