From 3a6dc04a7062a50f67f9c9cb56e7a042d0aff4c5 Mon Sep 17 00:00:00 2001
From: torset
Date: Fri, 20 Feb 2009 09:53:04 +0000
Subject: change prototypes of spec : if real on input, return real, not
 complex

---
 src/matrixOperations/spec2/dspec2a.c         | 69 +++++++++++++++-------------
 src/matrixOperations/spec2/sspec2a.c         | 12 ++---
 src/matrixOperations/spec2/testDoubleSpec2.c | 38 +++++----------
 src/matrixOperations/spec2/testFloatSpec2.c  | 22 ++-------
 4 files changed, 59 insertions(+), 82 deletions(-)

diff --git a/src/matrixOperations/spec2/dspec2a.c b/src/matrixOperations/spec2/dspec2a.c
index 47a4a445..ff4b0cb5 100644
--- a/src/matrixOperations/spec2/dspec2a.c
+++ b/src/matrixOperations/spec2/dspec2a.c
@@ -18,7 +18,7 @@
 #include "stdio.h"
 
 
-void dspec2a(double* in, int rows,doubleComplex* eigenvalues,doubleComplex* eigenvectors){
+void dspec2a(double* in, int rows,double* eigenvalues,double* eigenvectors){
 	int i=0, j=0, ij=0, ij1=0;
 	int symmetric=0;
 	int INFO=0;
@@ -30,9 +30,10 @@ void dspec2a(double* in, int rows,doubleComplex* eigenvalues,doubleComplex* eige
 	double* pdblRightvectors;
 	double* inCopy;
 	
+	/* FIXME : malloc here */
 	inCopy = malloc((uint)(rows*rows) * sizeof(double));
 	outReal = malloc((uint)rows * sizeof(double));
-	outImag = malloc((uint)rows * sizeof(double));
+	outImag = NULL;
 	pdblLeftvectors=NULL;
 	pdblRightvectors=NULL;
 	
@@ -59,52 +60,56 @@ void dspec2a(double* in, int rows,doubleComplex* eigenvalues,doubleComplex* eige
 	/* apply lapack function according to symmetry */
 	if(symmetric){	
 		C2F(dsyev)( "V", "U", &rows, inCopy, &rows, outReal, pdblWork, &iWorkSize, &INFO );
-		dzerosa(outImag,1,rows);
+	
+		/* Computation of eigenvectors */
+		for (i=0;i<rows*rows;i++) 	eigenvectors[i] = inCopy[i];
 	}
 	else {	
 		pdblRightvectors=malloc((uint)(rows*rows) * sizeof(double));
+		outImag = malloc((uint)rows * sizeof(double));
 		C2F(dgeev)( "N", "V", &rows, inCopy, &rows, outReal, outImag,
         			pdblLeftvectors, &rows, pdblRightvectors, &rows, pdblWork, &iWorkSize, &INFO );  
-	}
-	
-	/* Computation of eigenvalues */
-	zzerosa(eigenvalues,1,rows*rows);
-	for (i=0;i<rows;i++) eigenvalues[i+i*rows]=DoubleComplex(outReal[i],outImag[i]);
-	
 	
-	
-	/* Computation of eigenvectors */
-	j=0;
-	while (j<rows)
-	{
-		if (outImag[j]==0)
-		{	
-			for(i = 0 ; i < rows ; i++)
-			{
-				ij = i + j * rows;
-				if (symmetric)	eigenvectors[ij] = DoubleComplex(inCopy[ij],0);
-				else eigenvectors[ij] = DoubleComplex(pdblRightvectors[ij],0);
-			}
-			j = j + 1;
-		}
-		else
+		/* Computation of eigenvectors */        			
+        	j=0;
+		while (j<rows)
 		{
-			for(i = 0 ; i < rows ; i++)
+			if (outImag[j]==0)
+			{	
+				for(i = 0 ; i < rows ; i++)
+				{
+					ij = i + j * rows;
+					eigenvectors[ij] = pdblRightvectors[ij];
+				}
+				j = j + 1;
+			}
+			else
 			{
-				ij = i + j * rows;
-				ij1 = i + (j + 1) * rows;
-				eigenvectors[ij] = DoubleComplex(pdblRightvectors[ij],pdblRightvectors[ij1]);
-				eigenvectors[ij1] = DoubleComplex(pdblRightvectors[ij],-pdblRightvectors[ij1]);
+				for(i = 0 ; i < rows ; i++)
+				{
+					ij = i + j * rows;
+					ij1 = i + (j + 1) * rows;
+					eigenvectors[ij] = pdblRightvectors[ij];
+					eigenvectors[ij1] = pdblRightvectors[ij];
+				}
+				j = j + 2;
 			}
-			j = j + 2;
 		}
 	}
+	
+	/* Computation of eigenvalues */
+	dzerosa(eigenvalues,1,rows*rows);
+	for (i=0;i<rows;i++) eigenvalues[i+i*rows]=outReal[i];
+	
+
+	
 
 	free(inCopy);
 	free(outReal);
 	free(outImag);
-	free(pdblWork);
 	free(pdblLeftvectors);
 	free(pdblRightvectors);
+	free(pdblWork);
+
 }
 
diff --git a/src/matrixOperations/spec2/sspec2a.c b/src/matrixOperations/spec2/sspec2a.c
index 6b162b55..e168b540 100644
--- a/src/matrixOperations/spec2/sspec2a.c
+++ b/src/matrixOperations/spec2/sspec2a.c
@@ -15,22 +15,22 @@
 
 
 
-void sspec2a(float* in, int rows, floatComplex* eigenvalues,floatComplex* eigenvectors){
+void sspec2a(float* in, int rows, float* eigenvalues,float* eigenvectors){
 	/* As we use Lapack to find the eigenvalues, we must cast the float input into double
 	and the doubleComplex outputs of dspec2a into floatComplex*/
 	
 	int i;
 	double* dblin;
-	doubleComplex *dbleigenvalues,*dbleigenvectors;
+	double *dbleigenvalues,*dbleigenvectors;
 	
 	dblin=malloc((uint)(rows*rows)*sizeof(double));
-	dbleigenvalues=malloc((uint)(rows*rows)*sizeof(doubleComplex));	
-	dbleigenvectors=malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	dbleigenvalues=malloc((uint)(rows*rows)*sizeof(double));	
+	dbleigenvectors=malloc((uint)(rows*rows)*sizeof(double));
 	
 	for (i=0;i<rows*rows;i++) dblin[i]=(double)in[i];
 	
 	dspec2a(dblin,rows,dbleigenvalues,dbleigenvectors);
 	
-	for (i=0;i<rows*rows;i++) eigenvalues[i]=FloatComplex((float)zreals(dbleigenvalues[i]),(float)zimags(dbleigenvalues[i]));
-	for (i=0;i<rows*rows;i++) eigenvectors[i]=FloatComplex((float)zreals(dbleigenvectors[i]),(float)zimags(dbleigenvectors[i]));
+	for (i=0;i<rows*rows;i++) eigenvalues[i]=(float)dbleigenvalues[i];
+	for (i=0;i<rows*rows;i++) eigenvectors[i]=(float)dbleigenvectors[i];
 }
diff --git a/src/matrixOperations/spec2/testDoubleSpec2.c b/src/matrixOperations/spec2/testDoubleSpec2.c
index 0c9ce5f4..39550ca1 100644
--- a/src/matrixOperations/spec2/testDoubleSpec2.c
+++ b/src/matrixOperations/spec2/testDoubleSpec2.c
@@ -19,54 +19,42 @@
 static void dspec2aTest(void){
 	double in[4]={1,1,1,3};
 	double resultValuesR[4]={0.5857864376269050765700,0,0,3.4142135623730949234300};
-	double resultValuesI[4]={0,0,0,0};
-	double resultVectorsR[4]={- 0.9238795325112867384831,0.3826834323650897817792,
+		double resultVectorsR[4]={- 0.9238795325112867384831,0.3826834323650897817792,
 	0.3826834323650897817792,0.9238795325112867384831};
-	double resultVectorsI[4]={0,0,0,0};
-	
+		
 	
 	double in2[4]={1,1,-2,3};
 	double resultValues2R[4]={1.9999999999999997779554,0,0,1.9999999999999997779554};
-	double resultValues2I[4]={0.9999999999999997779554,0,0,-0.9999999999999997779554};
 	double resultVectors2R[4]={0.8164965809277261454824,- 0.4082482904638631282523,
 	0.8164965809277261454824,- 0.4082482904638631282523};
-	double resultVectors2I[4]={0,- 0.4082482904638629062077,0,0.4082482904638629062077};
 	
 	double in3[9]={0,-1,0,1,0,0,0,0,0};
 	double resultValues3R[9]={0};
-	double resultValues3I[9]={1,0,0,0,-1,0,0,0,0};
 	double resultVectors3R[9]={0.7071067811865474617150,0,0,0.7071067811865474617150,0,0,0,0,1};
-	double resultVectors3I[9]={0,0.7071067811865474617150,0,0,-0.7071067811865474617150,0,0,0,0};
+
 	
-	doubleComplex out1[4],out2[4],out3[9],out4[9];
+	double out1[4],out2[4],out3[9],out4[9];
 	
 	int i;
 	dspec2a(in3,3,out3,out4);
 	for(i=0;i<9;i++){
-		if (zreals(out3[i])>1e-16)  assert( fabs(zreals(out3[i])-resultValues3R[i]) / fabs(zreals(out3[i])) <1e-15);
-		else assert(1);
-		if (zimags(out3[i])>1e-16)  assert( fabs(zimags(out3[i])-resultValues3I[i]) / fabs(zimags(out3[i])) <1e-15);
+		if (out3[i]>1e-16)  assert( fabs(out3[i]-resultValues3R[i]) / fabs(out3[i]) <1e-15);
 		else assert(1);
 	}
 	for(i=0;i<9;i++){
-		if (zreals(out4[i])>1e-16)  assert( fabs(zreals(out4[i])-resultVectors3R[i]) / fabs(zreals(out4[i])) <1e-15);
-		else assert(1);
-		if (zimags(out4[i])>1e-16)  assert( fabs(zimags(out4[i])-resultVectors3I[i]) / fabs(zimags(out4[i])) <1e-15);
+		if (out4[i]>1e-16)  assert( fabs(out4[i]-resultVectors3R[i]) / fabs(out4[i]) <1e-15);
 		else assert(1);
 	}
 	
 
 	dspec2a(in,2,out1,out2);
+	for(i=0;i<4;i++) printf("%f\n",out1[i]);
 	for(i=0;i<4;i++){
-		if (zreals(out1[i])>1e-16)  assert( fabs(zreals(out1[i])-resultValuesR[i]) / fabs(zreals(out1[i])) <1e-15);
-		else assert(1);
-		if (zimags(out1[i])>1e-16)  assert( fabs(zimags(out1[i])-resultValuesI[i]) / fabs(zimags(out1[i])) <1e-16);
+		if (out1[i]>1e-16)  assert( fabs(out1[i]-resultValuesR[i]) / fabs(out1[i]) <1e-15);
 		else assert(1);
 	}
 	for(i=0;i<4;i++){
-		if (zreals(out2[i])>1e-16)  assert( fabs(zreals(out2[i])-resultVectorsR[i]) / fabs(zreals(out2[i])) <1e-15);
-		else assert(1);
-		if (zimags(out2[i])>1e-16)  assert( fabs(zimags(out2[i])-resultVectorsI[i]) / fabs(zimags(out2[i])) <1e-16);
+		if (out2[i]>1e-16)  assert( fabs(out2[i]-resultVectorsR[i]) / fabs(out2[i]) <1e-15);
 		else assert(1);
 	}
 	
@@ -74,15 +62,11 @@ static void dspec2aTest(void){
 	dspec2a(in2,2,out1,out2);
 
 	for(i=0;i<4;i++){
-		if (zreals(out1[i])>1e-16)  assert( fabs(zreals(out1[i])-resultValues2R[i]) / fabs(zreals(out1[i])) <1e-16);
-		else assert(1);
-		if (zimags(out1[i])>1e-16)  assert( fabs(zimags(out1[i])-resultValues2I[i]) / fabs(zimags(out1[i])) <1e-15);
+		if (out1[i]>1e-16)  assert( fabs(out1[i]-resultValues2R[i]) / fabs(out1[i]) <1e-16);
 		else assert(1);
 	}
 	for(i=0;i<4;i++){
-		if (zreals(out2[i])>1e-16)  assert( fabs(zreals(out2[i])-resultVectors2R[i]) / fabs(zreals(out2[i])) <1e-15);
-		else assert(1);
-		if (zimags(out2[i])>1e-16)  assert( fabs(zimags(out2[i])-resultVectors2I[i]) / fabs(zimags(out2[i])) <1e-15);
+		if (out2[i]>1e-16)  assert( fabs(out2[i]-resultVectors2R[i]) / fabs(out2[i]) <1e-15);
 		else assert(1);
 	}
 	
diff --git a/src/matrixOperations/spec2/testFloatSpec2.c b/src/matrixOperations/spec2/testFloatSpec2.c
index b8b4bbb6..fb2bc255 100644
--- a/src/matrixOperations/spec2/testFloatSpec2.c
+++ b/src/matrixOperations/spec2/testFloatSpec2.c
@@ -19,33 +19,25 @@
 static void sspec2aTest(void){
 	float in[4]={1.0f,1.0f,1.0f,3.0f};
 	float resultValuesR[4]={0.5857864376269050765700f,0,0,3.4142135623730949234300f};
-	float resultValuesI[4]={0.0f,0.0f,0.0f,0.0f};
 	float resultVectorsR[4]={- 0.9238795325112867384831f,0.3826834323650897817792f,
 	0.3826834323650897817792f,0.9238795325112867384831f};
-	float resultVectorsI[4]={0.0f,0.0f,0.0f,0.0f};
 	
 	float in2[4]={1.0f,1.0f,-2.0f,3.0f};
 	float resultValues2R[4]={1.9999999999999997779554f,0,0,1.9999999999999997779554f};
-	float resultValues2I[4]={0.9999999999999997779554f,0,0,-0.9999999999999997779554f};
 	float resultVectors2R[4]={0.8164965809277261454824f,- 0.4082482904638631282523f,
 	0.8164965809277261454824f,- 0.4082482904638631282523f};
-	float resultVectors2I[4]={0,- 0.4082482904638629062077f,0,0.4082482904638629062077f};
-	floatComplex eigenvalues[4],eigenvectors[4];
+	float eigenvalues[4],eigenvectors[4];
 
 	int i;
 		
 	sspec2a(in,2,eigenvalues,eigenvectors);
 
 	for(i=0;i<4;i++){
-		if (creals(eigenvalues[i])>1e-6)  assert( fabs(creals(eigenvalues[i])-resultValuesR[i]) / fabs(creals(eigenvalues[i])) <1e-16);
-		else assert(1);
-		if (cimags(eigenvalues[i])>1e-6)  assert( fabs(cimags(eigenvalues[i])-resultValuesI[i]) / fabs(cimags(eigenvalues[i])) <1e-16);
+		if (eigenvalues[i]>1e-6)  assert( fabs(eigenvalues[i]-resultValuesR[i]) / fabs(eigenvalues[i]) <1e-16);
 		else assert(1);
 	}
 	for(i=0;i<4;i++){
-		if (creals(eigenvectors[i])>1e-6)  assert( fabs(creals(eigenvectors[i])-resultVectorsR[i]) / fabs(creals(eigenvectors[i])) <1e-16);
-		else assert(1);
-		if (cimags(eigenvectors[i])>1e-6)  assert( fabs(cimags(eigenvectors[i])-resultVectorsI[i]) / fabs(cimags(eigenvectors[i])) <1e-16);
+		if (eigenvectors[i]>1e-6)  assert( fabs(eigenvectors[i]-resultVectorsR[i]) / fabs(eigenvectors[i]) <1e-16);
 		else assert(1);
 	}
 	
@@ -55,15 +47,11 @@ static void sspec2aTest(void){
 
 	
 	for(i=0;i<4;i++){
-		if (creals(eigenvalues[i])>1e-6)  assert( fabs(creals(eigenvalues[i])-resultValues2R[i]) / fabs(creals(eigenvalues[i])) <1e-16);
-		else assert(1);
-		if (cimags(eigenvalues[i])>1e-6)  assert( fabs(cimags(eigenvalues[i])-resultValues2I[i]) / fabs(cimags(eigenvalues[i])) <1e-16);
+		if (eigenvalues[i]>1e-6)  assert( fabs(eigenvalues[i]-resultValues2R[i]) / fabs(eigenvalues[i]) <1e-16);
 		else assert(1);
 	}
 	for(i=0;i<4;i++){
-		if (creals(eigenvectors[i])>1e-6)  assert( fabs(creals(eigenvectors[i])-resultVectors2R[i]) / fabs(creals(eigenvectors[i])) <1e-16);
-		else assert(1);
-		if (cimags(eigenvectors[i])>1e-6)  assert( fabs(cimags(eigenvectors[i])-resultVectors2I[i]) / fabs(cimags(eigenvectors[i])) <1e-16);
+		if (eigenvectors[i]>1e-6)  assert( fabs(eigenvectors[i]-resultVectors2R[i]) / fabs(eigenvectors[i]) <1e-16);
 		else assert(1);
 	}
 }
-- 
cgit