Unit test for the singular value decomposition

[[Imported from SVN: r7300]]

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

test/svdtest.cc

+#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