<commit_before>// Copyright (C) 2019 <NAME> <<EMAIL>>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.

#ifndef BK_LDLT_H
#define BK_LDLT_H

#include <Eigen/Core>
#include <vector>
#include <stdexcept>
#include <iostream>

namespace Spectra {


template <typename Scalar = double>
class BKLDLT
{
private:
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
    typedef Eigen::Map<Vector> MapVec;

    typedef typename Matrix::Index Index;

    typedef Eigen::Matrix<Index, Eigen::Dynamic, 1> IntVector;
    typedef Eigen::Ref<Matrix> GenericMatrix;
    typedef const Eigen::Ref<const Matrix> ConstGenericMatrix;

    Index m_n;
    Vector m_data;
    std::vector<Scalar*> m_mat;  // m_mat represents a lower-triangular matrix
                                 // m_mat[i] points to the head of the i-th column
    IntVector m_perm;

    bool m_computed;

    void compute_pointer()
    {
        m_mat.clear();
        m_mat.reserve(m_n);
        Scalar* head = m_data.data();

        for(Index i = 0; i < m_n; i++)
        {
            m_mat.push_back(head);
            head += (m_n - i);
        }
    }

    void copy_data(ConstGenericMatrix& mat, int uplo = Eigen::Lower)
    {
        if(uplo == Eigen::Lower)
        {
            for(Index j = 0; j < m_n; j++)
            {
                std::copy(&mat.coeffRef(j, j), &mat.coeffRef(0, j + 1), m_mat[j]);
            }
        } else {
            Scalar* dest = m_data.data();
            for(Index i = 0; i < m_n; i++)
            {
                for(Index j = i; j < m_n; j++, dest++)
                {
                    *dest = mat.coeff(i, j);
                }
            }
        }
    }

    void print_mat() const
    {
        Matrix mat = Matrix::Zero(m_n, m_n);
        for(Index j = 0; j < m_n; j++)
        {
            std::copy(m_mat[j], m_mat[j] + (m_n - j), &mat.coeffRef(j, j));
        }
        std::cout << mat << std::endl << std::endl;
    }

    // Working on the A[k:end, k:end] submatrix
    // Exchange k <-> r
    // Assume r >= k
    void pivoting_1x1(Index k, Index r)
    {
        // No permutation
        if(k == r)
        {
            m_perm[k] = r;
            return;
        }

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // A[k, k] <-> A[r, r]
        std::swap(matp[k][0], matp[r][0]);

        // A[(r+1):end, k] <-> A[(r+1):end, r]
        std::swap_ranges(matp[k] + (r + 1 - k), matp[k + 1], matp[r] + 1);

        // A[(k+1):(r-1), k] <-> A[r, (k+1):(r-1)]
        Scalar* src = matp[k] + 1;
        for(Index j = k + 1; j < r; j++, src++)
        {
            std::swap(*src, matp[j][r - j]);
        }

        m_perm[k] = r;
    }

    // Working on the A[k:end, k:end] submatrix
    // Exchange [k+1, k] <-> [r, p]
    // Assume p >= k, r >= k+1
    void pivoting_2x2(Index k, Index r, Index p)
    {
        pivoting_1x1(k, p);
        pivoting_1x1(k + 1, r);

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // A[k+1, k] <-> A[r, k]
        std::swap(matp[k][1], matp[k][r - k]);

        // Use negative signs to indicate a 2x2 block
        m_perm[k] = -m_perm[k];
        m_perm[k + 1] = -m_perm[k + 1];
    }

    // A[r1, c1:c2] <-> A[r2, c1:c2]
    // Assume r2 >= r1 > c2 >= c1
    void interchange_rows(Index r1, Index r2, Index c1, Index c2)
    {
        if(r1 == r2)
            return;

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        for(Index j = c1; j <= c2; j++)
        {
            std::swap(matp[j][r1 - j], matp[j][r2 - j]);
        }
    }

    // lambda = |A[r, k]| = max{|A[k+1, k]|, ..., |A[end, k]|}
    // Assume k < end
    Scalar find_lambda(Index k, Index& r)
    {
        using std::abs;

        const Scalar* head = m_mat[k];
        const Scalar* end = m_mat[k + 1];
        Scalar lambda = abs(head[1]);
        r = k + 1;
        for(const Scalar* ptr = head + 2; ptr < end; ptr++)
        {
            const Scalar abs_elem = abs(*ptr);
            if(lambda < abs_elem)
            {
                lambda = abs_elem;
                r = k + (ptr - head);
            }
        }
        return lambda;
    }

    // sigma = |A[p, r]| = max {|A[k, r]|, ..., |A[end, r]|} \ {A[r, r]}
    // Assume k < r < end
    Scalar find_sigma(Index k, Index r, Index& p)
    {
        using std::abs;

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // First search A[r+1, r], ...,  A[end, r], which has the same task as find_lambda()
        Scalar sigma = find_lambda(r, p);

        // Then search A[k, r], ..., A[r-1, r], which maps to A[r, k], ..., A[r, r-1]
        for(Index j = k; j < r; j++)
        {
            const Scalar abs_elem = abs(matp[j][r - j]);
            if(sigma < abs_elem)
            {
                sigma = abs_elem;
                p = j;
            }
        }

        return sigma;
    }

