Fixed errors in spec2

3 files changed, 39 insertions, 23 deletions

diff --git a/src/matrixOperations/includes/spec.h b/src/matrixOperations/includes/spec.h
index 0ab5e4b6..d00978e3 100644
--- a/src/matrixOperations/includes/spec.h
+++ b/src/matrixOperations/includes/spec.h
@@ -21,8 +21,8 @@
    If you want the eigenvalues and the eigenvectors, use spec2 */
 
 /* spec */
-#define sspecs(in) 	in
-#define dspecs(in) 	in
+#define sspecs(in) 	FloatComplex(in,0)
+#define dspecs(in) 	DoubleComplex(in,0)
 #define cspecs(in) 	in
 #define zspecs(in) 	in
 
@@ -34,10 +34,10 @@ void zspeca(doubleComplex* in, int rows,doubleComplex* out);
 
 
 /* spec2 */
-#define sspec2s(in,out)  	{out=in;return 1;}
-#define dspec2s(in,out)  	{out=in;return 1;}
-#define cspec2s(in,out)   	{out=in;return FloatComplex(1,0);}
-#define zspec2s(in,out)  	{out=in;return DoubleCompelx(1,0);}
+#define sspec2s(in,out)  	sspecs(1);*out=FloatComplex(in,0); 
+#define dspec2s(in,out)  	dspecs(1);*out=DoubleComplex(in,0); 
+#define cspec2s(in,out)   	cspecs(FloatComplex(1,0));*out=FloatComplex(creals(in),cimags(in)); 
+#define zspec2s(in,out)  	zspecs(DoubleComplex(1,0));*out=DoubleComplex(zreals(in),zimags(in)); 
 
 void sspec2a(float* in, int rows, floatComplex* eigenvalues,floatComplex* eigenvectors);
 void dspec2a(double* in, int rows, doubleComplex* eigenvalues,doubleComplex* eigenvectors);
diff --git a/src/matrixOperations/interfaces/int_spec.h b/src/matrixOperations/interfaces/int_spec.h
index b0293e66..a2d9e840 100644
--- a/src/matrixOperations/interfaces/int_spec.h
+++ b/src/matrixOperations/interfaces/int_spec.h
@@ -17,9 +17,9 @@
 
 /* Only eigenvalues */
 
-#define s0specs0(in)			sspecs(in)
+#define s0specc0(in)			sspecs(in)
 
-#define d0specd0(in)			dspecs(in)
+#define d0specz0(in)			dspecs(in)
 
 #define c0specc0(in)			cspecs(in)
 
@@ -35,20 +35,20 @@
 
 /* Eigenvalues and eigenvectors */
 
-#define s0specs0s0(in,out)			sspec2s(in,out)
+#define s0specc0c0(in,out)			sspec2s(in,out)
 
-#define d0specd0d0(in,out)			dspec2s(in,out)
+#define d0specz0z0(in,out)			dspec2s(in,out)
 
 #define c0specc0c0(in,out)			cspec2s(in,out)
 
 #define z0specz0z0(in,out)			zspec2s(in,out)
 
-#define s2specc2c2(in,size,out)		sspeca(in, size[0], out)
+#define s2specc2c2(in,size,out1,out2)		sspec2a(in, size[0], out2, out1)
 
-#define d2specz2z2(in,size,out)		dspeca(in, size[0], out)
+#define d2specz2z2(in,size,out1,out2)		dspec2a(in, size[0], out2, out1)
 
-#define c2specc2c2(in,size,out)		cspeca(in, size[0], out)
+#define c2specc2c2(in,size,out1,out2)		cspec2a(in, size[0], out2, out1)
 
-#define z2specz2z2(in,size,out)		zspeca(in, size[0], out)
+#define z2specz2z2(in,size,out1,out2)		zspec2a(in, size[0], out2, out1)
 
 #endif /* !__INT_SPEC_H__ */
diff --git a/src/matrixOperations/spec2/dspec2a.c b/src/matrixOperations/spec2/dspec2a.c
index 630cf55e..47a4a445 100644
--- a/src/matrixOperations/spec2/dspec2a.c
+++ b/src/matrixOperations/spec2/dspec2a.c
@@ -59,15 +59,36 @@ 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 );
-		for (i=0;i<rows*rows;i++) eigenvectors[i]=DoubleComplex(inCopy[i],0);
 		dzerosa(outImag,1,rows);
 	}
 	else {	
 		pdblRightvectors=malloc((uint)(rows*rows) * sizeof(double));
 		C2F(dgeev)( "N", "V", &rows, inCopy, &rows, outReal, outImag,
         			pdblLeftvectors, &rows, pdblRightvectors, &rows, pdblWork, &iWorkSize, &INFO );  
-        	      	
-        	for(j = 0 ; j < rows ; j+=2){
+	}
+	
+	/* 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
+		{
 			for(i = 0 ; i < rows ; i++)
 			{
 				ij = i + j * rows;
@@ -75,14 +96,9 @@ void dspec2a(double* in, int rows,doubleComplex* eigenvalues,doubleComplex* eige
 				eigenvectors[ij] = DoubleComplex(pdblRightvectors[ij],pdblRightvectors[ij1]);
 				eigenvectors[ij1] = DoubleComplex(pdblRightvectors[ij],-pdblRightvectors[ij1]);
 			}
+			j = j + 2;
 		}
-		
 	}
-	
-	
-	zzerosa(eigenvalues,1,rows*rows);
-	for (i=0;i<rows;i++) eigenvalues[i+i*rows]=DoubleComplex(outReal[i],outImag[i]);
-
 
 	free(inCopy);
 	free(outReal);