syev.hpp

Go to the documentation of this file.

//===========================================================================
/*!
 * 
 *
 * \brief      Contains the lapack bindings for the symmetric eigenvalue problem syev.
 *
 * \author      O. Krause
 * \date        2010
 *
 *
 * \par Copyright 1995-2015 Shark Development Team
 * 
 * <BR><HR>
 * This file is part of Shark.
 * <http://image.diku.dk/shark/>
 * 
 * Shark is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published 
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Shark 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
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with Shark.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
//===========================================================================
#ifndef REMORA_KERNELS_DEFAULT_SYEV_HPP
#define REMORA_KERNELS_DEFAULT_SYEV_HPP

#include "../../detail/traits.hpp"

namespace remora{ namespace bindings {
    
template <typename MatA, typename V>
void eigensort
(
    matrix_expression<MatA, cpu_tag>& matA,
    vector_expression<V, cpu_tag>& eigenValues
){
    REMORA_SIZE_CHECK(eigenValues().size() == matA().size1());
    REMORA_SIZE_CHECK(matA().size1() == matA().size2());

    std::size_t n = eigenValues().size();

    for (std::size_t i = 0; i < n - 1; i++)
    {
        std::size_t l = arg_max(subrange(eigenValues,i,n))+i;
        if (l != i) {
                //switch position of i's eigenvalue and the largest remaining eigenvalue
            std::swap(eigenValues()( l ),eigenValues()( i ));
            //switch postions of corresponding eigenvectors
            for (std::size_t j = 0; j < n; j++) {
                std::swap(matA()( j , l ), matA()( j , i ));
            }
        }
    }
}

template <typename MatA, typename V>
void syev(
    matrix_expression<MatA, cpu_tag>& vmatA,
    vector_expression<V, cpu_tag>& dvecA
) {
    REMORA_SIZE_CHECK(vmatA().size1() == vmatA().size2());
    REMORA_SIZE_CHECK(vmatA().size1() == dvecA().size());
    
    const std::size_t maxIterC = 50;
    std::size_t n = vmatA().size1();
    
    vector<double> odvecA(n,0.0);

    std::size_t j, k, l, m;
    double   b, c, f, g, h, hh, p, r, s, scale;


    //
    // reduction to tridiagonal form
    //
    for (std::size_t i = n; i-- > 1;) {
        h = 0.0;
        scale = 0.0;

        if (i > 1) {
            // scale row
            for (std::size_t k = 0; k < i; k++) {
                scale += std::fabs(vmatA()(i, k));
            }
        }

        if (scale == 0.0) {
            odvecA(i) = vmatA()(i, i-1);
        }
        else {
            for (k = 0; k < i; k++) {
                vmatA()(i, k) /= scale;
                h += vmatA()(i, k) * vmatA()(i, k);
            }

            f           = vmatA()(i, i-1);
            g           = f > 0 ? -std::sqrt(h) : std::sqrt(h);
            odvecA(i)       = scale * g;
            h          -= f * g;
            vmatA()(i, i-1) = f - g;
            f           = 0.0;

            for (j = 0; j < i; j++) {
                vmatA()(j, i) = vmatA()(i, j) / (scale * h);
                g = 0.0;

                // form element of a*u
                for (k = 0; k <= j; k++) {
                    g += vmatA()(j, k) * vmatA()(i, k);
                }

                for (k = j + 1; k < i; k++) {
                    g += vmatA()(k, j) * vmatA()(i, k);
                }

                // form element of p
                f += (odvecA(j) = g / h) * vmatA()(i, j);
            }

            hh = f / (h + h);

            // form reduced a
            for (j = 0; j < i; j++) {
                f     = vmatA()(i, j);
                g     = odvecA(j) - hh * f;
                odvecA(j) = g;

                for (k = 0; k <= j; k++) {
                    vmatA()(j, k) -= f * odvecA(k) + g * vmatA()(i, k);
                }
            }

            for (k = i; k--;) {
                vmatA()(i, k) *= scale;
            }
        }

        dvecA()(i) = h;
    }

    dvecA()(0) = odvecA(0) = 0.0;

    // accumulation of transformation matrices
    for (std::size_t i = 0; i < n; i++) {
        if (dvecA()(i)) {
            for (j = 0; j < i; j++) {
                g = 0.0;

                for (k = 0; k < i; k++) {
                    g += vmatA()(i, k) * vmatA()(k, j);
                }

                for (k = 0; k < i; k++) {
                    vmatA()(k, j) -= g * vmatA()(k, i);
                }
            }
        }

        dvecA()(i)     = vmatA()(i, i);
        vmatA()(i, i) = 1.0;

        for (j = 0; j < i; j++) {
            vmatA()(i, j) = vmatA()(j, i) = 0.0;
        }
    }

    //
    // eigenvalues from tridiagonal form
    //
    if (n <= 1) {
        return;
    }

    for (std::size_t i = 1; i < n; i++) {
        odvecA(i-1) = odvecA(i);
    }

    odvecA(n-1) = 0.0;

    for (l = 0; l < n; l++) {
        j = 0;

        do {
            // look for small sub-diagonal element
            for (m = l; m < n - 1; m++) {
                s = std::fabs(dvecA()(m)) + std::fabs(dvecA()(m+1));
                if (std::fabs(odvecA(m)) + s == s) {
                    break;
                }
            }

            p = dvecA()(l);

            if (m != l) {
                if (j++ == maxIterC)
                    throw std::invalid_argument("too many iterations in eigendecomposition");

                // form shift
                g = (dvecA()(l+1) - p) / (2.0 * odvecA(l));
                r = std::sqrt(g * g + 1.0);
                g = dvecA()(m) - p + odvecA(l) / (g + ((g) > 0 ? std::fabs(r) : -std::fabs(r)));
                s = c = 1.0;
                p = 0.0;

                for (std::size_t i = m; i-- > l;) {
                    f = s * odvecA(i);
                    b = c * odvecA(i);

                    if (std::fabs(f) >= std::fabs(g)) {
                        c       = g / f;
                        r       = std::sqrt(c * c + 1.0);
                        odvecA(i+1) = f * r;
                        s       = 1.0 / r;
                        c      *= s;
                    }
                    else {
                        s       = f / g;
                        r       = std::sqrt(s * s + 1.0);
                        odvecA(i+1) = g * r;
                        c       = 1.0 / r;
                        s      *= c;
                    }

                    g       = dvecA()(i+1) - p;
                    r       = (dvecA()(i) - g) * s + 2.0 * c * b;
                    p       = s * r;
                    dvecA()(i+1) = g + p;
                    g       = c * r - b;

                    // form vector
                    for (k = 0; k < n; k++) {
                        f           = vmatA()(k, i+1);
                        vmatA()(k, i+1) = s * vmatA()(k, i) + c * f;
                        vmatA()(k, i  ) = c * vmatA()(k, i) - s * f;
                    }
                }

                dvecA()(l) -= p;
                odvecA(l)  = g;
                odvecA(m)  = 0.0;
            }
        }
        while (m != l);
    }

    //
    // sorting eigenvalues
    //
    eigensort(vmatA, dvecA);

    //
    // normalizing eigenvectors
    //
    for (std::size_t j = n; j--;) {
        s = 0.0;
        for (std::size_t i = n; i--;) {
            s += vmatA()(i, j) * vmatA()(i, j);
        }
        s = std::sqrt(s);

        for (std::size_t i = n; i--;) {
            vmatA()(i, j) /= s;
        }
    }
}

/** @}*/

}}

#endif