symmetricschurdecomposition.cpp

/* -*- 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.
*/
 
#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
#include <vector>
 
namespace QuantLib {
 
    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<>());
        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];
            }
        }
    }
 
}