zlogma use now spec2

3 files changed, 54 insertions, 93 deletions

diff --git a/src/matrixOperations/logm/Makefile.am b/src/matrixOperations/logm/Makefile.am
index 85e8d7f2..cb69e8f0 100644
--- a/src/matrixOperations/logm/Makefile.am
+++ b/src/matrixOperations/logm/Makefile.am
@@ -42,23 +42,25 @@ libLogm_la_SOURCES = 	$(HEAD) \
 check_PROGRAMS = 	testDoubleLogm
 
 check_LDADD =		$(top_builddir)/type/libDoubleComplex.la \
-			$(top_builddir)/type/libFloatComplex.la \
-			$(top_builddir)/lib/lapack/libscilapack.la \
-			$(top_builddir)/lib/blas/libsciblas.la \
-			$(top_builddir)/elementaryFunctions/log/libLog.la \
-			$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
-			$(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
-			$(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
-			$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
-			$(top_builddir)/auxiliaryFunctions/conj/libConj.la \
-			$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
-			$(top_builddir)/auxiliaryFunctions/sign/libSign.la \
-			$(top_builddir)/operations/multiplication/libMultiplication.la \
-			$(top_builddir)/operations/addition/libAddition.la \
-			$(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
-			$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
-			$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
-			$(top_builddir)/matrixOperations/logm/libLogm.la
+				$(top_builddir)/type/libFloatComplex.la \
+				$(top_builddir)/lib/lapack/libscilapack.la \
+				$(top_builddir)/lib/blas/libsciblas.la \
+				$(top_builddir)/elementaryFunctions/log/libLog.la \
+				$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+				$(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
+				$(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
+				$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+				$(top_builddir)/auxiliaryFunctions/conj/libConj.la \
+				$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
+				$(top_builddir)/auxiliaryFunctions/sign/libSign.la \
+				$(top_builddir)/operations/multiplication/libMultiplication.la \
+				$(top_builddir)/operations/addition/libAddition.la \
+				$(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
+				$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
+				$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
+				$(top_builddir)/matrixOperations/spec2/libSpec2.la \
+				$(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
+				$(top_builddir)/matrixOperations/logm/libLogm.la
 
 			
 check_INCLUDES =	-I ./includes \
@@ -66,8 +68,7 @@ check_INCLUDES =	-I ./includes \
 			-I $(top_builddir)/elementaryFunctions/includes \
 			-I $(top_builddir)/operations/includes \
 			-I $(top_builddir)/auxiliaryFunctions/includes \
-			-I $(top_builddir)/matrixOperations/includes\
-			-I $(top_builddir)/matrixOperations/logm
+			-I $(top_builddir)/matrixOperations/includes
 
 testDoubleLogm_SOURCES = testDoubleLogm.c
 testDoubleLogm_LDADD = $(check_LDADD)
diff --git a/src/matrixOperations/logm/Makefile.in b/src/matrixOperations/logm/Makefile.in
index 93b2dc6f..cb74434a 100644
--- a/src/matrixOperations/logm/Makefile.in
+++ b/src/matrixOperations/logm/Makefile.in
@@ -210,31 +210,32 @@ libLogm_la_SOURCES = $(HEAD) \
 				 dlogma.c 
 
 check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
-			$(top_builddir)/type/libFloatComplex.la \
-			$(top_builddir)/lib/lapack/libscilapack.la \
-			$(top_builddir)/lib/blas/libsciblas.la \
-			$(top_builddir)/elementaryFunctions/log/libLog.la \
-			$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
-			$(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
-			$(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
-			$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
-			$(top_builddir)/auxiliaryFunctions/conj/libConj.la \
-			$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
-			$(top_builddir)/auxiliaryFunctions/sign/libSign.la \
-			$(top_builddir)/operations/multiplication/libMultiplication.la \
-			$(top_builddir)/operations/addition/libAddition.la \
-			$(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
-			$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
-			$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
-			$(top_builddir)/matrixOperations/logm/libLogm.la
+				$(top_builddir)/type/libFloatComplex.la \
+				$(top_builddir)/lib/lapack/libscilapack.la \
+				$(top_builddir)/lib/blas/libsciblas.la \
+				$(top_builddir)/elementaryFunctions/log/libLog.la \
+				$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+				$(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
+				$(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
+				$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+				$(top_builddir)/auxiliaryFunctions/conj/libConj.la \
+				$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
+				$(top_builddir)/auxiliaryFunctions/sign/libSign.la \
+				$(top_builddir)/operations/multiplication/libMultiplication.la \
+				$(top_builddir)/operations/addition/libAddition.la \
+				$(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
+				$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
+				$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
+				$(top_builddir)/matrixOperations/spec2/libSpec2.la \
+				$(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
+				$(top_builddir)/matrixOperations/logm/libLogm.la
 
 check_INCLUDES = -I ./includes \
 			-I $(top_builddir)/type \
 			-I $(top_builddir)/elementaryFunctions/includes \
 			-I $(top_builddir)/operations/includes \
 			-I $(top_builddir)/auxiliaryFunctions/includes \
-			-I $(top_builddir)/matrixOperations/includes\
-			-I $(top_builddir)/matrixOperations/logm
+			-I $(top_builddir)/matrixOperations/includes
 
 testDoubleLogm_SOURCES = testDoubleLogm.c
 testDoubleLogm_LDADD = $(check_LDADD)
diff --git a/src/matrixOperations/logm/zlogma.c b/src/matrixOperations/logm/zlogma.c
index 5a10151e..3a629bfd 100644
--- a/src/matrixOperations/logm/zlogma.c
+++ b/src/matrixOperations/logm/zlogma.c
@@ -20,6 +20,7 @@
 #include "matrixInversion.h"
 #include "max.h"
 #include "conj.h"
+#include "spec.h"
 
 void zlogma (doubleComplex* in, int rows, doubleComplex* out){
 
@@ -30,77 +31,37 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
 	*	logm = Vp * diag(log(diag(vp)) * inv(Vp) */
 	
 	
-	int i;
+	int i = 0,j = 0;
 	doubleComplex *eigenvalues, *eigenvectors, *tmp;
 	int hermitienne = 0;
-	int info = 0; /* Used by LAPACK */
-	int iWorkSize = 0; /* Used by LAPACK */
-	doubleComplex *pdblWork = NULL; /* Used by LAPACK */
-	doubleComplex *pdblRWork = NULL; /* Used by LAPACK */
-	doubleComplex * pdblLeftvectors = NULL; /* Used by LAPACK */
-	
-	
-	
+		
 	
 	/* Data initialization */
 	eigenvalues = malloc((uint)(rows*rows)*sizeof(doubleComplex));
 	eigenvectors = malloc((uint)(rows*rows)*sizeof(doubleComplex));
 	tmp = malloc((uint)(rows*rows)*sizeof(doubleComplex));
-	iWorkSize = max(1,2*rows);
-	pdblWork = malloc((uint)(iWorkSize)*sizeof(doubleComplex));
-	pdblRWork = malloc((uint)(2*rows)*sizeof(doubleComplex));
 	
-
-	/* Copy in in inCopy */
-	for (i=0;i<rows*rows;i++)
-		out[i] = in[i];
-
 	
 	
 	
 	
 	/* regarde si in est hermitienne */
-	ztransposea(in,rows,rows,tmp);
-	zconja(tmp,rows*rows,tmp);
-	for (i=0;i<rows*rows;i++){
-		if  (  (zreals(in[i])!=zreals(tmp[i])) || (zimags(in[i])!=zimags(tmp[i])) ) break;
+	for (i=0;i<rows;i++){
+		for(j=0;j<rows;j++)
+			if ( (zreals(in[i+j*rows])!=zreals(in[j+i*rows])) ||  
+			     (zimags(in[i+j*rows])!=zimags(in[j+i*rows])) )  break;
+		if (j!=rows) break;
 	}
-	if (i==rows*rows) hermitienne = 1;
+	if ((i==rows) && (j==rows)) hermitienne=1;
 		
 	/* trouver les valeurs propres vp ainsi que les vecteurs propres*/
-	if (hermitienne){
-	printf("wazaaaaaaa\n");
-		C2F(zheev)( "V", "U", &rows, out, &rows, eigenvalues,
-			    pdblWork, &iWorkSize, pdblRWork, &info );
-		/* stocker les valeurs propres comme il faut parce qu'elle sont mal rangées
-		   une fois sorties de zheev*/
-		for(i=rows-1;i>=0;i--) {
-			/* on le fait en descendant sinon on écrase l'element 0 alors qu'il doit etre réutilisé
-				---> erreur */
-			if ( ((i+i*rows)/rows)%2 == 0)
-				eigenvalues[i*rows+i] = DoubleComplex(zreals(eigenvalues[(i*rows+i)/(2*rows)]),0);
-			else 	eigenvalues[i*rows+i] = DoubleComplex(zimags(eigenvalues[(i*rows+i)/(2*rows)]),0);
-		}
-		for (i=1;i<rows;i++) eigenvalues[i]=DoubleComplex(0,0);
-		eigenvectors = out;
-
-	}
-	else {	
-	printf("passe\n");
-		C2F(zgeev)( "N", "V", &rows, out, &rows, eigenvalues, 
-			    pdblLeftvectors, &rows, eigenvectors, &rows, pdblWork, &iWorkSize,
-			    pdblRWork, &info );
-		for(i=0;i<rows;i++) {
-			eigenvalues[i*rows+i] = eigenvalues[i];
-			if ((i*rows+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
-		}
-	}
+	zspec2a(in,rows,eigenvalues,eigenvectors);
 
 
 	/* utiliser la formule suivante
 	*	logm = Vp * diag(log(diag(vp)) * inv(Vp) */
 	
-	
+
 	/* diag(log(diag(vp)) */
 	for (i=0;i<rows*rows;i++){
 		if ((i%(rows+1))==0) /* teste si i est sur la diagonale */
@@ -108,6 +69,8 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
 		else    eigenvalues[i] = DoubleComplex(0,0);
 	}
 	
+	
+	
 	/* Vp * diag(log(diag(vp)) */
 	zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
 	
@@ -122,13 +85,9 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
 	/* Vp * diag(log(diag(vp))*inv(Vp) */
 	zmulma(tmp, rows, rows, eigenvectors, rows, rows, out);
 	
-	
-	/* FIXME : pb with eigenvectors */
-	eigenvectors=NULL;
+
 	free(eigenvalues);
 	free(eigenvectors);		
 	free(tmp);
-	free(pdblWork);
-	free(pdblRWork);
 
 }