Subversion Repositories gelsvn

Rev

Rev 501 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 501 Rev 595
Line -... Line 1...
- 
 
 1 
/* ----------------------------------------------------------------------- *
- 
 
 2 
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
- 
 
 3 
 * Copyright (C) the authors and DTU Informatics
- 
 
 4 
 * For license and list of authors, see ../../doc/intro.pdf
- 
 
 5 
 * ----------------------------------------------------------------------- */
- 
 
 6 
 
1 
#include "eigensolution.h"
 7 
#include "eigensolution.h"
2 
 
 8 
 
3 
#include "Mat2x2f.h"
 9 
#include "Mat2x2f.h"
4 
#include "Mat3x3f.h"
 10 
#include "Mat3x3f.h"
5 
#include "Mat4x4f.h"
 11 
#include "Mat4x4f.h"
Line 18... Line 24...
18 
    const double EV_THRESH = 1e-6;
 24 
    const double EV_THRESH = 1e-6;
19 
}
 25 
}
20 
 
 26 
 
21 
namespace CGLA
 27 
namespace CGLA
22 
{
 28 
{
23 
		template <class MT>
 29 
    template <class MT>
24 
		int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
 30 
    int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
25 
		{
 31 
    {
26 
				L = MT(0);
 32 
        L = MT(0);
27 
 
 33
28 
				typedef typename MT::VectorType VT;
 34 
        typedef typename MT::VectorType VT;
29 
				MT A = Ap;
 35 
        MT A = Ap;
30 
				unsigned int n = s_min(MT::get_v_dim(), max_sol);
 36 
        unsigned int n = s_min(MT::get_v_dim(), max_sol);
31 
 
 37
- 
 
 38 
        gel_srand(0);
32 
				for(unsigned int i=0;i<n;++i)
 39 
        for(unsigned int i=0;i<n;++i)
33 
				{
 40 
        {
34 
						// Seed the eigenvector estimate
 41 
            // Seed the eigenvector estimate
35 
						VT q(1);
 42 
            VT q;
- 
 
 43 
            for (int j=0; j<MT::get_v_dim(); ++j)
- 
 
 44 
                q[j] = gel_rand()/static_cast<double>(GEL_RAND_MAX);
- 
 
 45
36 
						q.normalize();
 46 
            q.normalize();
37 
						double l=123.0,l_old;
 47 
            double l=123,l_old;
38 
 
 48
39 
						// As long as we haven't reached the max iterations and the
 49 
            // As long as we haven't reached the max iterations and the
40 
						// eigenvalue has not converged, do
 50 
            // eigenvalue has not converged, do
41 
						unsigned int k=0;
 51 
            unsigned int k=0;
42 
						do
 52 
            do
43 
						{
 53 
            {
44 
							const VT z = A * q;
 54 
                const VT z = A * q;
45 
							double z_len = length(z);
 55 
                double z_len = length(z);
46
56
47 
							if(z_len < EV_THRESH) return i;
 57 
                if(z_len < EV_THRESH) return i;
48
58
49 
							l_old = l;
 59 
                l_old = l;
50 
							l = dot(q, z)>0 ? z_len : -z_len;
 60 
                l = dot(q, z)>0 ? z_len : -z_len;
51 
							q = z/z_len;
 61 
                q = z/z_len;
52
62
53 
							if(++k==KMAX)
 63 
                if(++k==KMAX)
54 
								return i;
 64 
                    return i;
55 
						}
 65 
            }
56 
						while((fabs(l-l_old) > fabs(EV_THRESH * l)) || k<2);
 66 
            while((fabs(l-l_old) > fabs(EV_THRESH * l)) || k<2);
57
67
58 
						// Update the solution by adding the eigenvector to Q and
 68 
            // Update the solution by adding the eigenvector to Q and
59 
						// the eigenvalue to the diagonal of L.
 69 
            // the eigenvalue to the diagonal of L.
60 
						Q[i] = q;
 70 
            Q[i] = q;
61 
						L[i][i] = l;
 71 
            L[i][i] = l;
62 
 
 72
63 
						// Update A by subtracting the subspace represented by the
73 
            // Update A by subtracting the subspace represented by the
64 
						// eigensolution just found. This is called the method of
74 
            // eigensolution just found. This is called the method of
65 
						// deflation.
 75 
            // deflation.
66 
						MT B;
 76 
            MT B;
67 
						outer_product(q,q,B);
 77 
            outer_product(q,q,B);
68 
						A = A - l * B;
 78 
            A = A - l * B;
69 
				}
 79 
        }
70 
				return n;
 80 
        return n;
71 
		}
 81 
    }
72 
 
 82
73 
		/* There is no reason to put this template in a header file, since
83 
    /* There is no reason to put this template in a header file, since
74 
			 we will only use it on matrices defined in CGLA. Instead, we
84 
     we will only use it on matrices defined in CGLA. Instead, we
75 
			 explicitly instantiate the function for the square matrices
 85 
     explicitly instantiate the function for the square matrices
76 
			 of CGLA */
 86 
     of CGLA */
77 
		template int power_eigensolution<Mat2x2f>(const Mat2x2f&,
 87 
    template int power_eigensolution<Mat2x2f>(const Mat2x2f&,
78 
																							Mat2x2f&,Mat2x2f&,unsigned int);
 88 
                                              Mat2x2f&,Mat2x2f&,unsigned int);
79
89
80 
		template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
 90 
    template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
81 
																							Mat3x3f&,Mat3x3f&,unsigned int);
 91 
                                              Mat3x3f&,Mat3x3f&,unsigned int);
82 
		template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
 92 
    template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
83 
																							Mat4x4f&,Mat4x4f&,unsigned int);
 93 
                                              Mat4x4f&,Mat4x4f&,unsigned int);
84 
		template int power_eigensolution<Mat2x2d>(const Mat2x2d&,
 94 
    template int power_eigensolution<Mat2x2d>(const Mat2x2d&,
85 
																							Mat2x2d&,Mat2x2d&,unsigned int);
 95 
                                              Mat2x2d&,Mat2x2d&,unsigned int);
86
96
87 
		template int power_eigensolution<Mat3x3d>(const Mat3x3d&,
 97 
    template int power_eigensolution<Mat3x3d>(const Mat3x3d&,
88 
																							Mat3x3d&,Mat3x3d&,unsigned int);
 98 
                                              Mat3x3d&,Mat3x3d&,unsigned int);
89 
		template int power_eigensolution<Mat4x4d>(const Mat4x4d&,
 99 
    template int power_eigensolution<Mat4x4d>(const Mat4x4d&,
90 
																							Mat4x4d&,Mat4x4d&,unsigned int);
 100 
                                              Mat4x4d&,Mat4x4d&,unsigned int);
91 
}
 101 
}
92 
 
 102