ATLAS Offline Software
Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
Trk::AlSpaMat Class Reference

#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 ()
 
AlSpaMatoperator= (const double &)
 
AlSpaMatoperator= (const AlSpaMat &)
 
AlSpaMatoperator= (const AlSymMat &)
 
AlSpaMatoperator= (const AlMat &)
 
AlSpaMat operator+ (const AlSpaMat &) const
 
AlSpaMatoperator+= (const AlSpaMat &)
 
AlSpaMat operator- (const AlSpaMat &) const
 
AlSpaMatoperator-= (const AlSpaMat &)
 
AlMat operator* (const AlSymMatBase &) const
 
AlVec operator* (const AlVec &) const
 
AlMat operator* (const AlMat &) const
 
AlSpaMat operator* (const double &) const
 
AlSpaMatoperator*= (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 TMatrixDSparsemakeTMatrix () 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 datamapptrMap () 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 36 of file AlSpaMat.cxx.

37  : m_ptr_row(nullptr)
38  , m_ptr_col(nullptr)
39 {
40  m_matrix_type = 2;
41  m_size = 0;
42  m_nele = 0;
43 }

◆ AlSpaMat() [2/4]

Trk::AlSpaMat::AlSpaMat ( long int  N)

Definition at line 46 of file AlSpaMat.cxx.

47  : m_ptr_row(nullptr)
48  , m_ptr_col(nullptr)
49 {
50  m_matrix_type = 2;
51  m_size = N;
52  m_nele = 0;
53 }

◆ AlSpaMat() [3/4]

Trk::AlSpaMat::AlSpaMat ( const AlSpaMat m)

Definition at line 56 of file AlSpaMat.cxx.

57  : AlSymMatBase(m)
58  , m_ptr_row(nullptr)
59  , m_ptr_col(nullptr)
60 {
61  m_matrix_type = 2;
62  m_size = m.size();
63  copy(m);
64 }

◆ AlSpaMat() [4/4]

Trk::AlSpaMat::AlSpaMat ( const AlSymMat m)

Definition at line 67 of file AlSpaMat.cxx.

68  : m_ptr_row(nullptr)
69  , m_ptr_col(nullptr)
70 {
71  m_matrix_type = 2;
72  m_size = m.size();
73  copy(m);
74 }

◆ ~AlSpaMat()

Trk::AlSpaMat::~AlSpaMat ( )

Definition at line 77 of file AlSpaMat.cxx.

78 {
79  m_ptr_map.clear();
80  if( m_ptr_row != nullptr ) delete [] m_ptr_row;
81  if( m_ptr_col != nullptr ) delete [] m_ptr_col;
82 }

Member Function Documentation

◆ CheckMatVersion()

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

Implements Trk::AlSymMatBase.

Definition at line 705 of file AlSpaMat.cxx.

706 {
707  std::ifstream inmat((filename).c_str(), std::ios::binary);
708  if(inmat.fail())
709  return StatusCode::FAILURE;
710 
711  m_ptr_map.clear();
712 
713  int32_t msiz=0;
714  inmat.read((char*)&msiz, sizeof (msiz));
715 
716  float version=0.0;
717  inmat.read((char*)&version, sizeof (version));
718 
719  StdUnits = version>=2.0;
720 
721  inmat.close();
722 
723  return StatusCode::SUCCESS;
724 }

◆ copy() [1/3]

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

Definition at line 115 of file AlSpaMat.cxx.

116 {
117  if( size() != m.nrow() || size() != m.ncol() ) {
118  throw std::range_error( "AlSpaMat::copy: size do not match!" );
119  }
120 
121  // copy just the lower triangle:
122  m_ptr_map.clear();
123  m_nele=0;
124 
125  //Convert Storage System
126  for (int i = 0; i<m_size; i++)
127  for (int j = 0; j<=i; j++)
128  if (m.elemc(i,j) != 0.) {
129  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
130  m_nele++;
131  }
132 
133  }

◆ copy() [2/3]

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

Definition at line 85 of file AlSpaMat.cxx.

86 {
87  if( size() != m.size()) {
88  throw std::range_error( "AlSpaMat::copy: size do not match!" );
89  }
90  m_ptr_map.clear();
91  m_nele=m.m_nele;
92  m_ptr_map = m.m_ptr_map;
93 }

◆ copy() [3/3]

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

Definition at line 96 of file AlSpaMat.cxx.

97 {
98  if( size() != m.size()) {
99  throw std::range_error( "AlSpaMat::copy: size do not match!" );
100  }
101  m_ptr_map.clear();
102  m_nele=0;
103 
104  // Convert Storage System
105  for (int i = 0; i<m_size; i++)
106  for (int j = 0; j<=i; j++)
107  if (m.elemc(i,j) != 0.) {
108  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
109  m_nele++;
110  }
111 
112  }

◆ determinant()

double Trk::AlSpaMat::determinant ( )
finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 396 of file AlSpaMat.cxx.

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

◆ diagonalize()

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

Implements Trk::AlSymMatBase.

Definition at line 474 of file AlSpaMat.cxx.

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

◆ 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.

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

◆ elem() [2/2]

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

Definition at line 202 of file AlSpaMat.cxx.

203 {
204  // ATTENTION! the key value is returned:
205 #ifdef _DEBUG
206  if( i<0 ) {
207  throw std::out_of_range( "AlSpaMat::elem: Index 1 < zero! " );
208  }
209  if( i>=size() ) {
210  throw std::out_of_range( "AlSpaMat::elem: Index 1 too large! " );
211  }
212  if( j<0 ) {
213  throw std::out_of_range( "AlSpaMat::elem: Index 2 < zero! " );
214  }
215  if( j>=size() ) {
216  throw std::out_of_range( "AlSpaMat::elem: Index 2 too large! " );
217  }
218 #endif
219 
220  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
221  return key;
222 }

◆ elemc()

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

Implements Trk::AlSymMatBase.

Definition at line 172 of file AlSpaMat.cxx.

173 {
174 #ifdef _DEBUG
175  if( i<0 ) {
176  throw std::out_of_range( "AlSpaMat::elemc: Index 1 < zero! " );
177  }
178  if( i>=size() ) {
179  throw std::out_of_range( "AlSpaMat::elemc: Index 1 too large! " );
180  }
181  if( j<0 ) {
182  throw std::out_of_range( "AlSpaMat::elemc: Index 2 < zero! " );
183  }
184  if( j>=size() ) {
185  throw std::out_of_range( "AlSpaMat::elemc: Index 2 too large! " );
186  }
187 #endif
188  // try fast referencing:
190  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
191  pos=m_ptr_map.find(key);
192 
193  if( pos!=m_ptr_map.end() ) // already exists
194  return pos->second;
195 
196  // does not yet exist
197  // return zero:
198  return 0.;
199 }

◆ elemr()

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

Implements Trk::AlSymMatBase.

Definition at line 136 of file AlSpaMat.cxx.

137 {
138 #ifdef _DEBUG
139  if( i<0 ) {
140  throw std::out_of_range( "AlSpaMat::elemr: Index 1 < zero! " );
141  }
142  if( i>=size() ) {
143  throw std::out_of_range( "AlSpaMat::elemr: Index 1 too large! " );
144  }
145  if( j<0 ) {
146  throw std::out_of_range( "AlSpaMat::elemr: Index 2 < zero! " );
147  }
148  if( j>=size() ) {
149  throw std::out_of_range("AlSpaMat::elemr: Index 2 too large! " );;
150  }
151 #endif
152  // try fast referencing:
154  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
155  pos=m_ptr_map.find(key);
156  if( pos!=m_ptr_map.end() ) // already exists
157  return pos->second;
158  else { // does not yet exist
159  // create it with the null value:
160  m_nele++;
161  return (m_ptr_map.insert(std::make_pair(key, 0.))).first->second;
162  }
163  // simpler implementation of the above:
164  /*
165  m_ptr_map[key]=0;
166  m_nele = m_ptr_map.size();
167  return m_ptr_map[key];
168  */
169 }

◆ invert()

int Trk::AlSpaMat::invert ( )
finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 480 of file AlSpaMat.cxx.

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

◆ makeTMatrix()

TMatrixDSparse * Trk::AlSpaMat::makeTMatrix ( )
finaloverridevirtual

Implements Trk::AlSymMatBase.

Definition at line 854 of file AlSpaMat.cxx.

855 {
856 
857  long int nonZeroElements = m_ptr_map.size();
858  int *irow = new int[nonZeroElements];
859  int *icol = new int[nonZeroElements];
860  double *val = new double[nonZeroElements];
861 
862  long int i, j;
863  long int counter(0);
864  for (const datamap::value_type& p : m_ptr_map) {
865  i = p.first.first;
866  j = p.first.second;
867  *(val+counter)=p.second;
868  *(irow+counter)= i;
869  *(icol+counter)= j;
870  counter++;
871  }
872 
873  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,m_size-1,0,m_size-1);
874  myTMatrix->SetMatrixArray(nonZeroElements,irow,icol,val);
875 
876  delete [] irow;
877  delete [] icol;
878  delete [] val;
879  return myTMatrix;
880 }

