text: /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
  
  /*
   Copyright (C) 2003 Ferdinando Ametrano
   Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
  
   This file is part of QuantLib, a free-software/open-source library
   for financial quantitative analysts and developers - http://quantlib.org/
  
   QuantLib is free software: you can redistribute it and/or modify it
   under the terms of the QuantLib license.  You should have received a
   copy of the license along with this program; if not, please email
   <quantlib-dev@lists.sf.net>. The license is also available online at
   <http://quantlib.org/license.shtml>.
  
   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE.  See the license for more details.
  */
  
  /*! \file symmetricschurdecomposition.hpp
      \brief Eigenvalues/eigenvectors of a real symmetric matrix
  */
  
  #ifndef quantlib_math_jacobi_decomposition_h
  #define quantlib_math_jacobi_decomposition_h
  
  #include <ql/math/matrix.hpp>
  #include <vector>
  
  namespace QuantLib {
  
      //! symmetric threshold Jacobi algorithm.
      /*! Given a real symmetric matrix S, the Schur decomposition
          finds the eigenvalues and eigenvectors of S. If D is the
          diagonal matrix formed by the eigenvalues and U the
          unitarian matrix of the eigenvectors we can write the
          Schur decomposition as
          \f[ S = U \cdot D \cdot U^T \, ,\f]
          where \f$ \cdot \f$ is the standard matrix product
          and  \f$ ^T  \f$ is the transpose operator.
          This class implements the Schur decomposition using the
          symmetric threshold Jacobi algorithm. For details on the
          different Jacobi transfomations see "Matrix computation,"
          second edition, by Golub and Van Loan,
          The Johns Hopkins University Press
  
          \test the correctness of the returned values is tested by
                checking their properties.
      */
      class SymmetricSchurDecomposition {
        public:
          /*! \pre s must be symmetric */
          SymmetricSchurDecomposition(const Matrix &s);
          const Array& eigenvalues() const { return diagonal_; }
          const Matrix& eigenvectors() const { return eigenVectors_; }
        private:
          Array diagonal_;
          Matrix eigenVectors_;
          void jacobiRotate_(Matrix & m, Real rot, Real dil,
                             Size j1, Size k1, Size j2, Size k2) const;
      };
  
  
      // inline definitions
  
      //! This routines implements the Jacobi, a.k.a. Givens, rotation
      inline void SymmetricSchurDecomposition::jacobiRotate_(
                                        Matrix &m, Real rot, Real dil, Size j1,
                                        Size k1, Size j2, Size k2) const {
          Real x1, x2;
          x1 = m[j1][k1];
          x2 = m[j2][k2];
          m[j1][k1] = x1 - dil*(x2 + x1*rot);
          m[j2][k2] = x2 + dil*(x1 - x2*rot);
      }
  
      // implementation
  
      inline SymmetricSchurDecomposition::SymmetricSchurDecomposition(const Matrix & s)
      : diagonal_(s.rows()), eigenVectors_(s.rows(), s.columns(), 0.0) {
  
          QL_REQUIRE(s.rows() > 0 && s.columns() > 0, "null matrix given");
          QL_REQUIRE(s.rows()==s.columns(), "input matrix must be square");
  
          Size size = s.rows();
          for (Size q=0; q<size; q++) {
              diagonal_[q] = s[q][q];
              eigenVectors_[q][q] = 1.0;
          }
          Matrix ss = s;
  
          std::vector<Real> tmpDiag(diagonal_.begin(), diagonal_.end());
          std::vector<Real> tmpAccumulate(size, 0.0);
          Real threshold, epsPrec = 1e-15;
          bool keeplooping = true;
          Size maxIterations = 100, ite = 1;
          do {
              //main loop
              Real sum = 0;
              for (Size a=0; a<size-1; a++) {
                  for (Size b=a+1; b<size; b++) {
                      sum += std::fabs(ss[a][b]);
                  }
              }
  
              if (sum==0) {
                  keeplooping = false;
              } else {
                  /* To speed up computation a threshold is introduced to
                     make sure it is worthy to perform the Jacobi rotation
                  */
                  if (ite<5) threshold = 0.2*sum/(size*size);
                  else       threshold = 0.0;
  
                  Size j, k, l;
                  for (j=0; j<size-1; j++) {
                      for (k=j+1; k<size; k++) {
                          Real sine, rho, cosin, heig, tang, beta;
                          Real smll = std::fabs(ss[j][k]);
                          if(ite> 5 &&
                             smll<epsPrec*std::fabs(diagonal_[j]) &&
                             smll<epsPrec*std::fabs(diagonal_[k])) {
                                  ss[j][k] = 0;
                          } else if (std::fabs(ss[j][k])>threshold) {
                              heig = diagonal_[k]-diagonal_[j];
                              if (smll<epsPrec*std::fabs(heig)) {
                                  tang = ss[j][k]/heig;
                              } else {
                                  beta = 0.5*heig/ss[j][k];
                                  tang = 1.0/(std::fabs(beta)+
                                      std::sqrt(1+beta*beta));
                                  if (beta<0)
                                      tang = -tang;
                              }
                              cosin = 1/std::sqrt(1+tang*tang);
                              sine = tang*cosin;
                              rho = sine/(1+cosin);
                              heig = tang*ss[j][k];
                              tmpAccumulate[j] -= heig;
                              tmpAccumulate[k] += heig;
                              diagonal_[j] -= heig;
                              diagonal_[k] += heig;
                              ss[j][k] = 0.0;
                              for (l=0; l+1<=j; l++)
                                  jacobiRotate_(ss, rho, sine, l, j, l, k);
                              for (l=j+1; l<=k-1; l++)
                                  jacobiRotate_(ss, rho, sine, j, l, l, k);
                              for (l=k+1; l<size; l++)
                                  jacobiRotate_(ss, rho, sine, j, l, k, l);
                              for (l=0;   l<size; l++)
                                  jacobiRotate_(eigenVectors_,
                                                    rho, sine, l, j, l, k);
                          }
                      }
                  }
                  for (k=0; k<size; k++) {
                      tmpDiag[k] += tmpAccumulate[k];
                      diagonal_[k] = tmpDiag[k];
                      tmpAccumulate[k] = 0.0;
                  }
              }
          } while (++ite<=maxIterations && keeplooping);
  
          QL_ENSURE(ite<=maxIterations,
                    "Too many iterations (" << maxIterations << ") reached");
  
  
          // sort (eigenvalues, eigenvectors)
          std::vector<std::pair<Real, std::vector<Real> > > temp(size);
          std::vector<Real> eigenVector(size);
          Size row, col;
          for (col=0; col<size; col++) {
              std::copy(eigenVectors_.column_begin(col),
                        eigenVectors_.column_end(col), eigenVector.begin());
              temp[col] = std::make_pair(diagonal_[col], eigenVector);
          }
          std::sort(temp.begin(), temp.end(),
              std::greater<std::pair<Real, std::vector<Real> > >());
          Real maxEv = temp[0].first;
          for (col=0; col<size; col++) {
              // check for round-off errors
              diagonal_[col] =
                  (std::fabs(temp[col].first/maxEv)<1e-16 ? 0.0 :
                                                            temp[col].first);
              Real sign = 1.0;
              if (temp[col].second[0]<0.0)
                  sign = -1.0;
              for (row=0; row<size; row++) {
                  eigenVectors_[row][col] = sign * temp[col].second[row];
              }
          }
      }
      
  
  }
  
  
  #endif
  
  
  