<commit_before>/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, version 1.0 RC 1        *
*                (c) 2006-2011 MGH, INRIA, USTL, UJF, CNRS                    *
*                                                                             *
* This library 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 2.1 of the License, or (at     *
* your option) any later version.                                             *
*                                                                             *
* This library 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 this library; if not, write to the Free Software Foundation,     *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.          *
*******************************************************************************
*                               SOFA :: Modules                               *
*                                                                             *
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: <EMAIL>                             *
******************************************************************************/
#ifndef SOFA_COMPONENT_LINEARSOLVER_SPARSELDLSOLVERIMPL_H
#define SOFA_COMPONENT_LINEARSOLVER_SPARSELDLSOLVERIMPL_H

#include <sofa/core/behavior/LinearSolver.h>
#include <sofa/component/linearsolver/MatrixLinearSolver.h>

#ifdef SOFA_HAVE_METIS
extern "C" {
#include <metis.h>
}
#endif

namespace sofa
{

namespace component
{

namespace linearsolver
{

//defaut structure for a LDL factorization
template<class VecInt,class VecReal>
class SpaseLDLImplInvertData : public MatrixInvertData {
public :
    int n, P_nnz, L_nnz;
    VecInt P_rowind,P_colptr,L_rowind,L_colptr,LT_rowind,LT_colptr;
    VecInt perm, invperm;
    VecReal P_values,L_values,LT_values,invD;
    helper::vector<int> Parent;
};

template<class TMatrix, class TVector, class TThreadManager>
class SparseLDLSolverImpl : public sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector,TThreadManager>
{
public :
    SOFA_CLASS(SOFA_TEMPLATE3(SparseLDLSolverImpl,TMatrix,TVector,TThreadManager),SOFA_TEMPLATE3(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector,TThreadManager));
    typedef sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector, TThreadManager> Inherit;

public:
    typedef TMatrix Matrix;
    typedef TVector Vector;
    typedef TThreadManager ThreadManager;
    typedef typename TMatrix::Real Real;

protected :

    SparseLDLSolverImpl() : Inherit() {}

    template<class VecInt,class VecReal>
    void solve_cpu(Real * x,const Real * b,SpaseLDLImplInvertData<VecInt,VecReal> * data) {
        int n = data->n;
        const Real * invD = &data->invD[0];
        const int * perm = &data->perm[0];
        const int * L_colptr = &data->L_colptr[0];
        const int * L_rowind = &data->L_rowind[0];
        const Real * L_values = &data->L_values[0];
        const int * LT_colptr = &data->LT_colptr[0];
        const int * LT_rowind = &data->LT_rowind[0];
        const Real * LT_values = &data->LT_values[0];

        Tmp.clear();
        Tmp.fastResize(n);

        for (int j = 0 ; j < n ; j++) {
            Real acc = b[perm[j]];
            for (int p = LT_colptr [j] ; p < LT_colptr[j+1] ; p++) {
                acc -= LT_values[p] * Tmp[LT_rowind[p]];
            }
            Tmp[j] = acc;
        }

        for (int j = n-1 ; j >= 0 ; j--) {
            Tmp[j] *= invD[j];

            for (int p = L_colptr[j] ; p < L_colptr[j+1] ; p++) {
                Tmp[j] -= L_values[p] * Tmp[L_rowind[p]];
            }

            x[perm[j]] = Tmp[j];
        }
    }

