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


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 


32  void 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++ = gs*(h+g*tau); *aq++ = h+s*(gh*tau);


80  }


81  int ip = p+1; j = qip; ap += ip++; aq++;


82  while (j)


83  {


84  g = *ap; h = *aq;


85  *ap = gs*(h+g*tau); *aq++ = h+s*(gh*tau);


86  ap += ip++;


87  }


88  if (q < n1) // last loop is nonempty


89  {


90  int iq = q+1; j = niq; ap += ip++; aq += iq++;


91  for (;;)


92  {


93  g = *ap; h = *aq;


94  *ap = gs*(h+g*tau); *aq = h+s*(gh*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 


115  void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)


116  { REPORT SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }


117 


118  void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)


119  { REPORT Matrix V; Jacobi(X,D,A,V,false); }


120 


121  void 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  ///@}

