From b074501da9007198d67b37dd3527eeef2bb51dc8 Mon Sep 17 00:00:00 2001 From: simon Date: Mon, 11 Aug 2008 12:32:00 +0000 Subject: all test passed now but problem of precision --- src/matrixOperations/expm/cexpma.c | 22 ++++----- src/matrixOperations/expm/dexpma.c | 4 +- src/matrixOperations/expm/sexpma.c | 12 +++-- src/matrixOperations/expm/testMatrixExponential.c | 54 +++++++++++++++++++++-- src/matrixOperations/expm/zexpma.c | 14 ++++-- 5 files changed, 83 insertions(+), 23 deletions(-) (limited to 'src/matrixOperations') diff --git a/src/matrixOperations/expm/cexpma.c b/src/matrixOperations/expm/cexpma.c index 97d24b42..4bb4d836 100644 --- a/src/matrixOperations/expm/cexpma.c +++ b/src/matrixOperations/expm/cexpma.c @@ -24,9 +24,9 @@ void cexpma(floatComplex * in, floatComplex * out, int _iLeadDim) float fltExp = 0; float fltS = 0; float fltS2 = 0; - float fltCst = 0.5f; + double fltCst = 0.5; - floatComplex cfltCst ; + floatComplex *pfltMatrixA = NULL; floatComplex *pfltMatrixX = NULL; @@ -51,7 +51,6 @@ void cexpma(floatComplex * in, floatComplex * out, int _iLeadDim) pfltMatrixTemp = (floatComplex*)malloc(sizeof(floatComplex) * (unsigned int) iSquare); pfltMatrixTemp2 = (floatComplex*)malloc(sizeof(floatComplex) * (unsigned int) iSquare); - cfltCst = FloatComplex ( 0.5f , 0 ); /*// Scale A by power of 2 so that its norm is < 1/2 .*/ sfrexps( cinfnorma( in, _iLeadDim, _iLeadDim) , &fltExp); @@ -73,12 +72,13 @@ void cexpma(floatComplex * in, floatComplex * out, int _iLeadDim) ceyesa(pfltMatrixEye, _iLeadDim, _iLeadDim); - cmulma ( & cfltCst , 1 ,1, + /*cmulma ( & cfltCst , 1 ,1, pfltMatrixA , _iLeadDim, _iLeadDim, - pfltMatrixcA); + pfltMatrixcA);*/ /* cA = A * c */ - + for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) + pfltMatrixcA[iIndex1] = ctimess ( pfltMatrixA[iIndex1] , FloatComplex((float) fltCst , 0) ) ; /*E = Eye + cA*/ @@ -93,8 +93,7 @@ void cexpma(floatComplex * in, floatComplex * out, int _iLeadDim) for(iLoop1 = 2 ; iLoop1 <= iMax ; iLoop1++) { - fltCst *= (float) ((iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1))); - cfltCst = FloatComplex( fltCst , 0); + fltCst = fltCst * (iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1)); /*Temp = X */ /*C2F(zcopy)(&iSquare, pfltMatrixX, &iOne, pfltMatrixTemp, &iOne);*/ @@ -107,9 +106,12 @@ void cexpma(floatComplex * in, floatComplex * out, int _iLeadDim) pfltMatrixX ); /* cX = c * X */ - cmulma ( & cfltCst , 1 ,1, + /* cmulma ( & cfltCst , 1 ,1, pfltMatrixX , _iLeadDim, _iLeadDim, - pfltMatrixcX); + pfltMatrixcX);*/ + + for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) + pfltMatrixcX[iIndex1] = ctimess ( pfltMatrixX[iIndex1] , FloatComplex((float) fltCst , 0) ) ; /* E = E + cX */ diff --git a/src/matrixOperations/expm/dexpma.c b/src/matrixOperations/expm/dexpma.c index 9bfc4ea2..3bb3156b 100644 --- a/src/matrixOperations/expm/dexpma.c +++ b/src/matrixOperations/expm/dexpma.c @@ -59,7 +59,7 @@ void dexpma (double* in, double* out, int _iLeadDim){ for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) pdblMatrixA[iIndex1] = in[iIndex1] / dblS2 ; - + /* Pade approximation for exp(A)*/ /*X = A */ /*C2F(dcopy)(&iSquare, pdblMatrixA, &iOne, pdblMatrixX, &iOne );*/ @@ -89,7 +89,7 @@ void dexpma (double* in, double* out, int _iLeadDim){ for(iLoop1 = 2 ; iLoop1 <= iMax ; iLoop1++) { dblCst = dblCst * (iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1)); - + /*Temp = X*/ /*C2F(dcopy)(&iSquare, pdblMatrixX, &iOne, pdblMatrixTemp, &iOne);*/ for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) diff --git a/src/matrixOperations/expm/sexpma.c b/src/matrixOperations/expm/sexpma.c index 51f783e2..5a0b8bdf 100644 --- a/src/matrixOperations/expm/sexpma.c +++ b/src/matrixOperations/expm/sexpma.c @@ -25,6 +25,7 @@ void sexpma (float* in, float* out, int _iLeadDim){ float fltS = 0; float fltS2 = 0; float fltCst = 0.5f; + double dblCst = 0.5; float *pfltMatrixA = NULL;/*A'*/ float *pfltMatrixX = NULL;/*X*/ @@ -58,7 +59,7 @@ void sexpma (float* in, float* out, int _iLeadDim){ /*A = A./2^s*/ for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) - pfltMatrixA[iIndex1] = in[iIndex1] / fltS2 ; + pfltMatrixA[iIndex1] = in[iIndex1] / fltS2; /* Pade approximation for exp(A)*/ /*X = A */ @@ -88,8 +89,11 @@ void sexpma (float* in, float* out, int _iLeadDim){ for(iLoop1 = 2 ; iLoop1 <= iMax ; iLoop1++) { - fltCst *=(float) ( (iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1))); - + + + dblCst = dblCst * (iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1)); + + dblCst += 0 ; /*Temp = X*/ /*C2F(dcopy)(&iSquare, pfltMatrixX, &iOne, pfltMatrixTemp, &iOne);*/ for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) @@ -103,7 +107,7 @@ void sexpma (float* in, float* out, int _iLeadDim){ /*cX = c * X*/ for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) - pfltMatrixcX[iIndex1] = pfltMatrixX[iIndex1] * fltCst ; + pfltMatrixcX[iIndex1] = pfltMatrixX[iIndex1] * (float) dblCst ; /*E = E + cX*/ saddma ( out, iSquare , pfltMatrixcX , iSquare , out ) ; diff --git a/src/matrixOperations/expm/testMatrixExponential.c b/src/matrixOperations/expm/testMatrixExponential.c index 95056e74..361ec9b7 100644 --- a/src/matrixOperations/expm/testMatrixExponential.c +++ b/src/matrixOperations/expm/testMatrixExponential.c @@ -416,7 +416,7 @@ static void sexpmaTest (void ) { for ( i = 0 ; i < LEADDIM*LEADDIM ; i++ ) { printf ( "\t\t %d out : %e\tresult : %e\tassert : %e \n" , i , out[i] , result[i] , fabs ( out[i] - result[i] ) / fabs( out[i]) ) ; - /*assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 1e-6 ) ;*/ + assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 1e-5 ) ; } } @@ -443,6 +443,50 @@ static void dexpmaTest (void ) { } +static void cexpmaTest ( void) { + + int i = 0 ; + + float tRealIn [] = CRMATRIX_IN ; + float tImagIn [] = CIMATRIX_IN ; + + + + float tRealResult [] = CRMATRIX_RESULT ; + float tImagResult [] = CIMATRIX_RESULT ; + + floatComplex out[LEADDIM*LEADDIM ] ; + + floatComplex* in = FloatComplexMatrix ( tRealIn , tImagIn , LEADDIM*LEADDIM ); + floatComplex* Result = FloatComplexMatrix ( tRealResult , tImagResult ,LEADDIM*LEADDIM) ; + + cexpma ( in ,out , LEADDIM ) ; + + + /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/ + for ( i = 0 ; i < (LEADDIM*LEADDIM ) ; i++ ) + { + printf ( "\t\t %d out : %e\t + %e\t * i result : %e\t + %e\t * i assert : %e + %e \n" , + i ,creals(out[i]) , cimags(out[i]) , creals (Result[i]) , cimags (Result[i]), + fabs( creals(out[i]) - creals (Result[i]) ) / fabs (creals (out[i])) , + fabs( cimags(out[i]) - cimags (Result[i]) ) / fabs (cimags (out[i]))); + + if ( creals(out[i]) < 1e-14 && creals (Result[i]) < 1e-18 ) + assert ( 1 ) ; + else + assert ( fabs( creals(out[i]) - creals (Result[i]) ) / fabs (creals (out[i])) < 1e-4 ); + + + if ( cimags(out[i]) < 1e-14 && cimags (Result[i]) < 1e-18 ) + assert ( 1 ) ; + else + assert ( fabs( cimags(out[i]) - cimags (Result[i]) ) / fabs (cimags (out[i])) < 1e-4 ) ; + + } +} + + + static void zexpmaTest ( void) { int i = 0 ; @@ -474,13 +518,13 @@ static void zexpmaTest ( void) { if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 ) assert ( 1 ) ; else - assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-14 ); + assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 ); if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 ) assert ( 1 ) ; else - assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-14 ) ; + assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ; } } @@ -495,6 +539,10 @@ static int testExponential(void) { printf("\n\n\t>>>> Matrix Float Realt Tests\n"); sexpmaTest(); + printf("\n\n\n"); + printf("\t>>>> Matrix Float Complex Tests\n"); + cexpmaTest(); + printf("\n\n\n"); printf("\t>>>> Matrix Double Complex Tests\n"); zexpmaTest(); diff --git a/src/matrixOperations/expm/zexpma.c b/src/matrixOperations/expm/zexpma.c index 7ee85f17..9b378690 100644 --- a/src/matrixOperations/expm/zexpma.c +++ b/src/matrixOperations/expm/zexpma.c @@ -73,9 +73,12 @@ void zexpma(doubleComplex * in, doubleComplex * out, int _iLeadDim) zeyesa(pdblMatrixEye, _iLeadDim, _iLeadDim); - zmulma ( & zdblCst , 1 ,1, + /* zmulma ( & zdblCst , 1 ,1, pdblMatrixA , _iLeadDim, _iLeadDim, - pdblMatrixcA); + pdblMatrixcA);*/ + + for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) + pdblMatrixcA[iIndex1] = ztimess ( pdblMatrixA[iIndex1] , zdblCst ) ; /* cA = A * c */ @@ -107,9 +110,12 @@ void zexpma(doubleComplex * in, doubleComplex * out, int _iLeadDim) pdblMatrixX ); /* cX = c * X */ - zmulma ( & zdblCst , 1 ,1, + /* zmulma ( & zdblCst , 1 ,1, pdblMatrixX , _iLeadDim, _iLeadDim, - pdblMatrixcX); + pdblMatrixcX);*/ + + for ( iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++ ) + pdblMatrixcX[iIndex1] = ztimess ( pdblMatrixX[iIndex1] , zdblCst ) ; /* E = E + cX */ -- cgit