source: ntrip/branches/BNC_2.12/newmat/jacobi.cpp

Last change on this file was 2013, checked in by mervart, 13 years ago

* empty log message *

File size: 3.8 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file jacobi.cpp
5/// Eigen value decomposition using Jacobi method.
6
7
8// Copyright (C) 1991,2,3,4: R B Davies
9
10
11//#define WANT_STREAM
12
13
14#define WANT_MATH
15
16#include "include.h"
17#include "newmatap.h"
18#include "precisio.h"
19#include "newmatrm.h"
20
21#ifdef use_namespace
22namespace NEWMAT {
23#endif
24
25#ifdef DO_REPORT
26#define REPORT { static ExeCounter ExeCount(__LINE__,18); ++ExeCount; }
27#else
28#define REPORT {}
29#endif
30
31
32void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
33 Matrix& V, bool eivec)
34{
35 Real epsilon = FloatingPointPrecision::Epsilon();
36 Tracer et("Jacobi");
37 REPORT
38 int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.resize(n); A = X;
39 if (eivec) { REPORT V.resize(n,n); D = 1.0; V = D; }
40 B << A; D = B; Z = 0.0; A.Inject(Z);
41 bool converged = false;
42 for (int i=1; i<=50; i++)
43 {
44 Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
45 while (p--) sm += fabs(*a++); // have previously zeroed diags
46 if (sm==0.0) { REPORT converged = true; break; }
47 Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
48 for (p = 0; p < n; p++)
49 {
50 Real* ap1 = a + (p*(p+1))/2;
51 Real& zp = Z.element(p); Real& dp = D.element(p);
52 for (int q = p+1; q < n; q++)
53 {
54 Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
55 Real& zq = Z.element(q); Real& dq = D.element(q);
56 Real& apq = A.element(q,p);
57 Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);
58
59 if (i>4 && g < epsilon*adp && g < epsilon*adq) { REPORT apq = 0.0; }
60 else if (fabs(apq) > tresh)
61 {
62 REPORT
63 Real t; Real h = dq - dp; Real ah = fabs(h);
64 if (g < epsilon*ah) { REPORT t = apq / h; }
65 else
66 {
67 REPORT
68 Real theta = 0.5 * h / apq;
69 t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
70 if (theta<0.0) { REPORT t = -t; }
71 }
72 Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
73 Real tau = s / (1.0 + c); h = t * apq;
74 zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
75 int j = p;
76 while (j--)
77 {
78 g = *ap; h = *aq;
79 *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
80 }
81 int ip = p+1; j = q-ip; ap += ip++; aq++;
82 while (j--)
83 {
84 g = *ap; h = *aq;
85 *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
86 ap += ip++;
87 }
88 if (q < n-1) // last loop is non-empty
89 {
90 int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
91 for (;;)
92 {
93 g = *ap; h = *aq;
94 *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
95 if (!(--j)) break;
96 ap += ip++; aq += iq++;
97 }
98 }
99 if (eivec)
100 {
101 REPORT
102 RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
103 Rotate(VP, VQ, tau, s);
104 }
105 }
106 }
107 }
108 B = B + Z; D = B; Z = 0.0;
109 }
110 if (!converged) Throw(ConvergenceException(X));
111 if (eivec) SortSV(D, V, true);
112 else SortAscending(D);
113}
114
115void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)
116{ REPORT SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }
117
118void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)
119{ REPORT Matrix V; Jacobi(X,D,A,V,false); }
120
121void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, Matrix& V)
122{ REPORT SymmetricMatrix A; Jacobi(X,D,A,V,true); }
123
124
125#ifdef use_namespace
126}
127#endif
128
129
130///@}
Note: See TracBrowser for help on using the repository browser.