Trk::AlSpaMat Class Reference

contains the implementation for handling sparse matrices More...

#include <AlSpaMat.h>

Inheritance diagram for Trk::AlSpaMat:

Collaboration diagram for Trk::AlSpaMat:

Public Member Functions
	AlSpaMat ()
	AlSpaMat (long int N)
	AlSpaMat (const AlSpaMat &m)
	AlSpaMat (const AlSymMat &m)
	~AlSpaMat ()
AlSpaMat &	operator= (const double &)
AlSpaMat &	operator= (const AlSpaMat &)
AlSpaMat &	operator= (const AlSymMat &)
AlSpaMat &	operator= (const AlMat &)
AlSpaMat	operator+ (const AlSpaMat &) const
AlSpaMat &	operator+= (const AlSpaMat &)
AlSpaMat	operator- (const AlSpaMat &) const
AlSpaMat &	operator-= (const AlSpaMat &)
AlMat	operator* (const AlSymMatBase &) const
AlVec	operator* (const AlVec &) const
AlMat	operator* (const AlMat &) const
AlSpaMat	operator* (const double &) const
AlSpaMat &	operator*= (const double &)
int	SolveWithEigen (AlVec &RHS)
virtual void	RemoveModule (int) override final
virtual void	RemoveAlignPar (int, int) override final
void	RemoveDoF (int, int nelem=1)
virtual int	RemoveCollsRows (std::vector< int >) override final
virtual void	reSize (long int) override final
virtual void	SetPathBin (const std::string &) override final
virtual void	SetPathTxt (const std::string &) override final
virtual StatusCode	Write (const std::string &, bool, bool, double, float) override final
virtual StatusCode	CheckMatVersion (const std::string &, bool &) override final
virtual StatusCode	Read (const std::string &, int &, bool &, float &) override final
virtual StatusCode	ReadProjected (const std::string &, int &, bool &, float &) override final
virtual int	invert () override final
virtual int	diagonalize (char jobz, AlVec &w, AlMat &z) override final
virtual double &	elemr (long int, long int) override final
virtual double	elemc (long int, long int) const override final
long int	nele ()
virtual double	determinant () override final
void	elem (const indices &, long int &, long int &) const
const std::string &	pathBin () const
const std::string &	pathTxt () const
virtual TMatrixDSparse *	makeTMatrix () override final
AlSymMatBase_row	operator[] (long int)
AlSymMatBase_row_const	operator[] (long int) const
long int	nrow () const
long int	ncol () const
long int	size () const
int	matrix_type () const
const datamap *	ptrMap () const

Protected Member Functions
void	copy (const AlSpaMat &m)
void	copy (const AlSymMat &m)
void	copy (const AlMat &m)

Static Protected Member Functions
static indices	elem (long int, long int)

Protected Attributes
std::string	m_pathbin
std::string	m_pathtxt
int *	m_ptr_row = nullptr
int *	m_ptr_col = nullptr
int	m_matrix_type = 0
datamap	m_ptr_map
long int	m_size = 0
long int	m_nele = 0

Detailed Description

contains the implementation for handling sparse matrices

Definition at line 27 of file AlSpaMat.h.

Constructor & Destructor Documentation

AlSpaMat() [1/4]

Trk::AlSpaMat::AlSpaMat

(

)

Definition at line 43 of file AlSpaMat.cxx.

  : m_ptr_row(nullptr)
  , m_ptr_col(nullptr)
{
  m_matrix_type = 2;
  m_size = 0;
  m_nele = 0;
}

AlSpaMat() [2/4]

Trk::AlSpaMat::AlSpaMat

(

long int

N

)

Definition at line 53 of file AlSpaMat.cxx.

  : m_ptr_row(nullptr)
  , m_ptr_col(nullptr)
{
  m_matrix_type = 2;
  m_size = N;
  m_nele = 0;
}

AlSpaMat() [3/4]

Trk::AlSpaMat::AlSpaMat

(

const AlSpaMat &

m

)

Definition at line 63 of file AlSpaMat.cxx.

  : AlSymMatBase(m)
  , m_ptr_row(nullptr)
  , m_ptr_col(nullptr)
{
  m_matrix_type = 2;
  m_size = m.size();
  copy(m);
}

