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