shells/utilities/EigenDecomp.h

Go to the documentation of this file.

// -*- C++ -*- 
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//                                   Fehmi Cirak
//                        California Institute of Technology
//                           (C) 2004 All Rights Reserved
//
// <LicenseText>
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#ifndef EIGENDECOMP_H
#define EIGENDECOMP_H

#include "CholeskyDecomp.h"

#include <algorithm>
#include <numeric>
#include <limits>
#include <cassert>
#include <cmath>
#include <cstdlib>


namespace utilities {
    template<int ND>
    double forwardIteration(double km[ND][ND], double mm[ND][ND]); 
}


// implementation 
namespace utilities {

// helpers
namespace {
    template<int ND>
    void multMatrixVector(double km[ND][ND], double x[ND], double y[ND])
    {
        for (int i=0; i<ND; ++i) {
            double sum = 0.0;
            for (int k=0; k<ND; ++k) {
                sum += km[i][k]*x[k];
                }
            y[i] = sum;
        }       
        return;
    }   
    
    struct GenRandomNumbers {       
        GenRandomNumbers(int seed=5) {std::srand(seed);}            

        double operator()() {
            return (static_cast<double>(std::rand())/static_cast<double>(RAND_MAX));
        }
    };
}


template<int ND>
double forwardIteration(double km[ND][ND], double mm[ND][ND]) 
{
    //
    // compute the maximum eigenvalue of
    // (km - \rho mm) \phi = 0
    //  
    // note that mm has to be positive definite
    // 
    double x[ND];
    
    bool notConverged = true;    
    int step=0;
    const int maxStep=1000;
    
    std::generate_n(x, ND, GenRandomNumbers());
    
    double y[ND], ybar[ND];    
    multMatrixVector(km, x, y);
    
    double diag[ND];
    utilities::choleskyDecomposition(mm, diag);
        
    double rho=1.0;    
    while ((notConverged)&&(step<maxStep)) {    
        utilities::choleskySolve(mm, diag, y, x);       
        multMatrixVector(km, x, ybar);
        
        double rhoPr = rho;
        rho = std::inner_product(x, x+ND, ybar, 0.0);
        double tmp = std::inner_product(x, x+ND, y, 0.0);       
        rho /= tmp;
        
        if (std::fabs(rho-rhoPr)<(std::numeric_limits<double>().epsilon()*std::fabs(rho))) 
            notConverged = false;
            
        if (tmp<0.0) assert(false);
        tmp = 1./std::sqrt(tmp);        
        std::transform(ybar, ybar+ND, y, std::bind2nd(std::multiplies<double>(), tmp)); 
        
        ++step;
        assert(step<maxStep);
    }
    
    return rho;
}    

} // namespace utilities

#endif

---