AlSpaMat() [4/4]

Trk::AlSpaMat::AlSpaMat

(

const AlSymMat &

m

)

Definition at line 74 of file AlSpaMat.cxx.

  : m_ptr_row(nullptr)
  , m_ptr_col(nullptr)
{
  m_matrix_type = 2;
  m_size = m.size();
  copy(m);
}

~AlSpaMat()

Trk::AlSpaMat::~AlSpaMat

(

)

Definition at line 84 of file AlSpaMat.cxx.

{
  m_ptr_map.clear();
  if( m_ptr_row != nullptr )  delete [] m_ptr_row;
  if( m_ptr_col != nullptr )  delete [] m_ptr_col;
}

Member Function Documentation

CheckMatVersion()

StatusCode Trk::AlSpaMat::CheckMatVersion ( const std::string & filename, bool & StdUnits )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 712 of file AlSpaMat.cxx.

{
  std::ifstream inmat((filename).c_str(), std::ios::binary);
  if(inmat.fail())
    return StatusCode::FAILURE;
 
  m_ptr_map.clear();
 
  int32_t msiz=0;
  inmat.read(charAddress(msiz), sizeof (msiz));
 
  float version=0.0;
  inmat.read(charAddress(version), sizeof (version));
 
  StdUnits = version>=2.0;
 
  inmat.close();
 
  return StatusCode::SUCCESS;
}

copy() [1/3]

void Trk::AlSpaMat::copy ( const AlMat & m )

protected

Definition at line 122 of file AlSpaMat.cxx.

{
  if( size() != m.nrow() || size() != m.ncol() ) {
    throw std::range_error( "AlSpaMat::copy: size do not match!" );
  }
 
  // copy just the lower triangle:
  m_ptr_map.clear();
  m_nele=0;
 
  //Convert Storage System
  for (int i = 0; i<m_size; i++)
    for (int j = 0; j<=i; j++)
      if  (m.elemc(i,j) != 0.) {
        m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
        m_nele++;
      }
 
  }

copy() [2/3]

void Trk::AlSpaMat::copy ( const AlSpaMat & m )

protected

Definition at line 92 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error( "AlSpaMat::copy: size do not match!" );
  }
  m_ptr_map.clear();
  m_nele=m.m_nele;
  m_ptr_map = m.m_ptr_map;
}

copy() [3/3]

void Trk::AlSpaMat::copy ( const AlSymMat & m )

protected

Definition at line 103 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error( "AlSpaMat::copy: size do not match!" );
  }
  m_ptr_map.clear();
  m_nele=0;
 
  // Convert Storage System
  for (int i = 0; i<m_size; i++)
    for (int j = 0; j<=i; j++)
      if  (m.elemc(i,j) != 0.) {
        m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
        m_nele++;
      }
 
  }

determinant()

double Trk::AlSpaMat::determinant ( )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 403 of file AlSpaMat.cxx.

{
  throw std::invalid_argument( "AlSpaMat::determinant: not implemented!" );
}

diagonalize()

int Trk::AlSpaMat::diagonalize ( char jobz, AlVec & w, AlMat & z )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 481 of file AlSpaMat.cxx.

{
  throw std::invalid_argument( "AlSpaMat::diagonalize: not implemented!" );
}

elem() [1/2]

void Trk::AlSpaMat::elem ( const indices & key, long int & i, long int & j ) const

inline

Definition at line 94 of file AlSpaMat.h.

                                                                             {
  i = key.first;
  j = key.second;
  if(j>i) throw std::out_of_range( "AlSpaMat::elem: i<j ! " );
  return;
}

elem() [2/2]

indices Trk::AlSpaMat::elem ( long int i, long int j )

staticprotected

Definition at line 209 of file AlSpaMat.cxx.

{
  // ATTENTION! the key value is returned:
#ifdef _DEBUG
  if( i<0 ) {
    throw std::out_of_range( "AlSpaMat::elem: Index 1 < zero! " );
  }
  if( i>=size() ) {
    throw std::out_of_range( "AlSpaMat::elem: Index 1 too large! " );
  }
  if( j<0 ) {
    throw std::out_of_range( "AlSpaMat::elem: Index 2 < zero! " );
  }
  if( j>=size() ) {
    throw std::out_of_range( "AlSpaMat::elem: Index 2 too large! " );
  }
#endif
 
  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
  return key;
}

