jacobi.cpp
Go to the documentation of this file.
1 
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
22 namespace 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 
33  Matrix& V, bool eivec)
34 {
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 
116 { REPORT SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }
117 
119 { REPORT Matrix V; Jacobi(X,D,A,V,false); }
120 
122 { REPORT SymmetricMatrix A; Jacobi(X,D,A,V,true); }
123 
124 
125 #ifdef use_namespace
126 }
127 #endif
128 
129 
void resize(int)
Definition: newmat4.cpp:312
Real & element(int, int)
Definition: newmat6.cpp:748
void Jacobi(const SymmetricMatrix &X, DiagonalMatrix &D, SymmetricMatrix &A, Matrix &V, bool eivec)
Definition: jacobi.cpp:32
void Rotate(RectMatrixCol &U, RectMatrixCol &V, Real tau, Real s)
Definition: newmatrm.cpp:159
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
double Real
Definition: include.h:307
#define REPORT
Definition: jacobi.cpp:28
void SortAscending(GeneralMatrix &gm)
Definition: newmatap.h:137
Real square(Real x)
Definition: bandmat.cpp:672
int Nrows() const
Definition: newmat.h:494
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
FloatVector FloatVector * a
Diagonal matrix.
Definition: newmat.h:896
virtual void resize(int, int)
Definition: newmat4.cpp:289
static Real Epsilon()
Definition: precisio.h:231
Covergence failure exception.
Definition: newmat.h:1922
void SortSV(DiagonalMatrix &D, Matrix &U, bool ascending=false)
Definition: sort.cpp:194
const double epsilon
Definition: utils.h:126
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jun 7 2019 22:06:44