ATLAS Offline Software
AlSpaMat.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Original Author: Anthony Morley (15Jan2007)
6 // MKU (15Jan2007)
7 // PBdR (17Apr2007) interface ~identcal to AlSymMat
8 
9 #include "GaudiKernel/StatusCode.h"
11 
14 #include "TrkAlgebraUtils/AlVec.h"
15 #include "TrkAlgebraUtils/AlMat.h"
16 
17 #include <algorithm>
18 #include <cfloat>
19 #include <cmath>
20 #include <cstdint>
21 #include <fstream>
22 #include <iomanip>
23 #include <stdexcept>
24 
25 #include <TMatrixDSparse.h>
26 
27 
28 #include <Eigen/Core>
29 #include <Eigen/IterativeLinearSolvers>
30 #include <Eigen/SparseCholesky>
31 
32 namespace{
33  char *
34  charAddress(auto & v){
35  return reinterpret_cast<char *>(&v);
36  }
37 }
38 
39 
40 namespace Trk {
41 
42 //______________________________________________________________________________
44  : m_ptr_row(nullptr)
45  , m_ptr_col(nullptr)
46 {
47  m_matrix_type = 2;
48  m_size = 0;
49  m_nele = 0;
50 }
51 
52 //______________________________________________________________________________
54  : m_ptr_row(nullptr)
55  , m_ptr_col(nullptr)
56 {
57  m_matrix_type = 2;
58  m_size = N;
59  m_nele = 0;
60 }
61 
62 //______________________________________________________________________________
64  : AlSymMatBase(m)
65  , m_ptr_row(nullptr)
66  , m_ptr_col(nullptr)
67 {
68  m_matrix_type = 2;
69  m_size = m.size();
70  copy(m);
71 }
72 
73 //______________________________________________________________________________
75  : m_ptr_row(nullptr)
76  , m_ptr_col(nullptr)
77 {
78  m_matrix_type = 2;
79  m_size = m.size();
80  copy(m);
81 }
82 
83 //______________________________________________________________________________
85 {
86  m_ptr_map.clear();
87  if( m_ptr_row != nullptr ) delete [] m_ptr_row;
88  if( m_ptr_col != nullptr ) delete [] m_ptr_col;
89 }
90 
91 //______________________________________________________________________________
93 {
94  if( size() != m.size()) {
95  throw std::range_error( "AlSpaMat::copy: size do not match!" );
96  }
97  m_ptr_map.clear();
98  m_nele=m.m_nele;
99  m_ptr_map = m.m_ptr_map;
100 }
101 
102 //______________________________________________________________________________
104 {
105  if( size() != m.size()) {
106  throw std::range_error( "AlSpaMat::copy: size do not match!" );
107  }
108  m_ptr_map.clear();
109  m_nele=0;
110 
111  // Convert Storage System
112  for (int i = 0; i<m_size; i++)
113  for (int j = 0; j<=i; j++)
114  if (m.elemc(i,j) != 0.) {
115  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
116  m_nele++;
117  }
118 
119  }
120 
121 //______________________________________________________________________________
122 void AlSpaMat::copy(const AlMat& m)
123 {
124  if( size() != m.nrow() || size() != m.ncol() ) {
125  throw std::range_error( "AlSpaMat::copy: size do not match!" );
126  }
127 
128  // copy just the lower triangle:
129  m_ptr_map.clear();
130  m_nele=0;
131 
132  //Convert Storage System
133  for (int i = 0; i<m_size; i++)
134  for (int j = 0; j<=i; j++)
135  if (m.elemc(i,j) != 0.) {
136  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
137  m_nele++;
138  }
139 
140  }
141 
142 //______________________________________________________________________________
143 double& AlSpaMat::elemr(long int i,long int j)
144 {
145 #ifdef _DEBUG
146  if( i<0 ) {
147  throw std::out_of_range( "AlSpaMat::elemr: Index 1 < zero! " );
148  }
149  if( i>=size() ) {
150  throw std::out_of_range( "AlSpaMat::elemr: Index 1 too large! " );
151  }
152  if( j<0 ) {
153  throw std::out_of_range( "AlSpaMat::elemr: Index 2 < zero! " );
154  }
155  if( j>=size() ) {
156  throw std::out_of_range("AlSpaMat::elemr: Index 2 too large! " );;
157  }
158 #endif
159  // try fast referencing:
161  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
162  pos=m_ptr_map.find(key);
163  if( pos!=m_ptr_map.end() ) // already exists
164  return pos->second;
165  else { // does not yet exist
166  // create it with the null value:
167  m_nele++;
168  return (m_ptr_map.insert(std::make_pair(key, 0.))).first->second;
169  }
170  // simpler implementation of the above:
171  /*
172  m_ptr_map[key]=0;
173  m_nele = m_ptr_map.size();
174  return m_ptr_map[key];
175  */
176 }
177 
178 //______________________________________________________________________________
179 double AlSpaMat::elemc(long int i,long int j) const
180 {
181 #ifdef _DEBUG
182  if( i<0 ) {
183  throw std::out_of_range( "AlSpaMat::elemc: Index 1 < zero! " );
184  }
185  if( i>=size() ) {
186  throw std::out_of_range( "AlSpaMat::elemc: Index 1 too large! " );
187  }
188  if( j<0 ) {
189  throw std::out_of_range( "AlSpaMat::elemc: Index 2 < zero! " );
190  }
191  if( j>=size() ) {
192  throw std::out_of_range( "AlSpaMat::elemc: Index 2 too large! " );
193  }
194 #endif
195  // try fast referencing:
197  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
198  pos=m_ptr_map.find(key);
199 
200  if( pos!=m_ptr_map.end() ) // already exists
201  return pos->second;
202 
203  // does not yet exist
204  // return zero:
205  return 0.;
206 }
207 
208 //______________________________________________________________________________
209 indices AlSpaMat::elem(long int i,long int j)
210 {
211  // ATTENTION! the key value is returned:
212 #ifdef _DEBUG
213  if( i<0 ) {
214  throw std::out_of_range( "AlSpaMat::elem: Index 1 < zero! " );
215  }
216  if( i>=size() ) {
217  throw std::out_of_range( "AlSpaMat::elem: Index 1 too large! " );
218  }
219  if( j<0 ) {
220  throw std::out_of_range( "AlSpaMat::elem: Index 2 < zero! " );
221  }
222  if( j>=size() ) {
223  throw std::out_of_range( "AlSpaMat::elem: Index 2 too large! " );
224  }
225 #endif
226 
227  indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
228  return key;
229 }
230 
231 //______________________________________________________________________________
232 // cppcheck-suppress operatorEqVarError; m_ptr_row/col deliberately not copied.
234 {
235 
236  if ( this != &m ) {
237  m_size=m.size();
238  copy(m);
239  }
240  return *this;
241 }
242 
243 //______________________________________________________________________________
245 {
246  m_size=m.size();
247  copy(m);
248  return *this;
249 }
250 
251 //______________________________________________________________________________
253 {
254  if( m.nrow() != m.ncol() ) {
255  throw std::range_error( "AlSpaMat::=operator: allowed for square matrices only!" );
256  }
257 
258  m_size=m.nrow();
259  copy(m);
260  return *this;
261 }
262 
263 //______________________________________________________________________________
265 {
266  for (datamap::value_type& p : m_ptr_map)
267  p.second = d;
268 
269  return *this;
270 }
271 
272 //______________________________________________________________________________
274 {
275  if( size() != m.size()) {
276  throw std::range_error( "AlSpaMat::operator+: size do not match!" );
277  }
278 
279  AlSpaMat b(m);
280  for (const datamap::value_type& p : m_ptr_map)
281  b.m_ptr_map[p.first] += p.second;
282  b.m_nele = b.m_ptr_map.size();
283 
284  return b;
285 }
286 
287 //______________________________________________________________________________
289 {
290  if( size() != m.size()) {
291  throw std::range_error( "AlSpaMat::operator+=: size do not match!" );
292  }
293 
294  for (const datamap::value_type& p : m.m_ptr_map)
295  (*this).m_ptr_map[p.first] += p.second;
296  m_nele = m_ptr_map.size();
297 
298  return *this;
299 }
300 
301 //______________________________________________________________________________
303 {
304  if( size() != m.size()) {
305  throw std::range_error( "AlSpaMat::operator-: size do not match!" );
306  }
307 
308  AlSpaMat b(m);
309  for (const datamap::value_type& p : m_ptr_map)
310  b.m_ptr_map[p.first] -= p.second;
311  b.m_nele = b.m_ptr_map.size();
312 
313  return b;
314 }
315 
316 //______________________________________________________________________________
318 {
319  if( size() != m.size()) {
320  throw std::range_error( "AlSpaMat::operator-=: size do not match!" );
321  }
322 
323  for (const datamap::value_type& p : m.m_ptr_map)
324  (*this).m_ptr_map[p.first] -= p.second;
325  m_nele = m_ptr_map.size();
326 
327  return *this;
328 }
329 
330 //______________________________________________________________________________
332 {
333  if( size() != m.size() ) {
334  throw std::range_error( "AlSpaMat::operator*: size do not match!" );
335  }
336 
337  long int isiz(size());
338  AlMat b( isiz, isiz);
339 
340  for( long int i=0; i<isiz; i++ )
341  for( long int j=0; j<isiz; j++ )
342  for( int k=0; k<isiz; k++)
343  b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
344 
345  return b;
346 }
347 
348 //______________________________________________________________________________
350  if( size() != m.nrow() ) {
351  throw std::range_error( "AlSpaMat::operator*: size do not match!" );
352  }
353 
354  long int isiz(size());
355  long int icol(m.ncol());
356  AlMat b( isiz, icol);
357 
358  for( long int i=0; i<isiz; i++ )
359  for( long int j=0; j<icol; j++ )
360  for( long int k=0; k<isiz; k++)
361  b.elemr(i,j) += elemc(i,k)*m.elemc(k,j);
362 
363  return b;
364 }
365 
366 //______________________________________________________________________________
368 {
369  if( size() != v.size() ) {
370  throw std::range_error( "AlSpaMat::operator*: size do not match! " );
371  }
372 
373  long int isiz(size());
374  AlVec b(isiz);
375 
376  for( long int i=0; i<isiz; i++ )
377  for( long int j=0; j<isiz; j++ )
378  b[i] += elemc(i,j)*v[j];
379 
380  return b;
381 }
382 
383 //______________________________________________________________________________
385 {
386  for (datamap::value_type& p : m_ptr_map)
387  p.second *= d;
388 
389  return *this;
390 }
391 
392 //______________________________________________________________________________
393 AlSpaMat AlSpaMat::operator*(const double& d) const
394 {
395  AlSpaMat a(size());
396  for (const datamap::value_type& p : m_ptr_map)
397  a.m_ptr_map.emplace(p.first, p.second*d);
398 
399  return a;
400 }
401 
402 //______________________________________________________________________________
404 {
405  throw std::invalid_argument( "AlSpaMat::determinant: not implemented!" );
406 }
407 
409 
410  if(RHS.size() != size() ){
411  throw std::range_error( "AlSpaMat::SolveWithEigen vector size is incorrect" );
412  }
413 
414  Eigen::VectorXd eigenBigVector( RHS.size() );
415  for(int i(0); i< RHS.size(); ++i ){
416  eigenBigVector[i] = RHS[i];
417  }
418 
419  Eigen::SparseMatrix<double> eigenBigMatrix( m_size, m_size );
420 
421  using Triplet = Eigen::Triplet<double>;
422  std::vector<Triplet> tripletList;
423  tripletList.reserve(m_nele);
424  long int i, j;
425  for (const datamap::value_type& p : m_ptr_map) {
426  elem(p.first, i, j);
427  tripletList.emplace_back(i,j,p.second);
428  if(i!=j) tripletList.emplace_back(j,i,p.second);
429  }
430  eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
431 
432  // Eigen::CholmodSupernodalLLT is much much quicker (x50) so it would be great to move to that in the future
433  // requires an external package SuiteSparse
434  // BiCGSTAB was the fastest iterative solver in Eigen from a quick test
435  // SimplicialLDLT was the fastest direct solver in Eigen (x2 slower than BiCGSTAB )
436 
437  // Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
438  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
439 
440  solver.compute(eigenBigMatrix);
441  if(solver.info()!=Eigen::Success) {
442  // decomposition failed
443  throw std::domain_error("AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
444  return 1;
445  }
446 
447  Eigen::VectorXd x = solver.solve( eigenBigVector );
448  if(solver.info()!=Eigen::Success) {
449  // solving failed
450  throw std::domain_error("AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
451  return 2;
452  }
453 
454  //Copy results into vector
455  for(int i(0); i< RHS.size(); ++i ){
456  RHS[i] = x[i];
457  }
458 
459  //Check accuracy
460  Eigen::VectorXd residual = eigenBigMatrix * x - eigenBigVector;
461  double sumresidual = 0;
462  for( int i=0; i<residual.size(); ++i){
463  sumresidual += fabs(residual[i]);
464  }
465  sumresidual /= (double) residual.size();
466 
467  if( sumresidual > 1e-3 ){
468  throw std::overflow_error( "AlSpaMat::SolveWithEigen: your solution is no good! ");
469  return 3;
470  }
471 
472  return 0;
473 
474 
475 }
476 
477 
478 
479 //______________________________________________________________________________
480 //jlove int AlSpaMat::diagonalize(char jobz, AlVec& w, AlMat& z) {
482 {
483  throw std::invalid_argument( "AlSpaMat::diagonalize: not implemented!" );
484 }
485 
486 //______________________________________________________________________________
488 {
489  throw std::invalid_argument( "AlSpaMat::invert: not implemented!" );
490 }
491 
492 //______________________________________________________________________________
493 void AlSpaMat::RemoveDoF(int index, int nelem)
494 {
495  // index here is the first alignment parameter to remove, nelem number of consecutive ones
496 
497  m_size-=nelem;
498  long int n=index+nelem-1;
500  mapiterator pos_obs=m_ptr_map.end();
501  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); ++pos) {
502 
503  if( pos_obs!=m_ptr_map.end() )
504  m_ptr_map.erase(pos_obs);
505 
506  pos_obs = m_ptr_map.end();
507 
508  long int k, l;
509  elem(pos->first, k, l);
510 
511  if( (k>=index && k<=n) || (l>=index && l<=n) )
512  pos_obs = pos;
513  else {
514  if( k > n )
515  k -= nelem;
516  if( l > n )
517  l -= nelem;
518  if( k >= index || l >= index ) {
519  elemr(k, l) = pos->second;
520  pos_obs = pos;
521  }
522  }
523  }
524 
525  if( pos_obs!=m_ptr_map.end() )
526  m_ptr_map.erase(pos_obs);
527  m_nele = m_ptr_map.size();
528 }
529 
530 //______________________________________________________________________________
531 int AlSpaMat::RemoveCollsRows(std::vector<int> indices)
532 {
533  int n = indices.size();
534  if (n==0) {
535  return m_size;
536  }
537  if (n>m_size) {
538  throw std::invalid_argument( "AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
539  return m_size;
540  }
541 
542  // first sort the list of indices descending
543  // maybe not the fastest way to do that but it works
544  for (int i=0;i<n-1;i++)
545  for (int j=i+1; j<n; j++)
546  if (indices[j] > indices[i]) {
547  int itmp = indices[i];
548  indices[i] = indices[j];
549  indices[j] = itmp;
550  }
551 
552  // remove rows and columns starting from largest indices
553  for (int i=0;i<n;i++) {
554  if (indices[i] > m_size-1) {
555  throw std::out_of_range( "AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
556  }
557  RemoveDoF(indices[i]);
558  }
559 
560  return m_size;
561 }
562 
563 //______________________________________________________________________________
564 void AlSpaMat::RemoveAlignPar(int index, int control)
565 {
566  // Dynamic shiftRow
567  int shiftRow = 1;
568  int counterRow = 0;
569 
570  // Dynamic shiftCol
571  int shiftCol;
572  int counterCol = 0;
573 
574  index = index-control;
575 
577 
578  for(int row=index; row<(size()-shiftRow); row++) {
579  shiftCol = 0;
580  counterCol = 0;
581 
582  for(int col=0; col<row+1; col++) {
583 
584  if (col==index) {
585  shiftCol = 1;
586  counterCol=0;
587  }
588 
589  pos = m_ptr_map.find(std::make_pair(row+shiftRow,col+shiftCol));
590  if( pos!=m_ptr_map.end() ) { // exists
591  elemr(row,col) = pos->second;
592  m_ptr_map.erase(pos);
593  }
594  else {
595  pos = m_ptr_map.find(std::make_pair(row,col));
596  if( pos!=m_ptr_map.end() )
597  m_ptr_map.erase(pos);
598  }
599  counterCol++;
600  if (counterCol==5-control) {
601  counterCol=0;
602  shiftCol++;
603  }
604  }
605 
606  counterRow++;
607  if (counterRow==5-control) {
608  counterRow=0;
609  shiftRow++;
610  }
611  }
612 
613  m_nele = m_ptr_map.size();
614 }
615 
616 //______________________________________________________________________________
618 {
619  // index here is module number within matrix (size=6*max_index)
620  const int shift=6;
621  const int ind=shift*index;
622 
623  RemoveDoF(ind, shift);
624 }
625 
626 //______________________________________________________________________________
627 void AlSpaMat::reSize(long int n)
628 {
629  // balanced binary tree!
630  if( n!=m_size ) {
631  // make a temporary copy:
632  AlSpaMat m(*this);
633  m_ptr_map.clear();
634  m_size = n;
635  long int i, j;
636  for (const datamap::value_type& p : m.m_ptr_map) {
637  m.elem(p.first, i, j);
638  if( i<n && j<n )
639  m_ptr_map.insert(p);
640  }
641  }
642 
643  m_nele = m_ptr_map.size();
644 }
645 
646 //______________________________________________________________________________
647 void AlSpaMat::SetPathBin(const std::string &path)
648 {
649  m_pathbin.assign(path);
650 }
651 
652 //______________________________________________________________________________
653 void AlSpaMat::SetPathTxt(const std::string &path)
654 {
655  m_pathtxt.assign(path);
656 }
657 
658 //______________________________________________________________________________
659 //jlove StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
660 // bool square, double scale, float version){
661 StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
662  bool /*square*/, double, float version)
663 {
664  std::ofstream outmat;
665 
666  int32_t msizz = 1000000+size();
667  int32_t nelem = m_ptr_map.size();
668 
669  if(binary) {
670  outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
671  if(outmat.fail())
672  return StatusCode::FAILURE;
673  // add 10^6 to the size to distinguish
674  // from dense format
675  outmat.write(charAddress(msizz), sizeof (msizz));
676  outmat.write(charAddress(version), sizeof (version));
677  outmat.write(charAddress(nelem), sizeof (nelem));
678  }
679  else {
680  outmat.open((m_pathtxt+filename).c_str());
681  if(outmat.fail())
682  return StatusCode::FAILURE;
683  outmat.setf(std::ios::fixed);
684  outmat.setf(std::ios::showpoint);
685  outmat.precision(3);
686  outmat << "msizz: " << std::setw(6) << msizz << std::endl;
687  outmat << "nelem: " << std::setw(6) << nelem << std::endl;
688  outmat << "AlSpaMat version: " << std::setw(6) << version << std::endl;
689  }
690 
691  double melem=0;
692  long int ii, jj;
693  int32_t i, j;
694 
695  for (const datamap::value_type& p : m_ptr_map) {
696  melem = p.second;
697  elem(p.first, ii, jj); i=ii; j=jj; // just a type conversion
698  if(binary) {
699  outmat.write(charAddress((i)), sizeof (i));
700  outmat.write(charAddress((j)), sizeof (j));
701  outmat.write(charAddress((melem)), sizeof (melem));
702  }
703  else
704  outmat << std::setw(6) << i << std::setw(6) << j << std::setw(18) << melem << std::endl;
705  }
706  outmat.close();
707 
708  return StatusCode::SUCCESS;
709 }
710 
711 //______________________________________________________________________________
712 StatusCode AlSpaMat::CheckMatVersion(const std::string& filename, bool &StdUnits)
713 {
714  std::ifstream inmat((filename).c_str(), std::ios::binary);
715  if(inmat.fail())
716  return StatusCode::FAILURE;
717 
718  m_ptr_map.clear();
719 
720  int32_t msiz=0;
721  inmat.read(charAddress(msiz), sizeof (msiz));
722 
723  float version=0.0;
724  inmat.read(charAddress(version), sizeof (version));
725 
726  StdUnits = version>=2.0;
727 
728  inmat.close();
729 
730  return StatusCode::SUCCESS;
731 }
732 
733 //______________________________________________________________________________
734 StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang, float &version)
735 {
736  bool stdUnits = true;
737  if (StatusCode::SUCCESS != CheckMatVersion(m_pathbin+filename, stdUnits))
738  return StatusCode::FAILURE;
739 
740  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
741  if(inmat.fail())
742  return StatusCode::FAILURE;
743 
744  m_ptr_map.clear();
745 
746  int32_t msiz=0;
747  int32_t nelem;
748  inmat.read(charAddress(msiz), sizeof (msiz));
749  if( msiz>999999 )
750  dofs = msiz-1000000;
751  else
752  dofs = abs(msiz);
753  m_size = dofs;
754 
755  if (stdUnits)
756  inmat.read(charAddress(version), sizeof (version));
757 
758  double melem=0;
759  int32_t i, j;
760 
761  if(msiz>999999) { // sparse format
762  triang=false;
763  inmat.read(charAddress(nelem), sizeof (nelem));
764  m_nele=nelem;
765  for(int k=0; k<nelem; k++) {
766  inmat.read(charAddress(i), sizeof (i));
767  inmat.read(charAddress(j), sizeof (j));
768  inmat.read(charAddress(melem), sizeof (melem));
769  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
770  }
771  }
772  else if(msiz>0) { // square format
773  triang=false;
774  for(int32_t i=0; i<msiz; i++) {
775  for(int32_t j=0; j<msiz; j++) {
776  inmat.read(charAddress(melem), sizeof (melem));
777  if( i>=j && melem!=0. )
778  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
779  }
780  }
781  }
782  else { // triangular format
783  triang=true;
784  msiz = (-1)*msiz;
785  for( int32_t i=0; i<msiz; i++) {
786  for( int32_t j=0; j<=i; j++) {
787  inmat.read(charAddress(melem), sizeof (melem));
788  if( melem!=0. )
789  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
790  }
791  }
792  }
793 
794  m_nele = m_ptr_map.size();
795  inmat.close();
796  return StatusCode::SUCCESS;
797 }
798 
799 //______________________________________________________________________________
800 StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
801  bool &triang, float &version)
802 {
803  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
804  if(inmat.fail())
805  return StatusCode::FAILURE;
806 
807  m_ptr_map.clear();
808 
809  int32_t msiz=0;
810  int32_t nelem;
811  inmat.read(charAddress(msiz), sizeof (msiz));
812  if( msiz>999999 )
813  dofs = msiz-1000000;
814  else
815  dofs = std::abs(msiz);
816  m_size = dofs;
817 
818  inmat.read(charAddress(version), sizeof (version));
819 
820  double melem=0;
821  int32_t i, j;
822 
823  if(msiz>999999) { // sparse format
824  triang=false;
825  inmat.read(charAddress(nelem), sizeof (nelem));
826  m_nele=nelem;
827  for(int k=0; k<nelem; k++) {
828  inmat.read(charAddress(i), sizeof (i));
829  inmat.read(charAddress(j), sizeof (j));
830  inmat.read(charAddress(melem), sizeof (melem));
831  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
832  }
833  }
834  else if(msiz>0) { // square format
835  triang=false;
836  for(int32_t i=0; i<msiz; i++) {
837  for(int32_t j=0; j<msiz; j++) {
838  inmat.read(charAddress(melem), sizeof (melem));
839  if( i>=j && melem!=0. )
840  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
841  }
842  }
843  }
844  else { // triangular format
845  triang=true;
846  msiz = (-1)*msiz;
847  for( int32_t i=0; i<msiz; i++) {
848  for( int32_t j=0; j<=i; j++) {
849  inmat.read(charAddress(melem), sizeof (melem));
850  if( melem!=0. )
851  m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
852  }
853  }
854  }
855 
856  m_nele = m_ptr_map.size();
857  inmat.close();
858  return StatusCode::SUCCESS;
859 }
860 
862 {
863 
864  long int nonZeroElements = m_ptr_map.size();
865  int *irow = new int[nonZeroElements];
866  int *icol = new int[nonZeroElements];
867  double *val = new double[nonZeroElements];
868 
869  long int i, j;
870  long int counter(0);
871  for (const datamap::value_type& p : m_ptr_map) {
872  i = p.first.first;
873  j = p.first.second;
874  *(val+counter)=p.second;
875  *(irow+counter)= i;
876  *(icol+counter)= j;
877  counter++;
878  }
879 
880  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,m_size-1,0,m_size-1);
881  myTMatrix->SetMatrixArray(nonZeroElements,irow,icol,val);
882 
883  delete [] irow;
884  delete [] icol;
885  delete [] val;
886  return myTMatrix;
887 }
888 
889 
890 
891 } // end namespace Trk
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
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
Trk::AlSpaMat::SetPathBin
virtual void SetPathBin(const std::string &) override final
Definition: AlSpaMat.cxx:647
Trk::AlSpaMat::ReadProjected
virtual StatusCode ReadProjected(const std::string &, int &, bool &, float &) override final
Definition: AlSpaMat.cxx:800
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
Trk::AlSpaMat::RemoveCollsRows
virtual int RemoveCollsRows(std::vector< int >) override final
Definition: AlSpaMat.cxx:531
Trk::AlSpaMat::Write
virtual StatusCode Write(const std::string &, bool, bool, double, float) override final
Definition: AlSpaMat.cxx:661
Trk::AlVec
Definition: AlVec.h:23
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
AlSpaMat.h
index
Definition: index.py:1
Trk::AlSymMatBase
Definition: AlSymMatBase.h:38
hist_file_dump.d
d
Definition: hist_file_dump.py:142
Trk::AlSpaMat::elemr
virtual double & elemr(long int, long int) override final
Definition: AlSpaMat.cxx:143
Trk::indices
std::pair< long int, long int > indices
Definition: AlSymMatBase.h:24
Trk::AlSpaMat::copy
void copy(const AlSpaMat &m)
Definition: AlSpaMat.cxx:92
Trk::AlSymMatBase::m_ptr_map
datamap m_ptr_map
Definition: AlSymMatBase.h:101
Trk::AlSpaMat::operator*
AlMat operator*(const AlSymMatBase &) const
Definition: AlSpaMat.cxx:331
Trk::AlSpaMat::SolveWithEigen
int SolveWithEigen(AlVec &RHS)
Definition: AlSpaMat.cxx:408
AlVec.h
Trk::AlSpaMat::elemc
virtual double elemc(long int, long int) const override final
Definition: AlSpaMat.cxx:179
keylayer_zslicemap.row
row
Definition: keylayer_zslicemap.py:155
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
Trk::AlSpaMat::determinant
virtual double determinant() override final
Definition: AlSpaMat.cxx:403
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
Trk::AlSpaMat::m_pathbin
std::string m_pathbin
Definition: AlSpaMat.h:85
AlSymMat.h
Trk::AlSpaMat::CheckMatVersion
virtual StatusCode CheckMatVersion(const std::string &, bool &) override final
Definition: AlSpaMat.cxx:712
Trk::AlSpaMat::operator=
AlSpaMat & operator=(const double &)
Definition: AlSpaMat.cxx:264
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::AlSymMat
Definition: AlSymMat.h:26
TMatrixDSparse
class TMatrixTSparse< double > TMatrixDSparse
Definition: AlSpaMat.h:14
Trk::AlSpaMat::RemoveAlignPar
virtual void RemoveAlignPar(int, int) override final
Definition: AlSpaMat.cxx:564
Trk::AlSymMatBase::size
long int size() const
Definition: AlSymMatBase.h:145
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::AlVec::size
int size() const
Definition: AlVec.h:109
Trk::AlSpaMat::makeTMatrix
virtual TMatrixDSparse * makeTMatrix() override final
Definition: AlSpaMat.cxx:861
Trk::AlSpaMat::diagonalize
virtual int diagonalize(char jobz, AlVec &w, AlMat &z) override final
Definition: AlSpaMat.cxx:481
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::AlSpaMat::invert
virtual int invert() override final
Definition: AlSpaMat.cxx:487
AlMat.h
Trk::AlSpaMat::AlSpaMat
AlSpaMat()
Definition: AlSpaMat.cxx:43
Trk::AlMat
Definition: AlMat.h:27
Trk::AlSpaMat::~AlSpaMat
~AlSpaMat()
Definition: AlSpaMat.cxx:84
Trk::AlSpaMat::reSize
virtual void reSize(long int) override final
Definition: AlSpaMat.cxx:627
Trk::AlSymMatBase::m_nele
long int m_nele
Definition: AlSymMatBase.h:104
Trk::AlSpaMat::RemoveDoF
void RemoveDoF(int, int nelem=1)
Definition: AlSpaMat.cxx:493
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Trk::AlSpaMat::operator+
AlSpaMat operator+(const AlSpaMat &) const
Definition: AlSpaMat.cxx:273
EventPrimitives.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
Trk::mapiterator
datamap::iterator mapiterator
Definition: AlSymMatBase.h:29
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
python.PyAthena.v
v
Definition: PyAthena.py:154
get_generator_info.version
version
Definition: get_generator_info.py:33
Trk::AlSpaMat::operator+=
AlSpaMat & operator+=(const AlSpaMat &)
Definition: AlSpaMat.cxx:288
DeMoScan.index
string index
Definition: DeMoScan.py:362
a
TList * a
Definition: liststreamerinfos.cxx:10
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
Trk::AlSpaMat::operator*=
AlSpaMat & operator*=(const double &)
Definition: AlSpaMat.cxx:384
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:23
Trk::AlSpaMat::SetPathTxt
virtual void SetPathTxt(const std::string &) override final
Definition: AlSpaMat.cxx:653
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
Trk::AlSpaMat::operator-=
AlSpaMat & operator-=(const AlSpaMat &)
Definition: AlSpaMat.cxx:317
Trk::AlSpaMat
Definition: AlSpaMat.h:27
Trk::v
@ v
Definition: ParamDefs.h:78
Trk::AlSpaMat::RemoveModule
virtual void RemoveModule(int) override final
Definition: AlSpaMat.cxx:617
fitman.k
k
Definition: fitman.py:528
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Trk::AlSpaMat::Read
virtual StatusCode Read(const std::string &, int &, bool &, float &) override final
Definition: AlSpaMat.cxx:734
Trk::AlSpaMat::operator-
AlSpaMat operator-(const AlSpaMat &) const
Definition: AlSpaMat.cxx:302
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37