elemc()

double Trk::AlSpaMat::elemc ( long int i, long int j ) const

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 179 of file AlSpaMat.cxx.

{
#ifdef _DEBUG
  if( i<0 ) {
    throw std::out_of_range( "AlSpaMat::elemc: Index 1 < zero! " );
  }
  if( i>=size() ) {
    throw std::out_of_range( "AlSpaMat::elemc: Index 1 too large! " );
  }
  if( j<0 ) {
    throw std::out_of_range( "AlSpaMat::elemc: Index 2 < zero! " );
  }
  if( j>=size() ) {
    throw std::out_of_range( "AlSpaMat::elemc: Index 2 too large! " );
  }
#endif
  // try fast referencing:
  const_mapiterator pos;
  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
  pos=m_ptr_map.find(key);
 
  if( pos!=m_ptr_map.end() )    // already exists
    return pos->second;
 
  // does not yet exist
  // return zero:
  return 0.;
}

elemr()

double & Trk::AlSpaMat::elemr ( long int i, long int j )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 143 of file AlSpaMat.cxx.

{
#ifdef _DEBUG
  if( i<0 ) {
    throw std::out_of_range( "AlSpaMat::elemr: Index 1 < zero! " );
  }
  if( i>=size() ) {
    throw std::out_of_range( "AlSpaMat::elemr: Index 1 too large! " );
  }
  if( j<0 ) {
    throw std::out_of_range( "AlSpaMat::elemr: Index 2 < zero! " );
  }
  if( j>=size() ) {
    throw std::out_of_range("AlSpaMat::elemr: Index 2 too large! " );;
  }
#endif
  // try fast referencing:
  mapiterator pos;
  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
  pos=m_ptr_map.find(key);
  if( pos!=m_ptr_map.end() )    // already exists
    return pos->second;
  else {                      // does not yet exist
    // create it with the null value:
    m_nele++;
    return (m_ptr_map.insert(std::make_pair(key, 0.))).first->second;
  }
  // simpler implementation of the above:
  /*
    m_ptr_map[key]=0;
    m_nele = m_ptr_map.size();
    return m_ptr_map[key];
  */
}

invert()

int Trk::AlSpaMat::invert ( )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 487 of file AlSpaMat.cxx.

{
  throw std::invalid_argument( "AlSpaMat::invert: not implemented!" );
}

makeTMatrix()

TMatrixDSparse * Trk::AlSpaMat::makeTMatrix ( )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 861 of file AlSpaMat.cxx.

{
 
  long int nonZeroElements = m_ptr_map.size();
  int    *irow = new int[nonZeroElements];
  int    *icol = new int[nonZeroElements];
  double *val  = new double[nonZeroElements];
 
  long int      i, j;
  long int counter(0);
  for (const datamap::value_type& p : m_ptr_map) {
    i = p.first.first;
    j = p.first.second;
    *(val+counter)=p.second;
    *(irow+counter)= i;
    *(icol+counter)= j;
     counter++;
  }
 
  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,m_size-1,0,m_size-1);
  myTMatrix->SetMatrixArray(nonZeroElements,irow,icol,val);
 
  delete [] irow;
  delete [] icol;
  delete [] val;
  return myTMatrix;
}

matrix_type()

int Trk::AlSymMatBase::matrix_type ( ) const

inlineinherited

Definition at line 148 of file AlSymMatBase.h.

{ return m_matrix_type; }

ncol()

long int Trk::AlSymMatBase::ncol ( ) const

inlineinherited

Definition at line 147 of file AlSymMatBase.h.

{ return m_size; }

nele()

long int Trk::AlSpaMat::nele ( )

inline

Definition at line 101 of file AlSpaMat.h.

                               {
  // first check for intrinsic consistency:
  if( int(m_ptr_map.size()) != m_nele ) throw std::range_error( "AlSpaMat::m_nele has been corrupted!" );
  m_nele = m_ptr_map.size();
  return m_nele;
}

