change prototypes of matrixPow : if real on input, return real, not complex

3 files changed, 26 insertions, 19 deletions

diff --git a/src/matrixOperations/powm/dpowma.c b/src/matrixOperations/powm/dpowma.c
index 5d67189b..aecf3aac 100644
--- a/src/matrixOperations/powm/dpowma.c
+++ b/src/matrixOperations/powm/dpowma.c
@@ -15,19 +15,22 @@
 #include "pow.h"
 #include "matrixTranspose.h"
 #include "conj.h"
+#include "zeros.h"
 #include "matrixInversion.h"
 #include "matrixMultiplication.h"
 #include <stdio.h>
 
-void dpowma(double* in, int rows, double power, doubleComplex* out){
+void dpowma(double* in, int rows, double power, double* out){
 	int i=0, j=0;
 	int symmetric=0;
 	doubleComplex *eigenvalues,*eigenvectors,*tmp;
+	double* ZEROS;
 	
 	/* Data initialization */
 	eigenvalues = malloc((uint)(rows*rows)*sizeof(doubleComplex));
 	eigenvectors = malloc((uint)(rows*rows)*sizeof(doubleComplex));
 	tmp = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	ZEROS=malloc((uint)(rows*rows)*sizeof(double));
 	
 	/* symmetric test*/
 	for(i=0;i<rows;i++) {
@@ -40,8 +43,11 @@ void dpowma(double* in, int rows, double power, doubleComplex* out){
 	if ((i==rows)&&(j==rows)) symmetric=1;
 	
 	
+	
+	dzerosa(ZEROS,rows,rows);
+	tmp = DoubleComplexMatrix(in,ZEROS,rows*rows);
 	/* find eigenvalues and eigenvectors */
-	dspec2a(in, rows, eigenvalues,eigenvectors);
+	zspec2a(tmp, rows, eigenvalues,eigenvectors);
 /*	for (i=0;i<rows*rows;i++) printf("%f+%f*i\n",zreals(eigenvalues[i]),zimags(eigenvalues[i])); */
 	/* make operation on eigenvalues and eigenvectors */
 	for (i=0;i<rows;i++)
@@ -55,7 +61,9 @@ void dpowma(double* in, int rows, double power, doubleComplex* out){
 	}
 	else zinverma(eigenvectors, eigenvalues, rows);
 	
-	zmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	zmulma(tmp, rows, rows, eigenvalues, rows, rows, eigenvectors);
+	
+	for (i=0;i<rows*rows;i++) out[i]=zreals(eigenvectors[i]);
 	
 	free(eigenvalues);
 	free(eigenvectors);
diff --git a/src/matrixOperations/powm/spowma.c b/src/matrixOperations/powm/spowma.c
index 634cd449..eb48bda2 100644
--- a/src/matrixOperations/powm/spowma.c
+++ b/src/matrixOperations/powm/spowma.c
@@ -15,18 +15,21 @@
 #include "pow.h"
 #include "matrixTranspose.h"
 #include "conj.h"
+#include "zeros.h"
 #include "matrixInversion.h"
 #include "matrixMultiplication.h"
 
-void spowma(float* in, int rows, float power, floatComplex* out){
+void spowma(float* in, int rows, float power, float* out){
 	int i=0, j=0;
 	int symmetric=0;
 	floatComplex *eigenvalues,*eigenvectors,*tmp;
+	float* ZEROS;
 	
 	/* Data initialization */
 	eigenvalues = malloc((uint)(rows*rows)*sizeof(floatComplex));
 	eigenvectors = malloc((uint)(rows*rows)*sizeof(floatComplex));
 	tmp = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	ZEROS = malloc((uint)(rows*rows)*sizeof(float));
 	
 	/* symmetric test*/
 	for(i=0;i<rows;i++) {
@@ -38,9 +41,11 @@ void spowma(float* in, int rows, float power, floatComplex* out){
 	
 	if ((i==rows)&&(j==rows)) symmetric=1;
 	
+	szerosa(ZEROS,rows,rows);
+	tmp = FloatComplexMatrix(in,ZEROS,rows*rows);
 	
 	/* find eigenvalues and eigenvectors */
-	sspec2a(in, rows, eigenvalues,eigenvectors);
+	cspec2a(tmp, rows, eigenvalues,eigenvectors);
 	
 	/* make operation on eigenvalues and eigenvectors */
 	for (i=0;i<rows;i++)
@@ -54,7 +59,10 @@ void spowma(float* in, int rows, float power, floatComplex* out){
 	}
 	else cinverma(eigenvectors, eigenvalues, rows);
 	
-	cmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	cmulma(tmp, rows, rows, eigenvalues, rows, rows, eigenvectors);
+	
+	for (i=0;i<rows*rows;i++) out[i]=creals(eigenvectors[i]);
+	
 	
 	free(eigenvalues);
 	free(eigenvectors);
diff --git a/src/matrixOperations/powm/testDoublePowm.c b/src/matrixOperations/powm/testDoublePowm.c
index 05cad5e5..1eb5a5e6 100644
--- a/src/matrixOperations/powm/testDoublePowm.c
+++ b/src/matrixOperations/powm/testDoublePowm.c
@@ -21,9 +21,7 @@ static void dpowmaTest(void){
 	double expand1=2.2;
 	double result1R[4]={  27.93459280052221771484 ,    23.580294119266994812278  ,
 	18.864235295413593007652 ,   32.650651624375619519469   };
-	double result1I[4]={  3.6611113731522362257920 , - 3.6611113731522362257920 , 
-  - 2.9288890985217883589087  ,  2.9288890985217883589087 };
-	doubleComplex out1[4];
+	double out1[4];
 	int i;
 	
 	double in2[16]={    2.5358983855694532394409 ,   9.0725262500345706939697,    0.0026536155492067337036,    3.9639251008629798889160 , 
@@ -35,28 +33,21 @@ static void dpowmaTest(void){
 	24728.411825244897045195 ,   18392.823733925368287601  ,  18631.05868385956637212   ,  19357.84707477861229563   ,
 	16169.682243927050876664 ,   12258.542785024719705689  ,  12630.164466338968850323  ,  12827.915677254180991440  ,
 	13742.841851328515986097 ,   10198.0420642120679986    ,  10658.784670951883526868  ,  10839.51135004585739807  };
-	double result2I[16]={ - 7.1981835972120027378196  ,  1.9386514637886893552832,  - 17.692616672339234185074  ,  24.561537532538231687340 , 
-	- 2.2418859631076406557781 ,   0.6037961445855435371755,  - 5.5103941755046683681485,    7.649730724813480264857  , 
-	- 4.865855522250573272913  ,   1.310496989059492634056 ,  - 11.95992230200565309417 ,    16.603201547139228466676  ,
-	16.00935601900000193609 ,  - 4.3117212921047043394651  ,  39.34984366402868971591  , - 54.626892107189902958453 };
-	doubleComplex out2[16];
+	double out2[16];
 	
 	
 	dpowma(in1, 2, expand1, out1);
 	dpowma(in2, 4, expand2, out2);
 	
 	for (i=0;i<4;i++) {
-		assert( fabs(zreals(out1[i])-result1R[i]) / fabs(zreals(out1[i])) <1e-15);
-		assert( fabs(zimags(out1[i])-result1I[i]) / fabs(zimags(out1[i])) <1e-15);
+		assert( fabs(out1[i]-result1R[i]) / fabs(out1[i]) <1e-15);
 	}
 	
 /*	
 	FIXME : assert 1e-11 maybe due to spec2
 */	
 	for (i=0;i<16;i++) {
-		printf("out[%d] = %1.16f+%1.16f*i --- result = %1.16f+%1.16f*i\n",i,zreals(out2[i]),zimags(out2[i]),result2R[i],result2I[i]);		
-		assert( fabs(zreals(out2[i])-result2R[i]) / fabs(zreals(out2[i])) <1e-14);
-		assert( fabs(zimags(out2[i])-result2I[i]) / fabs(zimags(out2[i])) <1e-11);
+		assert( fabs(out2[i]-result2R[i]) / fabs(out2[i]) <1e-14);
 	}	
 	
 }