◆ matrix_type()

int Trk::AlSymMatBase::matrix_type ( ) const
inlineinherited

Definition at line 148 of file AlSymMatBase.h.

148 { return m_matrix_type; }

◆ ncol()

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

Definition at line 147 of file AlSymMatBase.h.

147 { return m_size; }

◆ nele()

long int Trk::AlSpaMat::nele ( )
inline

Definition at line 101 of file AlSpaMat.h.

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

◆ nrow()

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

Definition at line 146 of file AlSymMatBase.h.

146 { return m_size; }

◆ operator*() [1/4]

AlMat Trk::AlSpaMat::operator* ( const AlMat m) const

Definition at line 342 of file AlSpaMat.cxx.

342  {
343  if( size() != m.nrow() ) {
344  throw std::range_error( "AlSpaMat::operator*: size do not match!" );
345  }
346 
347  long int isiz(size());
348  long int icol(m.ncol());
349  AlMat b( isiz, icol);
350 
351  for( long int i=0; i<isiz; i++ )
352  for( long int j=0; j<icol; j++ )
353  for( long int k=0; k<isiz; k++)
354  b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
355 
356  return b;
357 }

◆ operator*() [2/4]

AlMat Trk::AlSpaMat::operator* ( const AlSymMatBase m) const

Definition at line 324 of file AlSpaMat.cxx.