nrow()

long int Trk::AlSymMatBase::nrow ( ) const

inlineinherited

Definition at line 146 of file AlSymMatBase.h.

{ return m_size; }

operator*() [1/4]

AlMat Trk::AlSpaMat::operator*

(

const AlMat &

m

)

const

Definition at line 349 of file AlSpaMat.cxx.

                                               {
  if( size() != m.nrow() ) {
    throw std::range_error(  "AlSpaMat::operator*: size do not match!" );
  }
 
  long int  isiz(size());
  long int  icol(m.ncol());
  AlMat b( isiz, icol);
 
  for( long int i=0; i<isiz; i++ )
    for( long int j=0; j<icol; j++ )
      for( long int k=0; k<isiz; k++)
        b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
 
  return b;
}

operator*() [2/4]

AlMat Trk::AlSpaMat::operator*

(

const AlSymMatBase &

m

)

const

Definition at line 331 of file AlSpaMat.cxx.

{
  if( size() != m.size() ) {
    throw std::range_error(  "AlSpaMat::operator*: size do not match!" );
  }
 
  long int  isiz(size());
  AlMat b( isiz, isiz);
 
  for( long int i=0; i<isiz; i++ )
    for( long int j=0; j<isiz; j++ )
      for( int k=0; k<isiz; k++)
        b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
 
  return b;
}

operator*() [3/4]

AlVec Trk::AlSpaMat::operator*

(

const AlVec &

v

)

const

Definition at line 367 of file AlSpaMat.cxx.

{
  if( size() != v.size() ) {
    throw std::range_error(  "AlSpaMat::operator*: size do not match! " );
  }
 
  long int  isiz(size());
  AlVec b(isiz);
 
  for( long int i=0; i<isiz; i++ )
    for( long int j=0; j<isiz; j++ )
      b[i] += elemc(i,j)*v[j];
 
  return b;
}

operator*() [4/4]

AlSpaMat Trk::AlSpaMat::operator*

(

const double &

d

)

const

Definition at line 393 of file AlSpaMat.cxx.

{
  AlSpaMat a(size());
  for (const datamap::value_type& p : m_ptr_map)
    a.m_ptr_map.emplace(p.first, p.second*d);
 
  return a;
}

operator*=()

AlSpaMat & Trk::AlSpaMat::operator*=

(

const double &

d

)

Definition at line 384 of file AlSpaMat.cxx.

{
  for (datamap::value_type& p : m_ptr_map)
    p.second *= d;
 
  return *this;
}

operator+()

AlSpaMat Trk::AlSpaMat::operator+

(

const AlSpaMat &

m

)

const

Definition at line 273 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error(  "AlSpaMat::operator+: size do not match!" );
  }
 
  AlSpaMat b(m);
  for (const datamap::value_type& p : m_ptr_map)
    b.m_ptr_map[p.first] += p.second;
  b.m_nele = b.m_ptr_map.size();
 
  return b;
}

operator+=()

AlSpaMat & Trk::AlSpaMat::operator+=

(

const AlSpaMat &

m

)

Definition at line 288 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error( "AlSpaMat::operator+=: size do not match!" );
  }
 
  for (const datamap::value_type& p : m.m_ptr_map)
    (*this).m_ptr_map[p.first] += p.second;
  m_nele = m_ptr_map.size();
 
  return *this;
}

operator-()

AlSpaMat Trk::AlSpaMat::operator-

(

const AlSpaMat &

m

)

const

Definition at line 302 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error(  "AlSpaMat::operator-: size do not match!" );
  }
 
  AlSpaMat b(m);
  for (const datamap::value_type& p : m_ptr_map)
    b.m_ptr_map[p.first] -= p.second;
  b.m_nele = b.m_ptr_map.size();
 
  return b;
}

operator-=()

AlSpaMat & Trk::AlSpaMat::operator-=

(

const AlSpaMat &

m

)

Definition at line 317 of file AlSpaMat.cxx.

