summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/matrixOperations/logm/Makefile.am39
-rw-r--r--src/matrixOperations/logm/Makefile.in39
-rw-r--r--src/matrixOperations/logm/zlogma.c69
3 files changed, 54 insertions, 93 deletions
diff --git a/src/matrixOperations/logm/Makefile.am b/src/matrixOperations/logm/Makefile.am
index 85e8d7f2..cb69e8f0 100644
--- a/src/matrixOperations/logm/Makefile.am
+++ b/src/matrixOperations/logm/Makefile.am
@@ -42,23 +42,25 @@ libLogm_la_SOURCES = $(HEAD) \
check_PROGRAMS = testDoubleLogm
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
- $(top_builddir)/type/libFloatComplex.la \
- $(top_builddir)/lib/lapack/libscilapack.la \
- $(top_builddir)/lib/blas/libsciblas.la \
- $(top_builddir)/elementaryFunctions/log/libLog.la \
- $(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
- $(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
- $(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
- $(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
- $(top_builddir)/auxiliaryFunctions/conj/libConj.la \
- $(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
- $(top_builddir)/auxiliaryFunctions/sign/libSign.la \
- $(top_builddir)/operations/multiplication/libMultiplication.la \
- $(top_builddir)/operations/addition/libAddition.la \
- $(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
- $(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
- $(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
- $(top_builddir)/matrixOperations/logm/libLogm.la
+ $(top_builddir)/type/libFloatComplex.la \
+ $(top_builddir)/lib/lapack/libscilapack.la \
+ $(top_builddir)/lib/blas/libsciblas.la \
+ $(top_builddir)/elementaryFunctions/log/libLog.la \
+ $(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+ $(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
+ $(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
+ $(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+ $(top_builddir)/auxiliaryFunctions/conj/libConj.la \
+ $(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
+ $(top_builddir)/auxiliaryFunctions/sign/libSign.la \
+ $(top_builddir)/operations/multiplication/libMultiplication.la \
+ $(top_builddir)/operations/addition/libAddition.la \
+ $(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
+ $(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
+ $(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
+ $(top_builddir)/matrixOperations/spec2/libSpec2.la \
+ $(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
+ $(top_builddir)/matrixOperations/logm/libLogm.la
check_INCLUDES = -I ./includes \
@@ -66,8 +68,7 @@ check_INCLUDES = -I ./includes \
-I $(top_builddir)/elementaryFunctions/includes \
-I $(top_builddir)/operations/includes \
-I $(top_builddir)/auxiliaryFunctions/includes \
- -I $(top_builddir)/matrixOperations/includes\
- -I $(top_builddir)/matrixOperations/logm
+ -I $(top_builddir)/matrixOperations/includes
testDoubleLogm_SOURCES = testDoubleLogm.c
testDoubleLogm_LDADD = $(check_LDADD)
diff --git a/src/matrixOperations/logm/Makefile.in b/src/matrixOperations/logm/Makefile.in
index 93b2dc6f..cb74434a 100644
--- a/src/matrixOperations/logm/Makefile.in
+++ b/src/matrixOperations/logm/Makefile.in
@@ -210,31 +210,32 @@ libLogm_la_SOURCES = $(HEAD) \
dlogma.c
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
- $(top_builddir)/type/libFloatComplex.la \
- $(top_builddir)/lib/lapack/libscilapack.la \
- $(top_builddir)/lib/blas/libsciblas.la \
- $(top_builddir)/elementaryFunctions/log/libLog.la \
- $(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
- $(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
- $(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
- $(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
- $(top_builddir)/auxiliaryFunctions/conj/libConj.la \
- $(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
- $(top_builddir)/auxiliaryFunctions/sign/libSign.la \
- $(top_builddir)/operations/multiplication/libMultiplication.la \
- $(top_builddir)/operations/addition/libAddition.la \
- $(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
- $(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
- $(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
- $(top_builddir)/matrixOperations/logm/libLogm.la
+ $(top_builddir)/type/libFloatComplex.la \
+ $(top_builddir)/lib/lapack/libscilapack.la \
+ $(top_builddir)/lib/blas/libsciblas.la \
+ $(top_builddir)/elementaryFunctions/log/libLog.la \
+ $(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+ $(top_builddir)/elementaryFunctions/log1p/libLog1p.la \
+ $(top_builddir)/elementaryFunctions/lnp1m1/libLnp1m1.la \
+ $(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+ $(top_builddir)/auxiliaryFunctions/conj/libConj.la \
+ $(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
+ $(top_builddir)/auxiliaryFunctions/sign/libSign.la \
+ $(top_builddir)/operations/multiplication/libMultiplication.la \
+ $(top_builddir)/operations/addition/libAddition.la \
+ $(top_builddir)/matrixOperations/multiplication/libMatrixMultiplication.la \
+ $(top_builddir)/matrixOperations/transpose/libMatrixTranspose.la \
+ $(top_builddir)/matrixOperations/inversion/libMatrixInversion.la \
+ $(top_builddir)/matrixOperations/spec2/libSpec2.la \
+ $(top_builddir)/matrixOperations/zeros/libMatrixZeros.la \
+ $(top_builddir)/matrixOperations/logm/libLogm.la
check_INCLUDES = -I ./includes \
-I $(top_builddir)/type \
-I $(top_builddir)/elementaryFunctions/includes \
-I $(top_builddir)/operations/includes \
-I $(top_builddir)/auxiliaryFunctions/includes \
- -I $(top_builddir)/matrixOperations/includes\
- -I $(top_builddir)/matrixOperations/logm
+ -I $(top_builddir)/matrixOperations/includes
testDoubleLogm_SOURCES = testDoubleLogm.c
testDoubleLogm_LDADD = $(check_LDADD)
diff --git a/src/matrixOperations/logm/zlogma.c b/src/matrixOperations/logm/zlogma.c
index 5a10151e..3a629bfd 100644
--- a/src/matrixOperations/logm/zlogma.c
+++ b/src/matrixOperations/logm/zlogma.c
@@ -20,6 +20,7 @@
#include "matrixInversion.h"
#include "max.h"
#include "conj.h"
+#include "spec.h"
void zlogma (doubleComplex* in, int rows, doubleComplex* out){
@@ -30,77 +31,37 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
* logm = Vp * diag(log(diag(vp)) * inv(Vp) */
- int i;
+ int i = 0,j = 0;
doubleComplex *eigenvalues, *eigenvectors, *tmp;
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 */
-
-
-
+
/* Data initialization */
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*rows)*sizeof(doubleComplex));
-
- /* Copy in in inCopy */
- for (i=0;i<rows*rows;i++)
- out[i] = in[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;
+ for (i=0;i<rows;i++){
+ for(j=0;j<rows;j++)
+ if ( (zreals(in[i+j*rows])!=zreals(in[j+i*rows])) ||
+ (zimags(in[i+j*rows])!=zimags(in[j+i*rows])) ) break;
+ if (j!=rows) break;
}
- if (i==rows*rows) hermitienne = 1;
+ if ((i==rows) && (j==rows)) hermitienne=1;
/* trouver les valeurs propres vp ainsi que les vecteurs propres*/
- if (hermitienne){
- printf("wazaaaaaaa\n");
- C2F(zheev)( "V", "U", &rows, out, &rows, eigenvalues,
- pdblWork, &iWorkSize, pdblRWork, &info );
- /* 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 {
- printf("passe\n");
- C2F(zgeev)( "N", "V", &rows, out, &rows, eigenvalues,
- pdblLeftvectors, &rows, eigenvectors, &rows, pdblWork, &iWorkSize,
- pdblRWork, &info );
- for(i=0;i<rows;i++) {
- eigenvalues[i*rows+i] = eigenvalues[i];
- if ((i*rows+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
- }
- }
+ zspec2a(in,rows,eigenvalues,eigenvectors);
/* utiliser la formule suivante
* logm = Vp * diag(log(diag(vp)) * inv(Vp) */
-
+
/* diag(log(diag(vp)) */
for (i=0;i<rows*rows;i++){
if ((i%(rows+1))==0) /* teste si i est sur la diagonale */
@@ -108,6 +69,8 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
else eigenvalues[i] = DoubleComplex(0,0);
}
+
+
/* Vp * diag(log(diag(vp)) */
zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
@@ -122,13 +85,9 @@ void zlogma (doubleComplex* in, int rows, doubleComplex* out){
/* Vp * diag(log(diag(vp))*inv(Vp) */
zmulma(tmp, rows, rows, eigenvectors, rows, rows, out);
-
- /* FIXME : pb with eigenvectors */
- eigenvectors=NULL;
+
free(eigenvalues);
free(eigenvectors);
free(tmp);
- free(pdblWork);
- free(pdblRWork);
}
generated by cgit (git 2.34.1) at 2025-03-13 16:58:47 +0000