325 {
326  if( size() != m.size() ) {
327  throw std::range_error( "AlSpaMat::operator*: size do not match!" );
328  }
329 
330  long int isiz(size());
331  AlMat b( isiz, isiz);
332 
333  for( long int i=0; i<isiz; i++ )
334  for( long int j=0; j<isiz; j++ )
335  for( int k=0; k<isiz; k++)
336  b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
337 
338  return b;
339 }

◆ operator*() [3/4]

AlVec Trk::AlSpaMat::operator* ( const AlVec v) const

Definition at line 360 of file AlSpaMat.cxx.

361 {
362  if( size() != v.size() ) {
363  throw std::range_error( "AlSpaMat::operator*: size do not match! " );
364  }
365 
366  long int isiz(size());
367  AlVec b(isiz);
368 
369  for( long int i=0; i<isiz; i++ )
370  for( long int j=0; j<isiz; j++ )
371  b[i] += elemc(i,j)*v[j];
372 
373  return b;
374 }

◆ operator*() [4/4]

AlSpaMat Trk::AlSpaMat::operator* ( const double &  d) const

Definition at line 386 of file AlSpaMat.cxx.

387 {
388  AlSpaMat a(size());
389  for (const datamap::value_type& p : m_ptr_map)
390  a.m_ptr_map.emplace(p.first, p.second*d);
391 
392  return a;
393 }

◆ operator*=()

AlSpaMat & Trk::AlSpaMat::operator*= ( const double &  d)

Definition at line 377 of file AlSpaMat.cxx.

378 {
379  for (datamap::value_type& p : m_ptr_map)
380  p.second *= d;
381 
382  return *this;
383 }

◆ operator+()

AlSpaMat Trk::AlSpaMat::operator+ ( const AlSpaMat m) const