{
  if( size() != m.size()) {
    throw std::range_error(  "AlSpaMat::operator-=: size do not match!" );
  }
 
  for (const datamap::value_type& p : m.m_ptr_map)
    (*this).m_ptr_map[p.first] -= p.second;
  m_nele = m_ptr_map.size();
 
  return *this;
}

operator=() [1/4]

AlSpaMat & Trk::AlSpaMat::operator=

(

const AlMat &

m

)

Definition at line 252 of file AlSpaMat.cxx.

{
  if( m.nrow() != m.ncol() ) {
    throw std::range_error( "AlSpaMat::=operator: allowed for square matrices only!" );
  }
 
  m_size=m.nrow();
  copy(m);
  return *this;
}

operator=() [2/4]

AlSpaMat & Trk::AlSpaMat::operator=

(

const AlSpaMat &

m

)

Definition at line 233 of file AlSpaMat.cxx.

{
  
 if ( this != &m ) {
    m_size=m.size();
    copy(m);
  }
  return *this;
}

operator=() [3/4]

AlSpaMat & Trk::AlSpaMat::operator=

(

const AlSymMat &

m

)

Definition at line 244 of file AlSpaMat.cxx.

{
  m_size=m.size();
  copy(m);
  return *this;
}

operator=() [4/4]

AlSpaMat & Trk::AlSpaMat::operator=

(

const double &

d

)

Definition at line 264 of file AlSpaMat.cxx.

{
  for (datamap::value_type& p : m_ptr_map)
    p.second = d;
 
  return *this;
}

operator[]() [1/2]

AlSymMatBase::AlSymMatBase_row Trk::AlSymMatBase::operator[] ( long int r )

inlineinherited

Definition at line 109 of file AlSymMatBase.h.

                                                                        {
  AlSymMatBase_row b(*this,r);
  return b;
}

operator[]() [2/2]

AlSymMatBase::AlSymMatBase_row_const Trk::AlSymMatBase::operator[] ( long int r ) const

inlineinherited

Definition at line 114 of file AlSymMatBase.h.

                                                                                    {
  const AlSymMatBase_row_const b(*this,r);
  return b;
}

pathBin()

const std::string & Trk::AlSpaMat::pathBin ( ) const

inline

Definition at line 108 of file AlSpaMat.h.

                                                {
  return m_pathbin;
}

pathTxt()

const std::string & Trk::AlSpaMat::pathTxt ( ) const

inline

Definition at line 112 of file AlSpaMat.h.

                                                {
  return m_pathtxt;
}

ptrMap()

const datamap * Trk::AlSymMatBase::ptrMap ( ) const

inlineinherited

Definition at line 149 of file AlSymMatBase.h.

{ return &m_ptr_map; }

Read()

StatusCode Trk::AlSpaMat::Read ( const std::string & filename, int & dofs, bool & triang, float & version )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 734 of file AlSpaMat.cxx.

{
  bool stdUnits = true;
  if (StatusCode::SUCCESS != CheckMatVersion(m_pathbin+filename, stdUnits))
    return StatusCode::FAILURE;
 
  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
  if(inmat.fail())
    return StatusCode::FAILURE;
 
  m_ptr_map.clear();
 
  int32_t msiz=0;
  int32_t nelem;
  inmat.read(charAddress(msiz), sizeof (msiz));
  if( msiz>999999 )
    dofs = msiz-1000000;
  else
    dofs = abs(msiz);
  m_size = dofs;
 
  if (stdUnits)
    inmat.read(charAddress(version), sizeof (version));
 
  double melem=0;
  int32_t i, j;
 
  if(msiz>999999) { // sparse format
    triang=false;
    inmat.read(charAddress(nelem), sizeof (nelem));
    m_nele=nelem;
    for(int k=0; k<nelem; k++) {
      inmat.read(charAddress(i), sizeof (i));
      inmat.read(charAddress(j), sizeof (j));
      inmat.read(charAddress(melem), sizeof (melem));
      m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
    }
  }
  else if(msiz>0) { // square format
    triang=false;
    for(int32_t i=0; i<msiz; i++) {
      for(int32_t j=0; j<msiz; j++) {
        inmat.read(charAddress(melem), sizeof (melem));
        if( i>=j && melem!=0. )
          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
      }
    }
  }
  else { // triangular format
    triang=true;
    msiz = (-1)*msiz;
    for( int32_t i=0; i<msiz; i++) {
      for( int32_t j=0; j<=i; j++) {
        inmat.read(charAddress(melem), sizeof (melem));
        if( melem!=0. )
          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
      }
    }
  }
 
  m_nele = m_ptr_map.size();
  inmat.close();
  return StatusCode::SUCCESS;
}