    template<class VecInt,class VecReal>
    void factorize(TMatrix & M,SpaseLDLImplInvertData<VecInt,VecReal> * data) {
        Mfiltered.copyNonZeros(M);
        Mfiltered.compress();

        int * M_colptr = (int *) &Mfiltered.getRowBegin()[0];
        int * M_rowind = (int *) &Mfiltered.getColsIndex()[0];
        Real * M_values = (Real *) &Mfiltered.getColsValue()[0];

        bool new_factorization_needed = need_symbolic_factorization(data->P_colptr,data->P_rowind);

        // we test if the matrix has the same struct as previous factorized matrix
        if (new_factorization_needed) {
            sout << "RECOMPUTE NEW FACTORIZATION" << sendl;
            data->n = M.colSize();
            data->P_nnz = M_colptr[data->n];

            data->perm.clear();data->perm.fastResize(data->n);
            data->invperm.clear();data->invperm.fastResize(data->n);
            data->invD.clear();data->invD.fastResize(data->n);
            data->P_colptr.clear();data->P_colptr.fastResize(data->n+1);
            data->L_colptr.clear();data->L_colptr.fastResize(data->n+1);
            data->LT_colptr.clear();data->LT_colptr.fastResize(data->n+1);
            data->P_rowind.clear();data->P_rowind.fastResize(data->P_nnz);
            data->P_values.clear();data->P_values.fastResize(data->P_nnz);

            memcpy(&data->P_colptr[0],M_colptr,(data->n+1) * sizeof(int));
            memcpy(&data->P_rowind[0],M_rowind,data->P_nnz * sizeof(int));
            memcpy(&data->P_values[0],M_values,data->P_nnz * sizeof(Real));

            //ordering function
            LDL_ordering(data->n,M_colptr,M_rowind,&data->perm[0],&data->invperm[0]);

            data->Parent.clear();
            data->Parent.resize(data->n);

            //symbolic factorization
            LDL_symbolic(data->n,M_colptr,M_rowind,&data->L_colptr[0],&data->perm[0],&data->invperm[0],&data->Parent[0]);

            data->L_nnz = data->L_colptr[data->n];

            data->L_rowind.clear();data->L_rowind.fastResize(data->L_nnz);
            data->L_values.clear();data->L_values.fastResize(data->L_nnz);
            data->LT_rowind.clear();data->LT_rowind.fastResize(data->L_nnz);
            data->LT_values.clear();data->LT_values.fastResize(data->L_nnz);
        }

        Real * D = &data->invD[0];
        int * rowind = &data->L_rowind[0];
        int * colptr = &data->L_colptr[0];
        Real * values = &data->L_values[0];
        int * tran_rowind = &data->LT_rowind[0];
        int * tran_colptr = &data->LT_colptr[0];
        Real * tran_values = &data->LT_values[0];

        //Numeric Factorization
        LDL_numeric(data->n,M_colptr,M_rowind,M_values,colptr,rowind,values,D,&data->perm[0],&data->invperm[0],&data->Parent[0]);

        //inverse the diagonal
        for (int i=0;i<data->n;i++) D[i] = 1.0/D[i];

        if (new_factorization_needed) {
            //Compute transpose in tran_colptr, tran_rowind, tran_values, tran_D
            tran_countvec.clear();
            tran_countvec.resize(data->n);

            //First we count the number of value on each row.
            for (int j=0;j<data->L_nnz;j++) tran_countvec[rowind[j]]++;

            //Now we make a scan to build tran_colptr
            tran_colptr[0] = 0;
            for (int j=0;j<data->n;j++) tran_colptr[j+1] = tran_colptr[j] + tran_countvec[j];
        }

        //we clear tran_countvec becaus we use it now to stro hown many value are written on each line
        tran_countvec.clear();
        tran_countvec.resize(data->n);

        for (int j=0;j<data->n;j++) {
          for (int i=colptr[j];i<colptr[j+1];i++) {
            int line = rowind[i];<commit_msg>r9855/sofa : Add an option to compute the bloc diagonal of the LDL factorization in a dense format<commit_after>/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, version 1.0 RC 1        *
*                (c) 2006-2011 MGH, INRIA, USTL, UJF, CNRS                    *
*                                                                             *
* This library 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 2.1 of the License, or (at     *
* your option) any later version.                                             *
*                                                                             *
* This library 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 this library; if not, write to the Free Software Foundation,     *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.          *
*******************************************************************************
*                               SOFA :: Modules                               *
*                                                                             *
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: <EMAIL>                             *
******************************************************************************/
#ifndef SOFA_COMPONENT_LINEARSOLVER_SPARSELDLSOLVERIMPL_H
#define SOFA_COMPONENT_LINEARSOLVER_SPARSELDLSOLVERIMPL_H

#include <sofa/core/behavior/LinearSolver.h>
#include <sofa/component/linearsolver/MatrixLinearSolver.h>

#ifdef SOFA_HAVE_METIS
extern "C" {
#include <metis.h>
}
#endif

namespace sofa
{

namespace component
{

namespace linearsolver
{

//defaut structure for a LDL factorization
template<class VecInt,class VecReal>
class SpaseLDLImplInvertData : public MatrixInvertData {
public :
    int n, P_nnz, L_nnz, group;
    VecInt P_rowind,P_colptr,L_rowind,L_colptr,LT_rowind,LT_colptr;
    VecInt perm, invperm;
    VecReal P_values,L_values,LT_values,invD;
    VecReal Bdiag;
    helper::vector<int> Parent;
};

template<class TMatrix, class TVector, class TThreadManager>
class SparseLDLSolverImpl : public sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector,TThreadManager>
{
public :
    SOFA_CLASS(SOFA_TEMPLATE3(SparseLDLSolverImpl,TMatrix,TVector,TThreadManager),SOFA_TEMPLATE3(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector,TThreadManager));
    typedef sofa::component::linearsolver::MatrixLinearSolver<TMatrix,TVector, TThreadManager> Inherit;

public:
    typedef TMatrix Matrix;
    typedef TVector Vector;
    typedef TThreadManager ThreadManager;
    typedef typename TMatrix::Real Real;

protected :

    SparseLDLSolverImpl() : Inherit() {}

    template<class VecInt,class VecReal>
    void solve_cpu(Real * x,const Real * b,SpaseLDLImplInvertData<VecInt,VecReal> * data) {
        int n = data->n;
        const Real * invD = &data->invD[0];
        const int * perm = &data->perm[0];
        const int * L_colptr = &data->L_colptr[0];
        const int * L_rowind = &data->L_rowind[0];
        const Real * L_values = &data->L_values[0];
        const int * LT_colptr = &data->LT_colptr[0];
        const int * LT_rowind = &data->LT_rowind[0];
        const Real * LT_values = &data->LT_values[0];

        Tmp.clear();
        Tmp.fastResize(n);

        for (int j = 0 ; j < n ; j++) {
            Real acc = b[perm[j]];
            for (int p = LT_colptr [j] ; p < LT_colptr[j+1] ; p++) {
                acc -= LT_values[p] * Tmp[LT_rowind[p]];
            }
            Tmp[j] = acc;
        }

        for (int j = n-1 ; j >= 0 ; j--) {
            Tmp[j] *= invD[j];

            for (int p = L_colptr[j] ; p < L_colptr[j+1] ; p++) {
                Tmp[j] -= L_values[p] * Tmp[L_rowind[p]];
            }

            x[perm[j]] = Tmp[j];
        }
    }

    template<class VecInt,class VecReal>
    void factorize(TMatrix & M,SpaseLDLImplInvertData<VecInt,VecReal> * data,int group = 1) {
        Mfiltered.copyNonZeros(M);
        Mfiltered.compress();

        int * M_colptr = (int *) &Mfiltered.getRowBegin()[0];
        int * M_rowind = (int *) &Mfiltered.getColsIndex()[0];
        Real * M_values = (Real *) &Mfiltered.getColsValue()[0];

        bool new_factorization_needed = need_symbolic_factorization(data->P_colptr,data->P_rowind);

        // we test if the matrix has the same struct as previous factorized matrix
        if (new_factorization_needed) {
            sout << "RECOMPUTE NEW FACTORIZATION" << sendl;
            data->n = M.colSize();
            data->P_nnz = M_colptr[data->n];

            data->perm.clear();data->perm.fastResize(data->n);
            data->invperm.clear();data->invperm.fastResize(data->n);
            data->invD.clear();data->invD.fastResize(data->n);
            data->P_colptr.clear();data->P_colptr.fastResize(data->n+1);
            data->L_colptr.clear();data->L_colptr.fastResize(data->n+1);
            data->LT_colptr.clear();data->LT_colptr.fastResize(data->n+1);
            data->P_rowind.clear();data->P_rowind.fastResize(data->P_nnz);
            data->P_values.clear();data->P_values.fastResize(data->P_nnz);

            memcpy(&data->P_colptr[0],M_colptr,(data->n+1) * sizeof(int));
            memcpy(&data->P_rowind[0],M_rowind,data->P_nnz * sizeof(int));
            memcpy(&data->P_values[0],M_values,data->P_nnz * sizeof(Real));

            //ordering function
            LDL_ordering(data->n,M_colptr,M_rowind,&data->perm[0],&data->invperm[0]);

            data->Parent.clear();
            data->Parent.resize(data->n);

            //symbolic factorization
            LDL_symbolic(data->n,M_colptr,M_rowind,&data->L_colptr[0],&data->perm[0],&data->invperm[0],&data->Parent[0]);

            data->L_nnz = data->L_colptr[data->n];

            data->L_rowind.clear();data->L_rowind.fastResize(data->L_nnz);
            data->L_values.clear();data->L_values.fastResize(data->L_nnz);
            data->LT_rowind.clear();data->LT_rowind.fastResize(data->L_nnz);
            data->LT_values.clear();data->LT_values.fastResize(data->L_nnz);
        }

        Real * D = &data->invD[0];
        int * rowind = &data->L_rowind[0];
        int * colptr = &data->L_colptr[0];
        Real * values = &data->L_values[0];
        int * tran_rowind = &data->LT_rowind[0];
        int * tran_colptr = &data->LT_colptr[0];
        Real * tran_values = &data->LT_values[0];

        //Numeric Factorization
        LDL_numeric(data->n,M_colptr,M_rowind,M_values,colptr,rowind,values,D,&data->perm[0],&data->invperm[0],&data->Parent[0]);

        //inverse the diagonal
        for (int i=0;i<data->n;i++) D[i] = 1.0/D[i];

        // split the bloc diag in data->Bdiag

        data->group=group;
        if (group>1) {
            data->Bdiag.clear();
            data->Bdiag.resize(data->n * group);
            Real * diag = &data->Bdiag[0];

//            int ptr = 0;
            int begin = colptr[0];

            for (int j=0;j<data->n;j++) {
                int bn = (j/group) * group;

                for (int i=begin;i<colptr[j+1];i++) {
                    if (rowind[i]-bn<group) {
                        diag[j*group + rowind[i]-bn] = values[i];
                        diag[rowind[i]*group + j-bn] = values[i];
                    } else {
//                        rowind[ptr] = rowind[i];
//                        values[ptr] = values[i];
//                        ptr++;
                    }
                }

                begin = colptr[j+1];
//                colptr[j+1] = ptr;
            }
        }

        if (new_factorization_needed) {
            //Compute transpose in tran_colptr, tran_rowind, tran_values, tran_D
            tran_countvec.clear();
            tran_countvec.resize(data->n);

            //First we count the number of value on each row.
            for (int j=0;j<data->L_nnz;j++) tran_countvec[rowind[j]]++;

            //Now we make a scan to build tran_colptr
            tran_colptr[0] = 0;
            for (int j=0;j<data->n;j++) tran_colptr[j+1] = tran_colptr[j] + tran_countvec[j];
        }

        //we clear tran_countvec becaus we use it now to stro hown many value are written on each line
        tran_countvec.clear();
        tran_countvec.resize(data->n);

        for (int j=0;j<data->n;j++) {
          for (int i=colptr[j];i<colptr[j+1];i++) {
            int line = rowind[i];
{'fraction_non_alphanumeric': 0.1527632344386271, 'fraction_numerical': 0.009482257126236184, 'mean_word_length': 2.1560492013952635, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 213, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2978400', 'n_tokens_mistral': 5488, 'n_tokens_neox': 5259, 'n_words': 1059}