summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/matrixOperations/logm/Makefile.am1
-rw-r--r--src/matrixOperations/logm/Makefile.in1
-rw-r--r--src/matrixOperations/logm/clogma.c12
-rw-r--r--src/matrixOperations/logm/dlogma.c8
-rw-r--r--src/matrixOperations/logm/slogma.c12
-rw-r--r--src/matrixOperations/logm/testDoubleLogm.c79
-rw-r--r--src/matrixOperations/logm/zlogma.c83
7 files changed, 144 insertions, 52 deletions
diff --git a/src/matrixOperations/logm/Makefile.am b/src/matrixOperations/logm/Makefile.am
index df939f92..85e8d7f2 100644
--- a/src/matrixOperations/logm/Makefile.am
+++ b/src/matrixOperations/logm/Makefile.am
@@ -50,6 +50,7 @@ check_LDADD = $(top_builddir)/type/libDoubleComplex.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 \
diff --git a/src/matrixOperations/logm/Makefile.in b/src/matrixOperations/logm/Makefile.in
index e662a1a0..93b2dc6f 100644
--- a/src/matrixOperations/logm/Makefile.in
+++ b/src/matrixOperations/logm/Makefile.in
@@ -218,6 +218,7 @@ check_LDADD = $(top_builddir)/type/libDoubleComplex.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 \
diff --git a/src/matrixOperations/logm/clogma.c b/src/matrixOperations/logm/clogma.c
index f2614b0e..c4e69d44 100644
--- a/src/matrixOperations/logm/clogma.c
+++ b/src/matrixOperations/logm/clogma.c
@@ -15,20 +15,20 @@
#include "logm.h"
-void clogma (floatComplex* in, int size, floatComplex* out){
+void clogma (floatComplex* in, int rows, floatComplex* out){
doubleComplex *inCopy, *outCopy;
int i=0;
- inCopy = malloc((uint)(size*size)*sizeof(doubleComplex));
- outCopy = malloc((uint)(size*size)*sizeof(doubleComplex));
+ inCopy = malloc((uint)(rows*rows)*sizeof(doubleComplex));
+ outCopy = malloc((uint)(rows*rows)*sizeof(doubleComplex));
- for(i=0;i<size*size;i++)
+ for(i=0;i<rows*rows;i++)
inCopy[i]=DoubleComplex ((double)creals(in[i]), (double)cimags(in[i]));
- zlogma(inCopy,size,outCopy);
+ zlogma(inCopy,rows,outCopy);
- for(i=0;i<size*size;i++)
+ for(i=0;i<rows*rows;i++)
out[i]=FloatComplex( (float)zreals(outCopy[i]),(float)zimags(outCopy[i]));
diff --git a/src/matrixOperations/logm/dlogma.c b/src/matrixOperations/logm/dlogma.c
index 8a4986b5..2798e127 100644
--- a/src/matrixOperations/logm/dlogma.c
+++ b/src/matrixOperations/logm/dlogma.c
@@ -14,15 +14,15 @@
#include "logm.h"
-void dlogma (double* in, int size, doubleComplex* out){
+void dlogma (double* in, int rows, doubleComplex* out){
doubleComplex *inCpx;
int i;
- inCpx=malloc((uint)(size*size)*sizeof(doubleComplex));
+ inCpx=malloc((uint)(rows*rows)*sizeof(doubleComplex));
- for (i=0;i<size*size;i++) inCpx[i] = DoubleComplex(in[i],0);
+ for (i=0;i<rows*rows;i++) inCpx[i] = DoubleComplex(in[i],0);
- zlogma(inCpx, size, out);
+ zlogma(inCpx, rows, out);
free(inCpx);
}
diff --git a/src/matrixOperations/logm/slogma.c b/src/matrixOperations/logm/slogma.c
index 930bf752..a42929f0 100644
--- a/src/matrixOperations/logm/slogma.c
+++ b/src/matrixOperations/logm/slogma.c
@@ -14,18 +14,18 @@
#include "logm.h"
-void slogma (float* in, int size, floatComplex* out){
+void slogma (float* in, int rows, floatComplex* out){
doubleComplex *inCpx, *outCopy;
int i;
- inCpx=malloc((uint)(size*size)*sizeof(doubleComplex));
- outCopy=malloc((uint)(size*size)*sizeof(doubleComplex));
+ inCpx=malloc((uint)(rows*rows)*sizeof(doubleComplex));
+ outCopy=malloc((uint)(rows*rows)*sizeof(doubleComplex));
- for (i=0;i<size*size;i++) inCpx[i] = DoubleComplex(in[i],0);
+ for (i=0;i<rows*rows;i++) inCpx[i] = DoubleComplex(in[i],0);
- zlogma(inCpx, size, outCopy);
+ zlogma(inCpx, rows, outCopy);
- for(i=0;i<size*size;i++)
+ for(i=0;i<rows*rows;i++)
out[i]=FloatComplex( (float)zreals(outCopy[i]),(float)zimags(outCopy[i]));
free(inCpx);
diff --git a/src/matrixOperations/logm/testDoubleLogm.c b/src/matrixOperations/logm/testDoubleLogm.c
index 8b31b857..f54b2af9 100644
--- a/src/matrixOperations/logm/testDoubleLogm.c
+++ b/src/matrixOperations/logm/testDoubleLogm.c
@@ -16,16 +16,88 @@
#include <stdio.h>
#include "logm.h"
+
+/* FIXME : Untested*/
+
static void dlogmaTest(void){
-}
+ int i=0;
+
+ double in4[4] = {0.0683740368112921714783,0.5608486062847077846527,
+ 0.6623569373041391372681,0.7263506767340004444122 };
+ double in9[9] = {0.1985143842175602912903,0.5442573162727057933807,0.2320747897028923034668,
+ 0.2312237196601927280426,0.2164632631465792655945,0.8833887814544141292572,
+ 0.6525134947150945663452,0.3076090742833912372589,0.9329616213217377662659};
+ double in16[16] = {2.1460078610107302665710,3.126419968903064727783,3.6163610080257058143616,2.922266637906432151794,
+ 5.664248815737664699554,4.8264719732105731964111,3.3217189135029911994934,5.9350947011262178421021,
+ 5.0153415976092219352722,4.3685875833034515380859,2.6931248093023896217346,6.3257448654621839523315,
+ 4.051954015158116817474,9.1847078315913677215576,0.4373343335464596748352,4.818508932366967201233};
+ double inHer[16] = {2,3,4,5,3,1,7,9,4,7,2,4,5,9,4,7};
+
+ double result4R[4]={- 0.8770746698483730119378,0.5288031941140065583795,
+ 0.6245116064569415925689,- 0.2566930821175439358939};
+ double result4I[4]={2.3169161534469857599561,- 1.2719608554746495432397,
+ - 1.5021738258100643115256,0.8246765001428078001311};
+ double result9R[9]={- 0.8537183024765053751537,1.9622873186332121520081,- 0.8931518563294031931577,
+ - 0.6983367178614473536,- 1.1199049557115796638129,1.539392415892589571769,
+ 1.0815849580686147657,- 0.4078285487954845245362,0.1789837218054647349774,};
+ double result9I[9]={- 0.0000000000000006661338,0.0000000000000005551115,0.0000000000000001110223,
+ 0.0000000000000002844947,0.0000000000000004163336,- 0.0000000000000008743006,
+ 0.0000000000000001665335,- 0.0000000000000004996004,0.0000000000000003330669};
+ double result16R[16]={1.4069394184534917968676,- 1.7915236034592834091228,2.3099248310776610004780,0.5320777350783232328979,
+ 0.6986579270151970222358,0.9874373592109859654187,0.8902400742596190408307,0.6067035092100719895214,
+ 0.0466142343271527676007,0.3821963653954382200695,1.129789301050062144682,1.2310339071369393870725,
+ 0.424502387440284123521,3.4580528009566289249221,- 2.1337657344549150195689,1.0826836628584670663145};
+ double result16I[16]={1.4875734535892033427729,- 1.8146634352088595321106,0.0551064758054355952055,0.8617364610492967980093,
+ - 1.0892397305420760122985,1.3287434690517783142383,- 0.0403503858667347281575,- 0.6309857092212021179023,
+ - 1.3838992102024216813,1.6881931367553519862668,- 0.0512659110447117594145,- 0.8016790061501579689463,
+ 0.6500053973791579675634,- 0.7929296025459252605927,0.0240791515993044047406,0.3765416419935233571792};
+ double resultHerR[16]={- 0.0707459665791660696765,0.9744491331561414559914,0.0316695005563280007621,0.9186825501429138896015,
+ 0.9744491331561416780360,1.9318283614573175110962,- 0.0664515082056250649956,0.3546926741474874522631,
+ 0.0316695005563280562733,- 0.0664515082056250649956,0.3206743709780528472919,1.82590946141052867802,
+ 0.9186825501429141116461,0.3546926741474875632854,1.82590946141052867802,0.9537374500729456361370};
+ double resultHerI[16]={2.1383917599789858954296,- 0.5351429414718615884539,- 1.3567950682942089279948,0.1340361867044403687554,
+ - 0.5351429414718615884539,2.2766753224580762449136,- 0.6151069729626167381653,- 1.1421039303668867326280,
+ - 1.3567950682942089279948,- 0.6151069729626165161207,1.2861990290089324595613,0.4088471998565592624431,
+ 0.1340361867044403965110,- 1.1421039303668867326280,0.4088471998565592624431,0.5819191957335942966267};
+
+ doubleComplex out4[4],out9[9],out16[16],outHer[16];
+
+
+ dlogma(in4,2,out4);
+ dlogma(in9,3,out9);
+ dlogma(in16,4,out16);
+ dlogma(inHer,4,outHer);
+
+ for(i=0;i<4;i++) {
+ assert( fabs(zreals(out4[i])-result4R[i]) / fabs(zreals(out4[i])) <1e-15);
+ assert( fabs(zimags(out4[i])-result4I[i]) / fabs(zimags(out4[i])) <1e-15);
+ }
+
+ for(i=0;i<9;i++) {
+ assert( fabs(zreals(out9[i])-result9R[i]) / fabs(zreals(out9[i])) <1e-14);
+ if (zimags(out9[i])>1e-15) assert( fabs(zimags(out9[i])-result9I[i]) / fabs(zimags(out9[i])) <1e-16);
+ else assert(1);
+ }
+
+ for(i=0;i<16;i++) {
+ printf("out[%d] = %f + %f *i ---result = %f + %f *i\n",i,zreals(out16[i]),zimags(out16[i]),result16R[i],result16I[i]);
+ assert( fabs(zreals(out16[i])-result16R[i]) / fabs(zreals(out16[i])) <1e-13);
+ /*assert( fabs(zimags(out16[i])-result16I[i]) / fabs(zimags(out16[i])) <1e-16);*/
+ }
+
+ for(i=0;i<16;i++) {
+ printf("out[%d] = %f + %f *i ---result = %f + %f *i\n",i,zreals(outHer[i]),zimags(outHer[i]),resultHerR[i],resultHerI[i]);
+ assert( fabs(zreals(outHer[i])-resultHerR[i]) / fabs(zreals(outHer[i])) <1e-12);
+ assert( fabs(zimags(outHer[i])-resultHerI[i]) / fabs(zimags(outHer[i])) <1e-13);
+ }
+}
static void zlogmaTest(void){
int i;
double inD[9]={0.2113249,0.3303271 , 0.8497452 , 0.7560439, 0.6653811 , 0.6857310 ,0.0002211 , 0.6283918 , 0.8782165};
- doubleComplex *in, *out;
+ doubleComplex *in, out[9];
in=malloc((uint)9*sizeof(doubleComplex));
- out=malloc((uint)9*sizeof(doubleComplex));
for(i=0;i<9;i++) in[i]=DoubleComplex(inD[i],0);
@@ -33,7 +105,6 @@ static void zlogmaTest(void){
for(i=0;i<9;i++) printf("out[%d] = %f + %f *i\n",i,zreals(out[i]),zimags(out[i]));
free(in);
- free(out);
}
static int logmTest(void){
diff --git a/src/matrixOperations/logm/zlogma.c b/src/matrixOperations/logm/zlogma.c
index 37516822..5a10151e 100644
--- a/src/matrixOperations/logm/zlogma.c
+++ b/src/matrixOperations/logm/zlogma.c
@@ -19,8 +19,9 @@
#include "matrixTranspose.h"
#include "matrixInversion.h"
#include "max.h"
+#include "conj.h"
-void zlogma (doubleComplex* in, int size, doubleComplex* out){
+void zlogma (doubleComplex* in, int rows, doubleComplex* out){
/* Algo : */
/* trouver les valeurs propres vp */
@@ -31,54 +32,67 @@ void zlogma (doubleComplex* in, int size, doubleComplex* out){
int i;
doubleComplex *eigenvalues, *eigenvectors, *tmp;
- int symetrique = 0;
+ 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 */
- doubleComplex * inCopy = NULL; /* Used by LAPACK */
/* Data initialization */
- eigenvalues = malloc((uint)(size*size)*sizeof(doubleComplex));
- eigenvectors = malloc((uint)(size*size)*sizeof(doubleComplex));
- tmp = malloc((uint)(size*size)*sizeof(doubleComplex));
- iWorkSize = max(1,2*size);
+ 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*size)*sizeof(doubleComplex));
- inCopy = malloc((uint)(size*size)*sizeof(doubleComplex));
+ pdblRWork = malloc((uint)(2*rows)*sizeof(doubleComplex));
+
/* Copy in in inCopy */
- for (i=0;i<size*size;i++)
- inCopy[i] = in[i];
+ for (i=0;i<rows*rows;i++)
+ out[i] = in[i];
- /* regarde si in est symetrique */
- ztransposea(in,size,size,tmp);
- for (i=0;i<size*size;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;
}
- if (i==size*size) symetrique = 1;
+ if (i==rows*rows) hermitienne = 1;
/* trouver les valeurs propres vp ainsi que les vecteurs propres*/
- if (symetrique){
- C2F(zheev)( "V", "U", &size, inCopy, &size, eigenvalues,
+ if (hermitienne){
+ printf("wazaaaaaaa\n");
+ C2F(zheev)( "V", "U", &rows, out, &rows, eigenvalues,
pdblWork, &iWorkSize, pdblRWork, &info );
- eigenvectors = inCopy;
+ /* 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 {
- C2F(zgeev)( "N", "V", &size, inCopy, &size, eigenvalues,
- pdblLeftvectors, &size, eigenvectors, &size, pdblWork, &iWorkSize,
+ printf("passe\n");
+ C2F(zgeev)( "N", "V", &rows, out, &rows, eigenvalues,
+ pdblLeftvectors, &rows, eigenvectors, &rows, pdblWork, &iWorkSize,
pdblRWork, &info );
- for(i=0;i<size;i++) {
- eigenvalues[i*size+i] = eigenvalues[i];
- if ((i*size+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
+ for(i=0;i<rows;i++) {
+ eigenvalues[i*rows+i] = eigenvalues[i];
+ if ((i*rows+i)!=i) eigenvalues[i] = DoubleComplex(0,0);
}
}
@@ -88,28 +102,33 @@ void zlogma (doubleComplex* in, int size, doubleComplex* out){
/* diag(log(diag(vp)) */
- for (i=0;i<size*size;i++){
- if ((i%(size+1))==0) /* teste si i est sur la diagonale */
+ for (i=0;i<rows*rows;i++){
+ if ((i%(rows+1))==0) /* teste si i est sur la diagonale */
eigenvalues[i] = zlogs(eigenvalues[i]);
else eigenvalues[i] = DoubleComplex(0,0);
}
/* Vp * diag(log(diag(vp)) */
- zmulma(eigenvectors, size, size, eigenvalues, size, size, tmp);
+ zmulma(eigenvectors, rows, rows, eigenvalues, rows, rows, tmp);
- /* inv(Vp) */
- zinverma(eigenvectors, eigenvectors, size);
+ /* Vp' ou inv(Vp) */
+ if (hermitienne) {
+ /* we use eigenvalues as a temporary matrix cause it's useless now*/
+ ztransposea(eigenvectors,rows,rows,eigenvalues);
+ zconja(eigenvalues,rows*rows,eigenvectors);
+ }
+ else zinverma(eigenvectors, eigenvectors, rows);
/* Vp * diag(log(diag(vp))*inv(Vp) */
- zmulma(tmp, size, size, eigenvectors, size, size, out);
+ zmulma(tmp, rows, rows, eigenvectors, rows, rows, out);
+ /* FIXME : pb with eigenvectors */
+ eigenvectors=NULL;
free(eigenvalues);
- free(eigenvectors);
+ free(eigenvectors);
free(tmp);
free(pdblWork);
free(pdblRWork);
- free(inCopy);
-
}
generated by cgit (git 2.34.1) at 2025-03-18 02:37:16 +0000