ReadProjected()

StatusCode Trk::AlSpaMat::ReadProjected ( const std::string & filename, int & dofs, bool & triang, float & version )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 800 of file AlSpaMat.cxx.

{
  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
  if(inmat.fail())
    return StatusCode::FAILURE;
 
  m_ptr_map.clear();
 
  int32_t msiz=0;
  int32_t nelem;
  inmat.read(charAddress(msiz), sizeof (msiz));
  if( msiz>999999 )
    dofs = msiz-1000000;
  else
    dofs = std::abs(msiz);
  m_size = dofs;
 
  inmat.read(charAddress(version), sizeof (version));
 
  double melem=0;
  int32_t i, j;
 
  if(msiz>999999) { // sparse format
    triang=false;
    inmat.read(charAddress(nelem), sizeof (nelem));
    m_nele=nelem;
    for(int k=0; k<nelem; k++) {
      inmat.read(charAddress(i), sizeof (i));
      inmat.read(charAddress(j), sizeof (j));
      inmat.read(charAddress(melem), sizeof (melem));
      m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
    }
  }
  else if(msiz>0) { // square format
    triang=false;
    for(int32_t i=0; i<msiz; i++) {
      for(int32_t j=0; j<msiz; j++) {
        inmat.read(charAddress(melem), sizeof (melem));
        if( i>=j && melem!=0. )
          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
      }
    }
  }
  else { // triangular format
    triang=true;
    msiz = (-1)*msiz;
    for( int32_t i=0; i<msiz; i++) {
      for( int32_t j=0; j<=i; j++) {
        inmat.read(charAddress(melem), sizeof (melem));
        if( melem!=0. )
          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
      }
    }
  }
 
  m_nele = m_ptr_map.size();
  inmat.close();
  return StatusCode::SUCCESS;
}

RemoveAlignPar()

void Trk::AlSpaMat::RemoveAlignPar ( int index, int control )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 564 of file AlSpaMat.cxx.

{
  // Dynamic shiftRow
  int shiftRow = 1;
  int counterRow = 0;
 
  // Dynamic shiftCol
  int shiftCol;
  int counterCol = 0;
 
  index = index-control;
 
  mapiterator pos;
 
  for(int row=index; row<(size()-shiftRow); row++) {
    shiftCol = 0;
    counterCol = 0;
 
    for(int col=0; col<row+1; col++) {
 
      if (col==index) {
        shiftCol = 1;
        counterCol=0;
      }
 
      pos = m_ptr_map.find(std::make_pair(row+shiftRow,col+shiftCol));
      if( pos!=m_ptr_map.end() ) {    // exists
        elemr(row,col) = pos->second;
        m_ptr_map.erase(pos);
      }
      else {
        pos = m_ptr_map.find(std::make_pair(row,col));
        if( pos!=m_ptr_map.end() )
          m_ptr_map.erase(pos);
      }
      counterCol++;
      if (counterCol==5-control) {
        counterCol=0;
        shiftCol++;
      }
    }
 
    counterRow++;
    if (counterRow==5-control) {
      counterRow=0;
      shiftRow++;
    }
  }
 
  m_nele = m_ptr_map.size();
}

RemoveCollsRows()

int Trk::AlSpaMat::RemoveCollsRows ( std::vector< int > indices )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 531 of file AlSpaMat.cxx.

