7 files changed, 341 insertions, 131 deletions
diff --git a/src/matrixOperations/powm/Makefile.am b/src/matrixOperations/powm/Makefile.am
index 4215c21e..d3ce9bc7 100644
--- a/src/matrixOperations/powm/Makefile.am
+++ b/src/matrixOperations/powm/Makefile.am
@@ -30,7 +30,7 @@ libMatrixPow_la_SOURCES = 	$(HEAD) \
 				cpowma.c\
 				zpowma.c
 
-check_PROGRAMS = 	testDoubleMatrixPow
+check_PROGRAMS = 	testDoubleMatrixPow	testFloatMatrixPow
 
 check_LDADD =	$(top_builddir)/type/libDoubleComplex.la \
 			$(top_builddir)/type/libFloatComplex.la \
@@ -50,6 +50,8 @@ check_LDADD =	$(top_builddir)/type/libDoubleComplex.la \
 			$(top_builddir)/elementaryFunctions/cosh/libCosh.la \
 			$(top_builddir)/elementaryFunctions/sin/libSin.la \
 			$(top_builddir)/elementaryFunctions/sinh/libSinh.la \
+			$(top_builddir)/matrixOperations/spec2/libSpec2.la \
+			$(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
 			$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
 			$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
 			$(top_builddir)/matrixOperations/cat/libMatrixConcatenation.la \
@@ -77,4 +79,8 @@ testDoubleMatrixPow_SOURCES = testDoublePowm.c
 testDoubleMatrixPow_LDADD = $(check_LDADD)
 testDoubleMatrixPow_CFLAGS = $(check_INCLUDES)
 
-TESTS =		testDoubleMatrixPow
+testFloatMatrixPow_SOURCES = testFloatPowm.c
+testFloatMatrixPow_LDADD = $(check_LDADD)
+testFloatMatrixPow_CFLAGS = $(check_INCLUDES)
+
+TESTS =	testDoubleMatrixPow	testFloatMatrixPow
diff --git a/src/matrixOperations/powm/Makefile.in b/src/matrixOperations/powm/Makefile.in
index 32b7f77b..0f801bb7 100644
--- a/src/matrixOperations/powm/Makefile.in
+++ b/src/matrixOperations/powm/Makefile.in
@@ -32,8 +32,9 @@ PRE_UNINSTALL = :
 POST_UNINSTALL = :
 build_triplet = @build@
 host_triplet = @host@
-check_PROGRAMS = testDoubleMatrixPow$(EXEEXT)
-TESTS = testDoubleMatrixPow$(EXEEXT)
+check_PROGRAMS = testDoubleMatrixPow$(EXEEXT) \
+	testFloatMatrixPow$(EXEEXT)
+TESTS = testDoubleMatrixPow$(EXEEXT) testFloatMatrixPow$(EXEEXT)
 subdir = matrixOperations/powm
 DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -69,6 +70,14 @@ testDoubleMatrixPow_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) \
 	$(testDoubleMatrixPow_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
 	$(LDFLAGS) -o $@
+am_testFloatMatrixPow_OBJECTS =  \
+	testFloatMatrixPow-testFloatPowm.$(OBJEXT)
+testFloatMatrixPow_OBJECTS = $(am_testFloatMatrixPow_OBJECTS)
+testFloatMatrixPow_DEPENDENCIES = $(check_LDADD)
+testFloatMatrixPow_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
+	$(LIBTOOLFLAGS) --mode=link $(CCLD) \
+	$(testFloatMatrixPow_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
+	$(LDFLAGS) -o $@
 DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)/includes
 depcomp = $(SHELL) $(top_srcdir)/config/depcomp
 am__depfiles_maybe = depfiles
@@ -81,9 +90,10 @@ CCLD = $(CC)
 LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
 	--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
 	$(LDFLAGS) -o $@
-SOURCES = $(libMatrixPow_la_SOURCES) $(testDoubleMatrixPow_SOURCES)
+SOURCES = $(libMatrixPow_la_SOURCES) $(testDoubleMatrixPow_SOURCES) \
+	$(testFloatMatrixPow_SOURCES)
 DIST_SOURCES = $(libMatrixPow_la_SOURCES) \
-	$(testDoubleMatrixPow_SOURCES)
+	$(testDoubleMatrixPow_SOURCES) $(testFloatMatrixPow_SOURCES)
 ETAGS = etags
 CTAGS = ctags
 DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
@@ -231,6 +241,8 @@ check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
 			$(top_builddir)/elementaryFunctions/cosh/libCosh.la \
 			$(top_builddir)/elementaryFunctions/sin/libSin.la \
 			$(top_builddir)/elementaryFunctions/sinh/libSinh.la \
+			$(top_builddir)/matrixOperations/spec2/libSpec2.la \
+			$(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
 			$(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
 			$(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
 			$(top_builddir)/matrixOperations/cat/libMatrixConcatenation.la \
@@ -257,6 +269,9 @@ check_INCLUDES = -I $(top_builddir)/type \
 testDoubleMatrixPow_SOURCES = testDoublePowm.c
 testDoubleMatrixPow_LDADD = $(check_LDADD)
 testDoubleMatrixPow_CFLAGS = $(check_INCLUDES)
+testFloatMatrixPow_SOURCES = testFloatPowm.c
+testFloatMatrixPow_LDADD = $(check_LDADD)
+testFloatMatrixPow_CFLAGS = $(check_INCLUDES)
 all: all-am
 
 .SUFFIXES:
@@ -329,6 +344,9 @@ clean-checkPROGRAMS:
 testDoubleMatrixPow$(EXEEXT): $(testDoubleMatrixPow_OBJECTS) $(testDoubleMatrixPow_DEPENDENCIES) 
 	@rm -f testDoubleMatrixPow$(EXEEXT)
 	$(testDoubleMatrixPow_LINK) $(testDoubleMatrixPow_OBJECTS) $(testDoubleMatrixPow_LDADD) $(LIBS)
+testFloatMatrixPow$(EXEEXT): $(testFloatMatrixPow_OBJECTS) $(testFloatMatrixPow_DEPENDENCIES) 
+	@rm -f testFloatMatrixPow$(EXEEXT)
+	$(testFloatMatrixPow_LINK) $(testFloatMatrixPow_OBJECTS) $(testFloatMatrixPow_LDADD) $(LIBS)
 
 mostlyclean-compile:
 	-rm -f *.$(OBJEXT)
@@ -341,6 +359,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixPow_la-spowma.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixPow_la-zpowma.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testDoubleMatrixPow-testDoublePowm.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testFloatMatrixPow-testFloatPowm.Po@am__quote@
 
 .c.o:
 @am__fastdepCC_TRUE@	$(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@@ -405,6 +424,20 @@ testDoubleMatrixPow-testDoublePowm.obj: testDoublePowm.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testDoubleMatrixPow_CFLAGS) $(CFLAGS) -c -o testDoubleMatrixPow-testDoublePowm.obj `if test -f 'testDoublePowm.c'; then $(CYGPATH_W) 'testDoublePowm.c'; else $(CYGPATH_W) '$(srcdir)/testDoublePowm.c'; fi`
 
+testFloatMatrixPow-testFloatPowm.o: testFloatPowm.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testFloatMatrixPow_CFLAGS) $(CFLAGS) -MT testFloatMatrixPow-testFloatPowm.o -MD -MP -MF $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Tpo -c -o testFloatMatrixPow-testFloatPowm.o `test -f 'testFloatPowm.c' || echo '$(srcdir)/'`testFloatPowm.c
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Tpo $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testFloatPowm.c' object='testFloatMatrixPow-testFloatPowm.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testFloatMatrixPow_CFLAGS) $(CFLAGS) -c -o testFloatMatrixPow-testFloatPowm.o `test -f 'testFloatPowm.c' || echo '$(srcdir)/'`testFloatPowm.c
+
+testFloatMatrixPow-testFloatPowm.obj: testFloatPowm.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testFloatMatrixPow_CFLAGS) $(CFLAGS) -MT testFloatMatrixPow-testFloatPowm.obj -MD -MP -MF $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Tpo -c -o testFloatMatrixPow-testFloatPowm.obj `if test -f 'testFloatPowm.c'; then $(CYGPATH_W) 'testFloatPowm.c'; else $(CYGPATH_W) '$(srcdir)/testFloatPowm.c'; fi`
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Tpo $(DEPDIR)/testFloatMatrixPow-testFloatPowm.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testFloatPowm.c' object='testFloatMatrixPow-testFloatPowm.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testFloatMatrixPow_CFLAGS) $(CFLAGS) -c -o testFloatMatrixPow-testFloatPowm.obj `if test -f 'testFloatPowm.c'; then $(CYGPATH_W) 'testFloatPowm.c'; else $(CYGPATH_W) '$(srcdir)/testFloatPowm.c'; fi`
+
 mostlyclean-libtool:
 	-rm -f *.lo
 
