diff --git a/test/Makefile.am b/test/Makefile.am index 36739c28ae0c47f201a351967f64827540d1109c..7109770fde2a252c9d8f5263e92e08ab2493f9e1 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -5,7 +5,8 @@ LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) AM_CPPFLAGS += $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) -Wall check_PROGRAMS = frameinvariancetest rotationtest fdcheck localgeodesicfefunctiontest \ - harmonicenergytest averagedistanceassemblertest unitvectortest + harmonicenergytest averagedistanceassemblertest unitvectortest \ + svdtest frameinvariancetest_SOURCES = frameinvariancetest.cc @@ -21,6 +22,8 @@ averagedistanceassemblertest_SOURCES = averagedistanceassemblertest.cc unitvectortest_SOURCES = unitvectortest.cc +svdtest_SOURCES = svdtest.cc + # don't follow the full GNU-standard # we need automake 1.5 AUTOMAKE_OPTIONS = foreign 1.5 diff --git a/test/svdtest.cc b/test/svdtest.cc new file mode 100644 index 0000000000000000000000000000000000000000..1d49302a4078d064b536fb399b46190b521ffc25 --- /dev/null +++ b/test/svdtest.cc @@ -0,0 +1,57 @@ +#include <iostream> + +#include <dune/gfe/svd.hh> + +#include "multiindex.hh" + +using namespace Dune; + +int main() +{ + const int N=3; + const int M=3; + + FieldMatrix<double,N,M> testMatrix; + + int nTestValues = 5; + double maxDiff = 0; + + MultiIndex<N*M> index(nTestValues); + int numIndices = index.cycle(); + + for (int i=0; i<numIndices; i++, ++index) { + + for (int j=0; j<N; j++) + for (int k=0; k<M; k++) { + testMatrix[j][k] = std::exp(double(2*index[j*M+k]-double(nTestValues))) - std::exp(double(nTestValues)); + //std::cout << std::exp(2*index[j*M+k]-double(nTestValues)) << std::endl; + } + //std::cout << "testMatrix:\n" << testMatrix << std::endl; + //exit(0); + FieldVector<double,N> w; + FieldMatrix<double,N,N> v; + + FieldMatrix<double,N,M> testMatrixBackup = testMatrix; + svdcmp(testMatrix,w,v); + + // Multiply the three matrices to see whether we get the original one back + FieldMatrix<double,N,M> product(0); + + for (int j=0; j<N; j++) + for (int k=0; k<M; k++) + for (int l=0; l<N; l++) + product[j][k] += testMatrix[j][l] * w[l] * v[k][l]; + + product -= testMatrixBackup; + //std::cout << product << std::endl; + if (maxDiff < product.infinity_norm()) { + maxDiff = product.infinity_norm(); + std::cout << "max diff: " << maxDiff << std::endl; + std::cout << "testMatrix:\n" << testMatrixBackup << std::endl; + std::cout << "difference:\n" << product << std::endl; + exit(0); + } + + } + +} \ No newline at end of file