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 index(N*M, 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);
        }
    
    }
    
}