system_interface_ldl.h

Go to the documentation of this file.

extern "C"{
#include <ldl.h>
}
#include <vector>
#include <map>
#include <algorithm>
#include <sparse_matrix.h>

class SystemLDL:SparseMatrix<double>{
private:
        typedef std::map<IndexType,int> mapType;
        std::vector<double> _B;
        std::vector<double> _X;
        
        std::vector<double> Lx,D,Y ;
    std::vector<int> Li,Lp,Parent,Lnz,Flag,Pattern;
        
        mapType _M;
        
public:

void Initalize(int dimension)
{
  _dimension=dimension;
  _Ap.resize(_dimension+1);
  _M.clear();
  _B.resize(_dimension);
  _X.resize(_dimension);
}

double &A(int row,int col)
{
        
        IndexType I=IndexType(row,col);
        mapType::const_iterator ci=_M.find(I);
        if (ci==_M.end())
        {
                std::swap<int>(I.first,I.second);
                ci=_M.find(I);
        }
        assert(ci!=_M.end());
        int index=(*ci).second;

        return(_Ax[index]);
}

double &B(int i)
{return(_B[i]);}

double &X(int i)
{return (_X[i]);}


void Solve()
{
        int d,i;

  //   /* factorize A into LDL' (P and Pinv not used) */
    ldl_symbolic (_dimension, &(*_Ap.begin()), &(*_Ai.begin()), &(*Lp.begin()), 
                &(*Parent.begin()), &(*Lnz.begin()), &(*Flag.begin()), NULL, NULL) ;

  //  printf ("Nonzeros in L, excluding diagonal: %d\n", Lp [_dimension]) ;

        d = ldl_numeric (_dimension, &(*_Ap.begin()), &(*_Ai.begin()), &(*_Ax.begin()), &(*Lp.begin())
                , &(*Parent.begin()), &(*Lnz.begin()), &(*Li.begin()), &(*Lx.begin()), 
                &(*D.begin()), &(*Y.begin()), &(*Pattern.begin()),
        &(*Flag.begin()), NULL, NULL) ;

    if (d == _dimension)
    {
        /* solve Ax=b, overwriting b with the solution x */
        ldl_lsolve (_dimension, &(*_B.begin()), &(*Lp.begin()), &(*Li.begin()), &(*Lx.begin())) ;
        ldl_dsolve (_dimension, &(*_B.begin()), (&*D.begin()) );
        ldl_ltsolve (_dimension, &(*_B.begin()), &(*Lp.begin()), &(*Li.begin()), &(*Lx.begin())) ;

        for (i = 0 ; i < _dimension ; i++) 
                _X[i]=_B[i];//printf ("x [%d] = %g\n", i, b [i]) ;
    }
         else
    {
                assert(0);
    }
}

bool IsSymmetric()
{return true;}

void Zero()
{
        for (int i=0;i<Size();i++)
                _Ax[i]=0;
}

int Size(){return (_dimension);}

void CreateSparse(std::vector<IndexType> Entries)
{
        _Ax.clear();
        _Ai.clear();
        int _nonzero=0;

        std::vector<IndexType>::iterator Vi;
        for (Vi=Entries.begin();Vi<Entries.end();Vi++)
        {
                assert((*Vi).first>=0);
                assert((*Vi).second>=0);
                if ((*Vi).first>(*Vi).second)
                        std::swap<int>((*Vi).first,(*Vi).second);
        }
        
        std::sort(Entries.begin(),Entries.end());
        std::vector<IndexType>::iterator Vend=std::unique(Entries.begin(),Entries.end());
        Entries.erase(Vend,Entries.end());
        
        _Ax.resize(Entries.size());
        _Ai.resize(Entries.size());
        _M.clear();

        int col=0;
        int i=0;
        Vi=Entries.begin();
        while (Vi<Entries.end())
        {
                col=(*Vi).first;
                _Ap[i]=_nonzero;
                //go to next colummn
                while ((col==(*Vi).first)&&(Vi<Entries.end()))
                {
                        IndexType I=IndexType((*Vi).first,(*Vi).second);
                        _M.insert(std::pair<IndexType,int>(I,_nonzero));
                        _Ai[_nonzero]=(*Vi).second;
                        _nonzero++;
                        Vi++;
                }
                i++;
        }

        _Ap[_dimension]=_nonzero;
        Lx.resize(_nonzero);
        D.resize(_dimension);
        Y.resize(_dimension);
    Li.resize(_nonzero);
        Lp.resize(_dimension+1);
        Parent.resize(_dimension);
        Lnz.resize(_dimension);
        Flag.resize(_dimension);
        Pattern.resize(_dimension);
}

};