shells/utilities/EigenDecomp.h

Go to the documentation of this file.
00001 // -*- C++ -*- 
00002 //
00003 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00004 //
00005 //                                   Fehmi Cirak
00006 //                        California Institute of Technology
00007 //                           (C) 2004 All Rights Reserved
00008 //
00009 // <LicenseText>
00010 //
00011 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00012 //
00013 #ifndef EIGENDECOMP_H
00014 #define EIGENDECOMP_H
00015 
00016 #include "CholeskyDecomp.h"
00017 
00018 #include <algorithm>
00019 #include <numeric>
00020 #include <limits>
00021 #include <cassert>
00022 #include <cmath>
00023 #include <cstdlib>
00024 
00025 
00026 namespace utilities {
00027     template<int ND>
00028     double forwardIteration(double km[ND][ND], double mm[ND][ND]); 
00029 }
00030 
00031 
00032 // implementation 
00033 namespace utilities {
00034 
00035 // helpers
00036 namespace {
00037     template<int ND>
00038     void multMatrixVector(double km[ND][ND], double x[ND], double y[ND])
00039     {
00040         for (int i=0; i<ND; ++i) {
00041             double sum = 0.0;
00042             for (int k=0; k<ND; ++k) {
00043                 sum += km[i][k]*x[k];
00044                 }
00045             y[i] = sum;
00046         }       
00047         return;
00048     }   
00049     
00050     struct GenRandomNumbers {       
00051         GenRandomNumbers(int seed=5) {std::srand(seed);}            
00052 
00053         double operator()() {
00054             return (static_cast<double>(std::rand())/static_cast<double>(RAND_MAX));
00055         }
00056     };
00057 }
00058 
00059 
00060 template<int ND>
00061 double forwardIteration(double km[ND][ND], double mm[ND][ND]) 
00062 {
00063     //
00064     // compute the maximum eigenvalue of
00065     // (km - \rho mm) \phi = 0
00066     //  
00067     // note that mm has to be positive definite
00068     // 
00069     double x[ND];
00070     
00071     bool notConverged = true;    
00072     int step=0;
00073     const int maxStep=1000;
00074     
00075     std::generate_n(x, ND, GenRandomNumbers());
00076     
00077     double y[ND], ybar[ND];    
00078     multMatrixVector(km, x, y);
00079     
00080     double diag[ND];
00081     utilities::choleskyDecomposition(mm, diag);
00082         
00083     double rho=1.0;    
00084     while ((notConverged)&&(step<maxStep)) {    
00085         utilities::choleskySolve(mm, diag, y, x);       
00086         multMatrixVector(km, x, ybar);
00087         
00088         double rhoPr = rho;
00089         rho = std::inner_product(x, x+ND, ybar, 0.0);
00090         double tmp = std::inner_product(x, x+ND, y, 0.0);       
00091         rho /= tmp;
00092         
00093         if (std::fabs(rho-rhoPr)<(std::numeric_limits<double>().epsilon()*std::fabs(rho))) 
00094             notConverged = false;
00095             
00096         if (tmp<0.0) assert(false);
00097         tmp = 1./std::sqrt(tmp);        
00098         std::transform(ybar, ybar+ND, y, std::bind2nd(std::multiplies<double>(), tmp)); 
00099         
00100         ++step;
00101         assert(step<maxStep);
00102     }
00103     
00104     return rho;
00105 }    
00106 
00107 } // namespace utilities
00108 
00109 #endif
Generated on Thu Jun 30 02:18:56 2016 for SFC Thin-Shell Finite Element Solver by 1.6.3