{
  int n = indices.size();
  if (n==0) {
    return m_size;
  }
  if (n>m_size) {
    throw std::invalid_argument( "AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
    return m_size;
  }
 
  // first sort the list of indices descending
  // maybe not the fastest way to do that but it works
  for (int i=0;i<n-1;i++)
    for (int j=i+1; j<n; j++)
      if (indices[j] > indices[i]) {
        int itmp = indices[i];
        indices[i] = indices[j];
        indices[j] = itmp;
      }
 
  // remove rows and columns starting from largest indices
  for (int i=0;i<n;i++) {
    if (indices[i] > m_size-1) {
      throw std::out_of_range( "AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
    }
    RemoveDoF(indices[i]);
  }
 
  return m_size;
}

RemoveDoF()

void Trk::AlSpaMat::RemoveDoF	(	int	index,
		int	nelem = 1 )

Definition at line 493 of file AlSpaMat.cxx.

{
  // index here is the first alignment parameter to remove, nelem number of consecutive ones
 
  m_size-=nelem;
  long int n=index+nelem-1;
  mapiterator pos;
  mapiterator pos_obs=m_ptr_map.end();
  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); ++pos) {
 
    if( pos_obs!=m_ptr_map.end() )
      m_ptr_map.erase(pos_obs);
 
    pos_obs = m_ptr_map.end();
 
    long int k, l;
    elem(pos->first, k, l);
 
    if( (k>=index && k<=n) || (l>=index && l<=n) )
      pos_obs = pos;
    else {
      if( k > n )
        k -= nelem;
      if( l > n )
        l -= nelem;
      if( k >= index || l >= index )  {
        elemr(k, l) = pos->second;
        pos_obs = pos;
      }
    }
  }
 
  if( pos_obs!=m_ptr_map.end() )
    m_ptr_map.erase(pos_obs);
  m_nele = m_ptr_map.size();
}

RemoveModule()

void Trk::AlSpaMat::RemoveModule ( int index )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 617 of file AlSpaMat.cxx.

{
  // index here is module number within matrix (size=6*max_index)
  const int shift=6;
  const int ind=shift*index;
 
  RemoveDoF(ind, shift);
}

reSize()

void Trk::AlSpaMat::reSize ( long int n )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 627 of file AlSpaMat.cxx.

{
  // balanced binary tree!
  if( n!=m_size ) {
    // make a temporary copy:
    AlSpaMat m(*this);
    m_ptr_map.clear();
    m_size = n;
    long int i, j;
    for (const datamap::value_type& p : m.m_ptr_map) {
      m.elem(p.first, i, j);
      if( i<n && j<n )
        m_ptr_map.insert(p);
    }
  }
 
  m_nele = m_ptr_map.size();
}

SetPathBin()

void Trk::AlSpaMat::SetPathBin ( const std::string & path )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 647 of file AlSpaMat.cxx.

{
  m_pathbin.assign(path);
}

SetPathTxt()

void Trk::AlSpaMat::SetPathTxt ( const std::string & path )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 653 of file AlSpaMat.cxx.

{
  m_pathtxt.assign(path);
}

size()

long int Trk::AlSymMatBase::size ( ) const

inlineinherited

Definition at line 145 of file AlSymMatBase.h.

{ return m_size; }

SolveWithEigen()

int Trk::AlSpaMat::SolveWithEigen

(

AlVec &

RHS

)