    // Generate permutations and apply to A
    // Return true if the resulting pivoting is 1x1, and false if 2x2
    bool permutate_mat(Index k, const Scalar& alpha)
    {
        // matp[j][i - j] -> A[i, j], i >= j<commit_msg>fix an edge case<commit_after>// Copyright (C) 2019 <NAME> <<EMAIL>>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.

#ifndef BK_LDLT_H
#define BK_LDLT_H

#include <Eigen/Core>
#include <vector>
#include <stdexcept>
#include <iostream>

namespace Spectra {


template <typename Scalar = double>
class BKLDLT
{
private:
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
    typedef Eigen::Map<Vector> MapVec;

    typedef typename Matrix::Index Index;

    typedef Eigen::Matrix<Index, Eigen::Dynamic, 1> IntVector;
    typedef Eigen::Ref<Matrix> GenericMatrix;
    typedef const Eigen::Ref<const Matrix> ConstGenericMatrix;

    Index m_n;
    Vector m_data;
    std::vector<Scalar*> m_mat;  // m_mat represents a lower-triangular matrix
                                 // m_mat[i] points to the head of the i-th column
    IntVector m_perm;

    bool m_computed;

    void compute_pointer()
    {
        m_mat.clear();
        m_mat.reserve(m_n);
        Scalar* head = m_data.data();

        for(Index i = 0; i < m_n; i++)
        {
            m_mat.push_back(head);
            head += (m_n - i);
        }
    }

    void copy_data(ConstGenericMatrix& mat, int uplo = Eigen::Lower)
    {
        if(uplo == Eigen::Lower)
        {
            for(Index j = 0; j < m_n; j++)
            {
                std::copy(&mat.coeffRef(j, j), &mat.coeffRef(0, j + 1), m_mat[j]);
            }
        } else {
            Scalar* dest = m_data.data();
            for(Index i = 0; i < m_n; i++)
            {
                for(Index j = i; j < m_n; j++, dest++)
                {
                    *dest = mat.coeff(i, j);
                }
            }
        }
    }

    void print_mat() const
    {
        Matrix mat = Matrix::Zero(m_n, m_n);
        for(Index j = 0; j < m_n; j++)
        {
            std::copy(m_mat[j], m_mat[j] + (m_n - j), &mat.coeffRef(j, j));
        }
        std::cout << mat << std::endl << std::endl;
    }

    // Working on the A[k:end, k:end] submatrix
    // Exchange k <-> r
    // Assume r >= k
    void pivoting_1x1(Index k, Index r)
    {
        // No permutation
        if(k == r)
        {
            m_perm[k] = r;
            return;
        }

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // A[k, k] <-> A[r, r]
        std::swap(matp[k][0], matp[r][0]);

        // A[(r+1):end, k] <-> A[(r+1):end, r]
        std::swap_ranges(matp[k] + (r + 1 - k), matp[k + 1], matp[r] + 1);

        // A[(k+1):(r-1), k] <-> A[r, (k+1):(r-1)]
        Scalar* src = matp[k] + 1;
        for(Index j = k + 1; j < r; j++, src++)
        {
            std::swap(*src, matp[j][r - j]);
        }

        m_perm[k] = r;
    }

    // Working on the A[k:end, k:end] submatrix
    // Exchange [k+1, k] <-> [r, p]
    // Assume p >= k, r >= k+1
    void pivoting_2x2(Index k, Index r, Index p)
    {
        pivoting_1x1(k, p);
        pivoting_1x1(k + 1, r);

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // A[k+1, k] <-> A[r, k]
        std::swap(matp[k][1], matp[k][r - k]);

        // Use negative signs to indicate a 2x2 block
        m_perm[k] = -m_perm[k];
        m_perm[k + 1] = -m_perm[k + 1];
    }

    // A[r1, c1:c2] <-> A[r2, c1:c2]
    // Assume r2 >= r1 > c2 >= c1
    void interchange_rows(Index r1, Index r2, Index c1, Index c2)
    {
        if(r1 == r2)
            return;

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        for(Index j = c1; j <= c2; j++)
        {
            std::swap(matp[j][r1 - j], matp[j][r2 - j]);
        }
    }

    // lambda = |A[r, k]| = max{|A[k+1, k]|, ..., |A[end, k]|}
    // Assume k < end
    Scalar find_lambda(Index k, Index& r)
    {
        using std::abs;

        const Scalar* head = m_mat[k];
        const Scalar* end = m_mat[k + 1];
        Scalar lambda = abs(head[1]);
        r = k + 1;
        for(const Scalar* ptr = head + 2; ptr < end; ptr++)
        {
            const Scalar abs_elem = abs(*ptr);
            if(lambda < abs_elem)
            {
                lambda = abs_elem;
                r = k + (ptr - head);
            }
        }
        return lambda;
    }

    // sigma = |A[p, r]| = max {|A[k, r]|, ..., |A[end, r]|} \ {A[r, r]}
    // Assume k < r < end
    Scalar find_sigma(Index k, Index r, Index& p)
    {
        using std::abs;

        // matp[j][i - j] -> A[i, j], i >= j
        Scalar** matp = &m_mat.front();

        // First search A[r+1, r], ...,  A[end, r], which has the same task as find_lambda()
        // If r == end, we skip this search
        Scalar sigma = Scalar(-1);
        if(r < m_n - 1)
            sigma = find_lambda(r, p);

        // Then search A[k, r], ..., A[r-1, r], which maps to A[r, k], ..., A[r, r-1]
        for(Index j = k; j < r; j++)
        {
            const Scalar abs_elem = abs(matp[j][r - j]);
            if(sigma < abs_elem)
            {
                sigma = abs_elem;
                p = j;
            }
        }

        return sigma;
    }

    // Generate permutations and apply to A
    // Return true if the resulting pivoting is 1x1, and false if 2x2
    bool permutate_mat(Index k, const Scalar& alpha)
    {
        // matp[j][i - j] -> A[i, j], i >= j
{'fraction_non_alphanumeric': 0.171360860450278, 'fraction_numerical': 0.014401604229331875, 'mean_word_length': 1.581039755351682, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 79, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5263018', 'n_tokens_mistral': 4233, 'n_tokens_neox': 4002, 'n_words': 1074}