Definition at line 266 of file AlSpaMat.cxx.

267 {
268  if( size() != m.size()) {
269  throw std::range_error( "AlSpaMat::operator+: size do not match!" );
270  }
271 
272  AlSpaMat b(m);
273  for (const datamap::value_type& p : m_ptr_map)
274  b.m_ptr_map[p.first] += p.second;
275  b.m_nele = b.m_ptr_map.size();
276 
277  return b;
278 }

◆ operator+=()

AlSpaMat & Trk::AlSpaMat::operator+= ( const AlSpaMat m)

Definition at line 281 of file AlSpaMat.cxx.

282 {
283  if( size() != m.size()) {
284  throw std::range_error( "AlSpaMat::operator+=: size do not match!" );
285  }
286 
287  for (const datamap::value_type& p : m.m_ptr_map)
288  (*this).m_ptr_map[p.first] += p.second;
289  m_nele = m_ptr_map.size();
290 
291  return *this;
292 }

◆ operator-()

AlSpaMat Trk::AlSpaMat::operator- ( const AlSpaMat m) const

Definition at line 295 of file AlSpaMat.cxx.

296 {
297  if( size() != m.size()) {
298  throw std::range_error( "AlSpaMat::operator-: size do not match!" );
299  }
300 
301  AlSpaMat b(m);
302  for (const datamap::value_type& p : m_ptr_map)
303  b.m_ptr_map[p.first] -= p.second;
304  b.m_nele = b.m_ptr_map.size();
305 
306  return b;
307 }

◆ operator-=()

AlSpaMat & Trk::AlSpaMat::operator-= ( const AlSpaMat m)

Definition at line 310 of file AlSpaMat.cxx.

311 {
312  if( size() != m.size()) {
313  throw std::range_error( "AlSpaMat::operator-=: size do not match!" );
314  }
315 
316  for (const datamap::value_type& p : m.m_ptr_map)
317  (*this).m_ptr_map[p.first] -= p.second;
318  m_nele = m_ptr_map.size();
319 
320  return *this;
321 }

◆ operator=() [1/4]

AlSpaMat & Trk::AlSpaMat::operator= ( const AlMat m)

Definition at line 245 of file AlSpaMat.cxx.

246 {
247  if( m.nrow() != m.ncol() ) {
248  throw std::range_error( "AlSpaMat::=operator: allowed for square matrices only!" );
249  }
250 
251  m_size=m.nrow();
252  copy(m);
253  return *this;
254 }

◆ operator=() [2/4]

AlSpaMat & Trk::AlSpaMat::operator= ( const AlSpaMat m)

Definition at line 226 of file AlSpaMat.cxx.

227 {
228 
229  if ( this != &m ) {
230  m_size=m.size();
231  copy(m);
232  }
233  return *this;
234 }

◆ operator=() [3/4]

AlSpaMat & Trk::AlSpaMat::operator= ( const AlSymMat m)

Definition at line 237 of file AlSpaMat.cxx.

238 {
239  m_size=m.size();
240  copy(m);
241  return *this;
242 }

◆ operator=() [4/4]

AlSpaMat & Trk::AlSpaMat::operator= ( const double &  d)

Definition at line 257 of file AlSpaMat.cxx.

258 {
259  for (datamap::value_type& p : m_ptr_map)
260  p.second = d;
261 
262  return *this;
263 }

◆ operator[]() [1/2]

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

Definition at line 109 of file AlSymMatBase.h.

109  {
110  AlSymMatBase_row b(*this,r);
111  return b;
112 }

◆ operator[]() [2/2]

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

Definition at line 114 of file AlSymMatBase.h.

114  {
115  const AlSymMatBase_row_const b(*this,r);
116  return b;
117 }

◆ pathBin()

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

Definition at line 108 of file AlSpaMat.h.

108  {
109  return m_pathbin;
110 }

◆ pathTxt()

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

Definition at line 112 of file AlSpaMat.h.

112  {
113  return m_pathtxt;
114 }

◆ ptrMap()

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

Definition at line 149 of file AlSymMatBase.h.