Definition at line 408 of file AlSpaMat.cxx.

                                      {
 
  if(RHS.size() != size() ){
    throw std::range_error(  "AlSpaMat::SolveWithEigen vector size is incorrect" );
  }
  
  Eigen::VectorXd eigenBigVector( RHS.size() );
  for(int i(0); i< RHS.size(); ++i ){
    eigenBigVector[i] = RHS[i];
  }
 
  Eigen::SparseMatrix<double> eigenBigMatrix( m_size, m_size );
  
  using Triplet = Eigen::Triplet<double>;
  std::vector<Triplet> tripletList;
  tripletList.reserve(m_nele);
  long int      i, j;
  for (const datamap::value_type& p : m_ptr_map) {
    elem(p.first, i, j);
    tripletList.emplace_back(i,j,p.second);
    if(i!=j) tripletList.emplace_back(j,i,p.second);
  }
  eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
 
  // Eigen::CholmodSupernodalLLT  is much much quicker (x50) so it would be great to move to that in the future
  // requires an external package SuiteSparse
  // BiCGSTAB was the fastest iterative solver in Eigen from a quick test
  // SimplicialLDLT was the fastest direct solver in Eigen (x2 slower than BiCGSTAB ) 
 
  // Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
 
  solver.compute(eigenBigMatrix);
  if(solver.info()!=Eigen::Success) {
    // decomposition failed
    throw std::domain_error("AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
    return 1;
  }
 
  Eigen::VectorXd x = solver.solve( eigenBigVector );
  if(solver.info()!=Eigen::Success) {
  // solving failed
    throw std::domain_error("AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
    return 2;
  }
  
  //Copy results into vector
  for(int i(0); i< RHS.size(); ++i ){
    RHS[i] = x[i];
  }
 
  //Check accuracy
  Eigen::VectorXd residual = eigenBigMatrix * x - eigenBigVector;
  double sumresidual = 0;
  for( int i=0; i<residual.size(); ++i){
    sumresidual += fabs(residual[i]);
  }
  sumresidual /= (double) residual.size();
  
  if( sumresidual > 1e-3 ){
    throw std::overflow_error( "AlSpaMat::SolveWithEigen: your solution is no good! ");
    return 3;
  }
 
  return 0;
 
 
}

Write()

StatusCode Trk::AlSpaMat::Write ( const std::string & filename, bool binary, bool , double , float version )

finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 661 of file AlSpaMat.cxx.

{
  std::ofstream outmat;
 
  int32_t msizz = 1000000+size();
  int32_t nelem = m_ptr_map.size();
 
  if(binary) {
    outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
    if(outmat.fail())
      return StatusCode::FAILURE;
    // add 10^6 to the size to distinguish
    // from dense format
    outmat.write(charAddress(msizz), sizeof (msizz));
    outmat.write(charAddress(version), sizeof (version));
    outmat.write(charAddress(nelem), sizeof (nelem));
  }
  else {
    outmat.open((m_pathtxt+filename).c_str());
    if(outmat.fail())
      return StatusCode::FAILURE;
    outmat.setf(std::ios::fixed);
    outmat.setf(std::ios::showpoint);
    outmat.precision(3);
    outmat << "msizz: " << std::setw(6) << msizz << std::endl;
    outmat << "nelem: " << std::setw(6) << nelem << std::endl;
    outmat << "AlSpaMat version: " << std::setw(6) << version << std::endl;
  }
 
  double melem=0;
  long int ii, jj;
  int32_t i, j;
 
  for (const datamap::value_type& p : m_ptr_map) {
    melem = p.second;
    elem(p.first, ii, jj);   i=ii;  j=jj;     // just a type conversion
    if(binary) {
      outmat.write(charAddress((i)), sizeof (i));
      outmat.write(charAddress((j)), sizeof (j));
      outmat.write(charAddress((melem)), sizeof (melem));
    }
    else
      outmat << std::setw(6) << i  << std::setw(6) << j << std::setw(18) << melem << std::endl;
  }
  outmat.close();
 
  return StatusCode::SUCCESS;
}

Member Data Documentation

m_matrix_type

int Trk::AlSymMatBase::m_matrix_type = 0

protectedinherited

Definition at line 100 of file AlSymMatBase.h.

m_nele

long int Trk::AlSymMatBase::m_nele = 0

protectedinherited

Definition at line 104 of file AlSymMatBase.h.

m_pathbin

std::string Trk::AlSpaMat::m_pathbin

protected

Definition at line 85 of file AlSpaMat.h.

m_pathtxt

std::string Trk::AlSpaMat::m_pathtxt

protected

Definition at line 86 of file AlSpaMat.h.

m_ptr_col

int* Trk::AlSpaMat::m_ptr_col = nullptr

protected

Definition at line 89 of file AlSpaMat.h.

m_ptr_map

datamap Trk::AlSymMatBase::m_ptr_map

protectedinherited

Definition at line 101 of file AlSymMatBase.h.

m_ptr_row

int* Trk::AlSpaMat::m_ptr_row = nullptr

protected

Definition at line 88 of file AlSpaMat.h.

m_size

long int Trk::AlSymMatBase::m_size = 0

protectedinherited

Definition at line 103 of file AlSymMatBase.h.

---

The documentation for this class was generated from the following files:

• AlSpaMat.h
• AlSpaMat.cxx