From 92d3899c14d576c20762d7860f5e2f137923da0c Mon Sep 17 00:00:00 2001
From: torset
Date: Mon, 16 Feb 2009 15:29:51 +0000
Subject: Miscellanous change - Major change : change size by rows for more
 comprehension

---
 src/matrixOperations/logm/zlogma.c | 83 +++++++++++++++++++++++---------------
 1 file changed, 51 insertions(+), 32 deletions(-)

(limited to 'src/matrixOperations/logm/zlogma.c')

diff --git a/src/matrixOperations/logm/zlogma.c b/src/matrixOperations/logm/zlogma.c
index 37516822..5a10151e 100644
--- a/src/matrixOperations/logm/zlogma.c
+++ b/src/matrixOperations/logm/zlogma.c
@@ -19,8 +19,9 @@
 #include "matrixTranspose.h"
 #include "matrixInversion.h"
 #include "max.h"
+#include "conj.h"
 
-void zlogma (doubleComplex* in, int size, doubleComplex* out){
+void zlogma (doubleComplex* in, int rows, doubleComplex* out){
 
 	/* Algo : */
 	/* trouver les valeurs propres vp */
@@ -31,54 +32,67 @@ void zlogma (doubleComplex* in, int size, doubleComplex* out){
 	
 	int i;
 	doubleComplex *eigenvalues, *eigenvectors, *tmp;
-	int symetrique = 0;
+	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 */
-	doubleComplex * inCopy = NULL; /* Used by LAPACK */
 	
 	
 	
 	
 	/* Data initialization */
-	eigenvalues = malloc((uint)(size*size)*sizeof(doubleComplex));
-	eigenvectors = malloc((uint)(size*size)*sizeof(doubleComplex));
-	tmp = malloc((uint)(size*size)*sizeof(doubleComplex));
-	iWorkSize = max(1,2*size);
+	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*size)*sizeof(doubleComplex));
-	inCopy = malloc((uint)(size*size)*sizeof(doubleComplex));
+	pdblRWork = malloc((uint)(2*rows)*sizeof(doubleComplex));
 	
+
 	/* Copy in in inCopy */
-	for (i=0;i<size*size;i++)
-		inCopy[i] = in[i];
+	for (i=0;i<rows*rows;i++)
+		out[i] = in[i];
 
 	
 	
 	
 	
-	/* regarde si in est symetrique */
-	ztransposea(in,size,size,tmp);
-	for (i=0;i<size*size;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;
 	}
-	if (i==size*size) symetrique = 1;
+	if (i==rows*rows) hermitienne = 1;
 		
 	/* trouver les valeurs propres vp ainsi que les vecteurs propres*/
-	if (symetrique){
-		C2F(zheev)( "V", "U", &size, inCopy, &size, eigenvalues,
+	if (hermitienne){
+	printf("wazaaaaaaa\n");
+		C2F(zheev)( "V", "U", &rows, out, &rows, eigenvalues,
 			    pdblWork, &iWorkSize, pdblRWork, &info );
-		eigenvectors = inCopy;
+		/* 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 {	
-		C2F(zgeev)( "N", "V", &size, inCopy, &size, eigenvalues, 
-			    pdblLeftvectors, &size, eigenvectors, &size, pdblWork, &iWorkSize,
+	printf("passe\n");
+		C2F(zgeev)( "N", "V", &rows, out, &rows, eigenvalues, 
+			    pdblLeftvectors, &rows, eigenvectors, &rows, pdblWork, &iWorkSize,
 			    pdblRWork, &info );
-		for(i=0;i<size;i++) {
-			eigenvalues[i*size+i] = eigenvalues[i];
-			if ((i*size+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
+		for(i=0;i<rows;i++) {
+			eigenvalues[i*rows+i] = eigenvalues[i];
+			if ((i*rows+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
 		}
 	}
 
@@ -88,28 +102,33 @@ void zlogma (doubleComplex* in, int size, doubleComplex* out){
 	
 	
 	/* diag(log(diag(vp)) */
-	for (i=0;i<size*size;i++){
-		if ((i%(size+1))==0) /* teste si i est sur la diagonale */
+	for (i=0;i<rows*rows;i++){
+		if ((i%(rows+1))==0) /* teste si i est sur la diagonale */
 			eigenvalues[i] = zlogs(eigenvalues[i]);
 		else    eigenvalues[i] = DoubleComplex(0,0);
 	}
 	
 	/* Vp * diag(log(diag(vp)) */
-	zmulma(eigenvectors, size, size, eigenvalues, size, size, tmp);
+	zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
 	
-	/* inv(Vp) */
-	zinverma(eigenvectors, eigenvectors, size);
+	/* Vp' ou inv(Vp) */
+	if (hermitienne) {
+		/* we use eigenvalues as a temporary matrix cause it's useless now*/
+		ztransposea(eigenvectors,rows,rows,eigenvalues);
+		zconja(eigenvalues,rows*rows,eigenvectors);
+	}
+	else	zinverma(eigenvectors, eigenvectors, rows);
 	
 	/* Vp * diag(log(diag(vp))*inv(Vp) */
-	zmulma(tmp, size, size, eigenvectors, size, size, out);
+	zmulma(tmp, rows, rows, eigenvectors, rows, rows, out);
 	
 	
+	/* FIXME : pb with eigenvectors */
+	eigenvectors=NULL;
 	free(eigenvalues);
-	free(eigenvectors);
+	free(eigenvectors);		
 	free(tmp);
 	free(pdblWork);
 	free(pdblRWork);
-	free(inCopy);
-	
 
 }
-- 
cgit