149 { 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 727 of file AlSpaMat.cxx.

728 {
729  bool stdUnits = true;
730  if (StatusCode::SUCCESS != CheckMatVersion(m_pathbin+filename, stdUnits))
731  return StatusCode::FAILURE;
732 
733  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
734  if(inmat.fail())
735  return StatusCode::FAILURE;
736 
737  m_ptr_map.clear();
738 
739  int32_t msiz=0;
740  int32_t nelem;
741  inmat.read((char*)&msiz, sizeof (msiz));
742  if( msiz>999999 )
743  dofs = msiz-1000000;
744  else
745  dofs = abs(msiz);
746  m_size = dofs;
747 
748  if (stdUnits)
749  inmat.read((char*)&version, sizeof (version));
750 
751  double melem=0;
752  int32_t i, j;
753 
754  if(msiz>999999) { // sparse format
755  triang=false;
756  inmat.read((char*)&nelem, sizeof (nelem));
757  m_nele=nelem;
758  for(int k=0; k<nelem; k++) {
759  inmat.read((char*)&i, sizeof (i));
760  inmat.read((char*)&j, sizeof (j));
761  inmat.read((char*)&melem, sizeof (melem));
762  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
763  }
764  }
765  else if(msiz>0) { // square format
766  triang=false;
767  for(int32_t i=0; i<msiz; i++) {
768  for(int32_t j=0; j<msiz; j++) {
769  inmat.read((char*)&melem, sizeof (melem));
770  if( i>=j && melem!=0. )
771  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
772  }
773  }
774  }
775  else { // triangular format
776  triang=true;
777  msiz = (-1)*msiz;
778  for( int32_t i=0; i<msiz; i++) {
779  for( int32_t j=0; j<=i; j++) {
780  inmat.read((char*)&melem, sizeof (melem));
781  if( melem!=0. )
782  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
783  }
784  }
785  }
786 
787  m_nele = m_ptr_map.size();
788  inmat.close();
789  return StatusCode::SUCCESS;
790 }

◆ ReadProjected()

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

Implements Trk::AlSymMatBase.

Definition at line 793 of file AlSpaMat.cxx.

795 {
796  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
797  if(inmat.fail())
798  return StatusCode::FAILURE;
799 
800  m_ptr_map.clear();
801 
802  int32_t msiz=0;
803  int32_t nelem;
804  inmat.read((char*)&msiz, sizeof (msiz));
805  if( msiz>999999 )
806  dofs = msiz-1000000;
807  else
808  dofs = abs(msiz);
809  m_size = dofs;
810 
811  inmat.read((char*)&version, sizeof (version));
812 
813  double melem=0;
814  int32_t i, j;
815 
816  if(msiz>999999) { // sparse format
817  triang=false;
818  inmat.read((char*)&nelem, sizeof (nelem));
819  m_nele=nelem;
820  for(int k=0; k<nelem; k++) {
821  inmat.read((char*)&i, sizeof (i));
822  inmat.read((char*)&j, sizeof (j));
823  inmat.read((char*)&melem, sizeof (melem));
824  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
825  }
826  }
827  else if(msiz>0) { // square format
828  triang=false;
829  for(int32_t i=0; i<msiz; i++) {
830  for(int32_t j=0; j<msiz; j++) {
831  inmat.read((char*)&melem, sizeof (melem));
832  if( i>=j && melem!=0. )
833  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
834  }
835  }
836  }
837  else { // triangular format
838  triang=true;
839  msiz = (-1)*msiz;
840  for( int32_t i=0; i<msiz; i++) {
841  for( int32_t j=0; j<=i; j++) {
842  inmat.read((char*)&melem, sizeof (melem));
843  if( melem!=0. )
844  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
845  }
846  }
847  }
848 
849  m_nele = m_ptr_map.size();
850  inmat.close();
851  return StatusCode::SUCCESS;
852 }

◆ RemoveAlignPar()

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

Implements Trk::AlSymMatBase.

Definition at line 557 of file AlSpaMat.cxx.

558 {
559  // Dynamic shiftRow
560  int shiftRow = 1;
561  int counterRow = 0;
562 
563  // Dynamic shiftCol
564  int shiftCol;
565  int counterCol = 0;
566 
567  index = index-control;
568 
570 
571  for(int row=index; row<(size()-shiftRow); row++) {
572  shiftCol = 0;
573  counterCol = 0;
574 
575  for(int col=0; col<row+1; col++) {
576 
577  if (col==index) {
578  shiftCol = 1;
579  counterCol=0;
580  }
581 
582  pos = m_ptr_map.find(std::make_pair(row+shiftRow,col+shiftCol));
583  if( pos!=m_ptr_map.end() ) { // exists
584  elemr(row,col) = pos->second;
585  m_ptr_map.erase(pos);
586  }
587  else {
588  pos = m_ptr_map.find(std::make_pair(row,col));
589  if( pos!=m_ptr_map.end() )
590  m_ptr_map.erase(pos);
591  }
592  counterCol++;
593  if (counterCol==5-control) {
594  counterCol=0;
595  shiftCol++;
596  }
597  }
598 
599  counterRow++;
600  if (counterRow==5-control) {
601  counterRow=0;
602  shiftRow++;
603  }
604  }
605 
606  m_nele = m_ptr_map.size();
607 }

◆ RemoveCollsRows()

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

Implements Trk::AlSymMatBase.

Definition at line 524 of file AlSpaMat.cxx.

525 {
526  int n = indices.size();
527  if (n==0) {
528  return m_size;
529  }
530  if (n>m_size) {
531  throw std::invalid_argument( "AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
532  return m_size;
533  }
534 
535  // first sort the list of indices descending
536  // maybe not the fastest way to do that but it works
537  for (int i=0;i<n-1;i++)
538  for (int j=i+1; j<n; j++)
539  if (indices[j] > indices[i]) {
540  int itmp = indices[i];
541  indices[i] = indices[j];
542  indices[j] = itmp;
543  }
544 
545  // remove rows and columns starting from largest indices
546  for (int i=0;i<n;i++) {
547  if (indices[i] > m_size-1) {
548  throw std::out_of_range( "AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
549  }
550  RemoveDoF(indices[i]);
551  }
552 
553  return m_size;
554 }

◆ RemoveDoF()

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

Definition at line 486 of file AlSpaMat.cxx.

487 {
488  // index here is the first alignment parameter to remove, nelem number of consecutive ones
489 
490  m_size-=nelem;
491  long int n=index+nelem-1;
493  mapiterator pos_obs=m_ptr_map.end();
494  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); ++pos) {
495 
496  if( pos_obs!=m_ptr_map.end() )
497  m_ptr_map.erase(pos_obs);
498 
499  pos_obs = m_ptr_map.end();
500 
501  long int k, l;
502  elem(pos->first, k, l);
503 
504  if( (k>=index && k<=n) || (l>=index && l<=n) )
505  pos_obs = pos;
506  else {
507  if( k > n )
508  k -= nelem;
509  if( l > n )
510  l -= nelem;
511  if( k >= index || l >= index ) {
512  elemr(k, l) = pos->second;
513  pos_obs = pos;
514  }
515  }
516  }
517 
518  if( pos_obs!=m_ptr_map.end() )
519  m_ptr_map.erase(pos_obs);
520  m_nele = m_ptr_map.size();
521 }

◆ RemoveModule()

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

Implements Trk::AlSymMatBase.

Definition at line 610 of file AlSpaMat.cxx.

611 {
612  // index here is module number within matrix (size=6*max_index)
613  const int shift=6;
614  const int ind=shift*index;
615 
616  RemoveDoF(ind, shift);
617 }

◆ reSize()

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

Implements Trk::AlSymMatBase.

Definition at line 620 of file AlSpaMat.cxx.

621 {
622  // balanced binary tree!
623  if( n!=m_size ) {
624  // make a temporary copy:
625  AlSpaMat m(*this);
626  m_ptr_map.clear();
627  m_size = n;
628  long int i, j;
629  for (const datamap::value_type& p : m.m_ptr_map) {
630  m.elem(p.first, i, j);
631  if( i<n && j<n )
632  m_ptr_map.insert(p);
633  }
634  }
635 
636  m_nele = m_ptr_map.size();
637 }

◆ SetPathBin()

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

Implements Trk::AlSymMatBase.

Definition at line 640 of file AlSpaMat.cxx.

641 {
642  m_pathbin.assign(path);
643 }

◆ SetPathTxt()

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

Implements Trk::AlSymMatBase.

Definition at line 646 of file AlSpaMat.cxx.

647 {
648  m_pathtxt.assign(path);
649 }

◆ size()

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

Definition at line 145 of file AlSymMatBase.h.

145 { return m_size; }

◆ SolveWithEigen()

int Trk::AlSpaMat::SolveWithEigen ( AlVec RHS)

Definition at line 401 of file AlSpaMat.cxx.

401  {
402 
403  if(RHS.size() != size() ){
404  throw std::range_error( "AlSpaMat::SolveWithEigen vector size is incorrect" );
405  }
406 
407  Eigen::VectorXd eigenBigVector( RHS.size() );
408  for(int i(0); i< RHS.size(); ++i ){
409  eigenBigVector[i] = RHS[i];
410  }
411 
412  Eigen::SparseMatrix<double> eigenBigMatrix( m_size, m_size );
413 
414  using Triplet = Eigen::Triplet<double>;
415  std::vector<Triplet> tripletList;
416  tripletList.reserve(m_nele);
417  long int i, j;
418  for (const datamap::value_type& p : m_ptr_map) {
419  elem(p.first, i, j);
420  tripletList.emplace_back(i,j,p.second);
421  if(i!=j) tripletList.emplace_back(j,i,p.second);
422  }
423  eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
424 
425  // Eigen::CholmodSupernodalLLT is much much quicker (x50) so it would be great to move to that in the future
426  // requires an external package SuiteSparse
427  // BiCGSTAB was the fastest iterative solver in Eigen from a quick test
428  // SimplicialLDLT was the fastest direct solver in Eigen (x2 slower than BiCGSTAB )
429 
430  // Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
431  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
432 
433  solver.compute(eigenBigMatrix);
434  if(solver.info()!=Eigen::Success) {
435  // decomposition failed
436  throw std::domain_error("AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
437  return 1;
438  }
439 
440  Eigen::VectorXd x = solver.solve( eigenBigVector );
441  if(solver.info()!=Eigen::Success) {
442  // solving failed
443  throw std::domain_error("AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
444  return 2;
445  }
446 
447  //Copy results into vector
448  for(int i(0); i< RHS.size(); ++i ){
449  RHS[i] = x[i];
450  }
451 
452  //Check accuracy
453  Eigen::VectorXd residual = eigenBigMatrix * x - eigenBigVector;
454  double sumresidual = 0;
455  for( int i=0; i<residual.size(); ++i){
456  sumresidual += fabs(residual[i]);
457  }
458  sumresidual /= (double) residual.size();
459 
460  if( sumresidual > 1e-3 ){
461  throw std::overflow_error( "AlSpaMat::SolveWithEigen: your solution is no good! ");
462  return 3;
463  }
464 
465  return 0;
466 
467 
468 }

◆ Write()

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

Implements Trk::AlSymMatBase.

Definition at line 654 of file AlSpaMat.cxx.

656 {
657  std::ofstream outmat;
658 
659  int32_t msizz = 1000000+size();
660  int32_t nelem = m_ptr_map.size();
661 
662  if(binary) {
663  outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
664  if(outmat.fail())
665  return StatusCode::FAILURE;
666  // add 10^6 to the size to distinguish
667  // from dense format
668  outmat.write((char*)&msizz, sizeof (msizz));
669  outmat.write((char*)&version, sizeof (version));
670  outmat.write((char*)&nelem, sizeof (nelem));
671  }
672  else {
673  outmat.open((m_pathtxt+filename).c_str());
674  if(outmat.fail())
675  return StatusCode::FAILURE;
676  outmat.setf(std::ios::fixed);
677  outmat.setf(std::ios::showpoint);
678  outmat.precision(3);
679  outmat << "msizz: " << std::setw(6) << msizz << std::endl;
680  outmat << "nelem: " << std::setw(6) << nelem << std::endl;
681  outmat << "AlSpaMat version: " << std::setw(6) << version << std::endl;
682  }
683 
684  double melem=0;
685  long int ii, jj;
686  int32_t i, j;
687 
688  for (const datamap::value_type& p : m_ptr_map) {
689  melem = p.second;
690  elem(p.first, ii, jj); i=ii; j=jj; // just a type conversion
691  if(binary) {
692  outmat.write((char*)&(i), sizeof (i));
693  outmat.write((char*)&(j), sizeof (j));
694  outmat.write((char*)&(melem), sizeof (melem));
695  }
696  else
697  outmat << std::setw(6) << i << std::setw(6) << j << std::setw(18) << melem << std::endl;
698  }
699  outmat.close();
700 
701  return StatusCode::SUCCESS;
702 }

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:
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
query_example.row
row
Definition: query_example.py:24
beamspotman.r
def r
Definition: beamspotman.py:676
Trk::AlSpaMat::m_ptr_row
int * m_ptr_row
Definition: AlSpaMat.h:88
Trk::const_mapiterator
datamap::const_iterator const_mapiterator
Definition: AlSymMatBase.h:30
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
index
Definition: index.py:1
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Trk::AlSpaMat::elemr
virtual double & elemr(long int, long int) override final
Definition: AlSpaMat.cxx:136
Trk::indices
std::pair< long int, long int > indices
Definition: AlSymMatBase.h:24
Trk::AlSpaMat::copy
void copy(const AlSpaMat &m)
Definition: AlSpaMat.cxx:85
Trk::AlSymMatBase::m_ptr_map
datamap m_ptr_map
Definition: AlSymMatBase.h:101
Trk::AlSpaMat::elemc
virtual double elemc(long int, long int) const override final
Definition: AlSpaMat.cxx:172
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
Trk::AlSpaMat::m_pathbin
std::string m_pathbin
Definition: AlSpaMat.h:85
Trk::AlSpaMat::CheckMatVersion
virtual StatusCode CheckMatVersion(const std::string &, bool &) override final
Definition: AlSpaMat.cxx:705
atlasStyleMacro.icol
int icol
Definition: atlasStyleMacro.py:13
Trk::AlSpaMat::elem
void elem(const indices &, long int &, long int &) const
Definition: AlSpaMat.h:94
Trk::AlSymMatBase::AlSymMatBase
AlSymMatBase()
Definition: AlSymMatBase.cxx:20
TMatrixDSparse
class TMatrixTSparse< double > TMatrixDSparse
Definition: AlSpaMat.h:14
Trk::AlSymMatBase::size
long int size() const
Definition: AlSymMatBase.h:145
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
Trk::AlSpaMat::AlSpaMat
AlSpaMat()
Definition: AlSpaMat.cxx:36
Trk::AlSymMatBase::m_nele
long int m_nele
Definition: AlSymMatBase.h:104
Trk::AlSpaMat::RemoveDoF
void RemoveDoF(int, int nelem=1)
Definition: AlSpaMat.cxx:486
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Trk::mapiterator
datamap::iterator mapiterator
Definition: AlSymMatBase.h:29
query_example.col
col
Definition: query_example.py:7
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
get_generator_info.version
version
Definition: get_generator_info.py:33
DeMoScan.index
string index
Definition: DeMoScan.py:364
a
TList * a
Definition: liststreamerinfos.cxx:10
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
Trk::AlSpaMat::m_ptr_col
int * m_ptr_col
Definition: AlSpaMat.h:89
Trk::AlSymMatBase::m_matrix_type
int m_matrix_type
Definition: AlSymMatBase.h:100
Trk::x
@ x
Definition: ParamDefs.h:55
test_pyathena.counter
counter
Definition: test_pyathena.py:15
Trk::AlSpaMat::m_pathtxt
std::string m_pathtxt
Definition: AlSpaMat.h:86
Trk::AlSymMatBase::m_size
long int m_size
Definition: AlSymMatBase.h:103
checkFileSG.ind
list ind
Definition: checkFileSG.py:118
Trk::v
@ v
Definition: ParamDefs.h:78
fitman.k
k
Definition: fitman.py:528
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37