diff --git a/src/matrixOperations/powm/cpowma.c b/src/matrixOperations/powm/cpowma.c
index 51095302..5fb8e365 100644
--- a/src/matrixOperations/powm/cpowma.c
+++ b/src/matrixOperations/powm/cpowma.c
@@ -11,15 +11,53 @@
  */
 
 #include "matrixPow.h"
-#include "logm.h"
-#include "matrixExponential.h"
-#include "multiplication.h"
+#include "spec.h"
+#include "pow.h"
+#include "matrixTranspose.h"
+#include "conj.h"
+#include "matrixInversion.h"
+#include "matrixMultiplication.h"
 
-void cpowma(floatComplex* in, int rows, floatComplex expand, floatComplex* out){
-	int i=0;
-	/* use the formula a^b=exp(b*ln(a)) */
-	clogma(in,rows,out);
-	for(i=0;i<rows*rows;i++) out[i]= cmuls(expand,out[i]);
-	cexpma(out,out,(int)rows);
-					
+void cpowma(floatComplex* in, int rows, floatComplex power, floatComplex* out){
+	int i=0, j=0;
+	int hermitian=0;
+	floatComplex *eigenvalues,*eigenvectors,*tmp;
+	
+	/* Data initialization */
+	eigenvalues = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	eigenvectors = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	tmp = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	
+	/* symmetric test*/
+	for(i=0;i<rows;i++) {
+		for (j=0;j<rows;j++)
+			if ((creals(in[i*rows+j])!=creals(in[j*rows+i])) || (cimags(in[i*rows+j])!=-cimags(in[j*rows+i]))) break;
+		
+		if (j!=rows) break;
+	}
+	
+	if ((i==rows)&&(j==rows)) hermitian=1;
+	
+	
+	/* find eigenvalues and eigenvectors */
+	cspec2a(in, rows, eigenvalues,eigenvectors);
+	
+	/* make operation on eigenvalues and eigenvectors */
+	for (i=0;i<rows;i++)
+			eigenvalues[i+i*rows]=cpows(eigenvalues[i+i*rows],power);
+			
+	cmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
+	
+	if (hermitian){		
+		ctransposea(eigenvectors, rows,rows, eigenvalues);
+		cconja(eigenvalues, rows*rows, eigenvalues);
+	}
+	else cinverma(eigenvectors, eigenvalues, rows);
+	
+	cmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	
+	free(eigenvalues);
+	free(eigenvectors);
+	free(tmp);
+	
 }
diff --git a/src/matrixOperations/powm/dpowma.c b/src/matrixOperations/powm/dpowma.c
index 076aa405..5d67189b 100644
--- a/src/matrixOperations/powm/dpowma.c
+++ b/src/matrixOperations/powm/dpowma.c
@@ -11,65 +11,54 @@
  */
 
 #include "matrixPow.h"
-#include "lapack.h"
-#include "eye.h"
+#include "spec.h"
+#include "pow.h"
+#include "matrixTranspose.h"
+#include "conj.h"
+#include "matrixInversion.h"
 #include "matrixMultiplication.h"
+#include <stdio.h>
 
-
-void dpowma(double* in, int size, double expand, double* out){
-
-#ifndef WITHOUT_BLAS		
-	switch ((int)expand){
-		case 0 : 
-			deyea(out,size,size);
-			break;
-		case 1 : 
-			{
-				int i;
-				for (i=0;i<size*size;i++) out[i]=in[i];
-			}
-			break;
-		default : 
-			{
-				int i=0,j=0,k=0;
-				int One=1;
-				
-				for (i=1; i<expand; i++){
-					for (j=0;j<size;j++) /*column*/ {
-						for (k=0;k<size;k++)/*row*/{
-							out[k+j*size]=ddot_(&size,in+k,&size,in+j*size,&One);
-						}
-					}
-				}
-			}
-			break;
+void dpowma(double* in, int rows, double power, doubleComplex* out){
+	int i=0, j=0;
+	int symmetric=0;
+	doubleComplex *eigenvalues,*eigenvectors,*tmp;
+	
+	/* Data initialization */
+	eigenvalues = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	eigenvectors = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	tmp = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	
+	/* symmetric test*/
+	for(i=0;i<rows;i++) {
+		for (j=0;j<rows;j++)
+			if (in[i*rows+j]!=in[j*rows+i]) break;
+		
+		if (j!=rows) break;
 	}
-#else
-	switch ((int)expand){
-		case 0 : 
-			deyea(out,size,size);
-			break;
-		case 1 : 
-			{
-				int i;
-				for (i=0;i<size*size;i++) out[i]=in[i];
-			}
-			break;
-		default : 
-			{
-				int i=0,j=0;	
-				double* Pow;
-				Pow=malloc((uint)(size*size)*sizeof(double));
-				for (i=0;i<size*size;i++) out[i]=in[i];
-				for (i=1; i<expand; i++){
-					for (j=0;j<size*size;j++) Pow[j]=out[j];
-					dmulma(Pow,size,size,in,size,size,out);
-				}
-			}
-			break;
+	
+	if ((i==rows)&&(j==rows)) symmetric=1;
+	
+	
+	/* find eigenvalues and eigenvectors */
+	dspec2a(in, 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++)
+			eigenvalues[i+i*rows]=zpows(eigenvalues[i+i*rows],DoubleComplex(power,0));
+			
+	zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
+	
+	if (symmetric){		
+		ztransposea(eigenvectors, rows,rows, eigenvalues);
+		zconja(eigenvalues, rows*rows, eigenvalues);
 	}
-					
-
-#endif
-					
+	else zinverma(eigenvectors, eigenvalues, rows);
+	
+	zmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	
+	free(eigenvalues);
+	free(eigenvectors);
+	free(tmp);
+	
 }
diff --git a/src/matrixOperations/powm/spowma.c b/src/matrixOperations/powm/spowma.c
index aa0f7043..634cd449 100644
--- a/src/matrixOperations/powm/spowma.c
+++ b/src/matrixOperations/powm/spowma.c
@@ -11,35 +11,53 @@
  */
 
 #include "matrixPow.h"
-#include "eye.h"
+#include "spec.h"
+#include "pow.h"
+#include "matrixTranspose.h"
+#include "conj.h"
+#include "matrixInversion.h"
 #include "matrixMultiplication.h"
 
-
-void spowma(float* in, int size, float expand, float* out){
-
-
-	switch ((int)expand){
-		case 0 : 
-			seyea(out,size,size);
-			break;
-		case 1 : 
-			{
-				int i;
-				for (i=0;i<size*size;i++) out[i]=in[i];
-			}
-			break;
-		default : 
-			{
-				int i=0,j=0;	
-				float* Pow;
-				Pow=malloc((uint)(size*size)*sizeof(float));
-				for (i=0;i<size*size;i++) out[i]=in[i];
-				for (i=1; i<expand; i++){
-					for (j=0;j<size*size;j++) Pow[j]=out[j];
-					smulma(Pow,size,size,in,size,size,out);
-				}
-			}
-			break;
+void spowma(float* in, int rows, float power, floatComplex* out){
+	int i=0, j=0;
+	int symmetric=0;
+	floatComplex *eigenvalues,*eigenvectors,*tmp;
+	
+	/* Data initialization */
+	eigenvalues = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	eigenvectors = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	tmp = malloc((uint)(rows*rows)*sizeof(floatComplex));
+	
+	/* symmetric test*/
+	for(i=0;i<rows;i++) {
+		for (j=0;j<rows;j++)
+			if (in[i*rows+j]!=in[j*rows+i]) break;
+		
+		if (j!=rows) break;
 	}
+	
+	if ((i==rows)&&(j==rows)) symmetric=1;
+	
+	
+	/* find eigenvalues and eigenvectors */
+	sspec2a(in, rows, eigenvalues,eigenvectors);
+	
+	/* make operation on eigenvalues and eigenvectors */
+	for (i=0;i<rows;i++)
+			eigenvalues[i+i*rows]=cpows(eigenvalues[i+i*rows],FloatComplex(power,0));
 			
+	cmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
+	
+	if (symmetric){		
+		ctransposea(eigenvectors, rows,rows, eigenvalues);
+		cconja(eigenvalues, rows*rows, eigenvalues);
+	}
+	else cinverma(eigenvectors, eigenvalues, rows);
+	
+	cmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	
+	free(eigenvalues);
+	free(eigenvectors);
+	free(tmp);
+	
 }
diff --git a/src/matrixOperations/powm/testDoublePowm.c b/src/matrixOperations/powm/testDoublePowm.c
index 6fdbec10..05cad5e5 100644
--- a/src/matrixOperations/powm/testDoublePowm.c
+++ b/src/matrixOperations/powm/testDoublePowm.c
@@ -17,38 +17,126 @@
 
 
 static void dpowmaTest(void){
-	double in[4]={1,5,4,2};
-	double out[4];
+	double in1[4]={1,5,4,2};
+	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];
 	int i;
 	
-	dpowma(in, 2, 2, out);
+	double in2[16]={    2.5358983855694532394409 ,   9.0725262500345706939697,    0.0026536155492067337036,    3.9639251008629798889160 , 
+	7.9845732506364583969116,    7.5407014600932598114014,   10.196942830458283424377 ,   8.2287722378969192504883 , 
+	10.538597775623202323914,    0.8204884417355060577393,    6.7301832754164934158325,    7.9482832476496696472168,  
+	8.7162081208080053329468 ,   2.3821726106107234954834 ,   6.5310877952724695205688,    2.784897476434707641602   };
+	double expand2 = 3.4683557949028909206390;
+	double result2R[16]={13801.893971410685480805 ,   9622.6108799100766191259 ,   10325.586569611912636901,   10694.791005280343597406  ,
+	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];
 	
-	for (i=0;i<4;i++) printf("out[%d] = %f\n",i,out[i]);
+	
+	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);
+	}
+	
+/*	
+	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);
+	}	
 	
 }
 
-/* FIXME : assert à10^-12 */
+/* FIXME : assert 1e-14 */
 static void zpowmaTest(void){
-	double inR[9]={1,2,3,4,5,6,7,8,9};
-	double inI[9]={1,2,3,4,5,6,7,8,9};
-	double resultR[9]={- 4.7115011361608578610571,- 2.0782061409646632732517,0.5550888542315330909105,
-	- 2.3202132490900626571317,- 2.4412168031527574640904,- 2.5622203572154611528333,
-	0.0710746379807356554181,- 2.80422746534086453352,- 5.6795295686624518438634};
-	double resultI[9]={- 12.188702380084603049681,- 4.0827818504168584823333,4.0231386792508754268738,
-	- 3.0919079733956360556135,- 2.5964710348850239540752,- 2.1010340963744131848046,
-	6.0048864332933264975622,- 1.1101602193531934226201,- 8.2252068719997026846613};
-	doubleComplex *in,out[9];
-	int i;
+	/* Tests 1 */
+	{
+		double inR[9]={1,2,3,4,5,6,7,8,9};
+		double inI[9]={1,2,3,4,5,6,7,8,9};
+		double resultR[9]={- 4.7115011361608578610571,- 2.0782061409646632732517,0.5550888542315330909105,
+		- 2.3202132490900626571317,- 2.4412168031527574640904,- 2.5622203572154611528333,
+		0.0710746379807356554181,- 2.80422746534086453352,- 5.6795295686624518438634};
+		double resultI[9]={- 12.188702380084603049681,- 4.0827818504168584823333,4.0231386792508754268738,
+		- 3.0919079733956360556135,- 2.5964710348850239540752,- 2.1010340963744131848046,
+		6.0048864332933264975622,- 1.1101602193531934226201,- 8.2252068719997026846613};
+		doubleComplex *in,out[9];
+		int i;
+	
+		in=DoubleComplexMatrix(inR,inI,9);
+	
+		zpowma(in, 3, DoubleComplex(1,1), out);
+	
+		for (i=0;i<9;i++) printf("out[%d] = %f+%f*i\n",i,zreals(out[i]),zimags(out[i]));
+	
+		for (i=0;i<9;i++){
+			assert( (fabs(zreals(out[i])-resultR[i])/ fabs(zreals(out[i])) ) <1e-14);
+			assert( (fabs(zimags(out[i])-resultI[i])/ fabs(zimags(out[i])) ) <1e-14);
+		}
+	}
+	
+	/* Tests 2 and 3 */
+	{
+		double in1R[4]={1,5,4,2};
+		double in1I[4]={0};
+		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];
+		int i;
+	
+		double in2R[16]={    2.5358983855694532394409 ,   9.0725262500345706939697,    0.0026536155492067337036,    3.9639251008629798889160 , 
+		7.9845732506364583969116,    7.5407014600932598114014,   10.196942830458283424377 ,   8.2287722378969192504883 , 
+		10.538597775623202323914,    0.8204884417355060577393,    6.7301832754164934158325,    7.9482832476496696472168,  
+		8.7162081208080053329468 ,   2.3821726106107234954834 ,   6.5310877952724695205688,    2.784897476434707641602   };
+		double in2I[16]={0};
+		double expand2 = 3.4683557949028909206390;
+		double result2R[16]={13801.893971410685480805 ,   9622.6108799100766191259 ,   10325.586569611912636901,   10694.791005280343597406  ,
+		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];
+		doubleComplex *in1,*in2;
 	
-	in=DoubleComplexMatrix(inR,inI,9);
 	
-	zpowma(in, 3, DoubleComplex(1,1), out);
+		in1=DoubleComplexMatrix(in1R,in1I,4);
+		in2=DoubleComplexMatrix(in2R,in2I,16);
+		
+		zpowma(in1, 2, DoubleComplex(expand1,0), out1);
+		zpowma(in2, 4, DoubleComplex(expand2,0), out2);
 	
-	for (i=0;i<9;i++) printf("out[%d] = %f+%f*i\n",i,zreals(out[i]),zimags(out[i]));
+		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);
+		}
 	
-	for (i=0;i<9;i++){
-		assert( (fabs(zreals(out[i])-resultR[i])/ fabs(zreals(out[i])) ) <1e-12);
-		assert( (fabs(zimags(out[i])-resultI[i])/ fabs(zimags(out[i])) ) <1e-13);
+	/*	
+		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);
+		}	
 	}
 	
 }
diff --git a/src/matrixOperations/powm/zpowma.c b/src/matrixOperations/powm/zpowma.c
index c8cac267..2def5f40 100644
--- a/src/matrixOperations/powm/zpowma.c
+++ b/src/matrixOperations/powm/zpowma.c
@@ -11,15 +11,53 @@
  */
 
 #include "matrixPow.h"
-#include "logm.h"
-#include "matrixExponential.h"
-#include "multiplication.h"
+#include "spec.h"
+#include "pow.h"
+#include "matrixTranspose.h"
+#include "conj.h"
+#include "matrixInversion.h"
+#include "matrixMultiplication.h"
 
-void zpowma(doubleComplex* in, int rows, doubleComplex expand, doubleComplex* out){
-	int i=0;
-	/* use the formula a^b=exp(b*ln(a)) */
-	zlogma(in,rows,out);
-	for(i=0;i<rows*rows;i++) out[i]= zmuls(expand,out[i]);
-	zexpma(out,out,(int)rows);
-					
+void zpowma(doubleComplex* in, int rows, doubleComplex power, doubleComplex* out){
+	int i=0, j=0;
+	int hermitian=0;
+	doubleComplex *eigenvalues,*eigenvectors,*tmp;
+	
+	/* Data initialization */
+	eigenvalues = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	eigenvectors = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	tmp = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+	
+	/* hermitian test*/
+	for(i=0;i<rows;i++) {
+		for (j=0;j<rows;j++)
+			if ((zreals(in[i*rows+j])!=zreals(in[j*rows+i])) || (zimags(in[i*rows+j])!=-zimags(in[j*rows+i]))) break;
+		
+		if (j!=rows) break;
+	}
+	
+	if ((i==rows)&&(j==rows)) hermitian=1;
+	
+	
+	/* find eigenvalues and eigenvectors */
+	zspec2a(in, rows, eigenvalues,eigenvectors);
+	
+	/* make operation on eigenvalues and eigenvectors */
+	for (i=0;i<rows;i++)
+			eigenvalues[i+i*rows]=zpows(eigenvalues[i+i*rows],power);
+			
+	zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
+	
+	if (hermitian){		
+		ztransposea(eigenvectors, rows,rows, eigenvalues);
+		zconja(eigenvalues, rows*rows, eigenvalues);
+	}
+	else zinverma(eigenvectors, eigenvalues, rows);
+	
+	zmulma(tmp, rows, rows, eigenvalues, rows, rows, out);
+	
+	free(eigenvalues);
+	free(eigenvectors);
+	free(tmp);
+	
 }