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 
33 namespace Trk {
34 
35 //______________________________________________________________________________
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 }
44 
45 //______________________________________________________________________________
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 }
54 
55 //______________________________________________________________________________
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 }
65 
66 //______________________________________________________________________________
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 }
75 
76 //______________________________________________________________________________
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 }
83 
84 //______________________________________________________________________________
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 }
94 
95 //______________________________________________________________________________
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  }
113 
114 //______________________________________________________________________________
115 void AlSpaMat::copy(const AlMat& m)
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  }
134 
135 //______________________________________________________________________________
136 double& AlSpaMat::elemr(long int i,long int j)
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 }
170 
171 //______________________________________________________________________________
172 double AlSpaMat::elemc(long int i,long int j) const
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 }
200 
201 //______________________________________________________________________________
202 indices AlSpaMat::elem(long int i,long int j)
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 }
223 
224 //______________________________________________________________________________
225 // cppcheck-suppress operatorEqVarError; m_ptr_row/col deliberately not copied.
227 {
228 
229  if ( this != &m ) {
230  m_size=m.size();
231  copy(m);
232  }
233  return *this;
234 }
235 
236 //______________________________________________________________________________
238 {
239  m_size=m.size();
240  copy(m);
241  return *this;
242 }
243 
244 //______________________________________________________________________________
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 }
255 
256 //______________________________________________________________________________
258 {
259  for (datamap::value_type& p : m_ptr_map)
260  p.second = d;
261 
262  return *this;
263 }
264 
265 //______________________________________________________________________________
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 }
279 
280 //______________________________________________________________________________
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 }
293 
294 //______________________________________________________________________________
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 }
308 
309 //______________________________________________________________________________
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 }
322 
323 //______________________________________________________________________________
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 }
340 
341 //______________________________________________________________________________
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 }
358 
359 //______________________________________________________________________________
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 }
375 
376 //______________________________________________________________________________
378 {
379  for (datamap::value_type& p : m_ptr_map)
380  p.second *= d;
381 
382  return *this;
383 }
384 
385 //______________________________________________________________________________
386 AlSpaMat AlSpaMat::operator*(const double& d) const
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 }
394 
395 //______________________________________________________________________________
397 {
398  throw std::invalid_argument( "AlSpaMat::determinant: not implemented!" );
399 }
400 
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 }
469 
470 
471 
472 //______________________________________________________________________________
473 //jlove int AlSpaMat::diagonalize(char jobz, AlVec& w, AlMat& z) {
475 {
476  throw std::invalid_argument( "AlSpaMat::diagonalize: not implemented!" );
477 }
478 
479 //______________________________________________________________________________
481 {
482  throw std::invalid_argument( "AlSpaMat::invert: not implemented!" );
483 }
484 
485 //______________________________________________________________________________
486 void AlSpaMat::RemoveDoF(int index, int nelem)
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 }
522 
523 //______________________________________________________________________________
524 int AlSpaMat::RemoveCollsRows(std::vector<int> indices)
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 }
555 
556 //______________________________________________________________________________
557 void AlSpaMat::RemoveAlignPar(int index, int control)
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 }
608 
609 //______________________________________________________________________________
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 }
618 
619 //______________________________________________________________________________
620 void AlSpaMat::reSize(long int n)
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 }
638 
639 //______________________________________________________________________________
640 void AlSpaMat::SetPathBin(const std::string &path)
641 {
642  m_pathbin.assign(path);
643 }
644 
645 //______________________________________________________________________________
646 void AlSpaMat::SetPathTxt(const std::string &path)
647 {
648  m_pathtxt.assign(path);
649 }
650 
651 //______________________________________________________________________________
652 //jlove StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
653 // bool square, double scale, float version){
654 StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
655  bool /*square*/, double, float version)
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 }
703 
704 //______________________________________________________________________________
705 StatusCode AlSpaMat::CheckMatVersion(const std::string& filename, bool &StdUnits)
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 }
725 
726 //______________________________________________________________________________
727 StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang, float &version)
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 }
791 
792 //______________________________________________________________________________
793 StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
794  bool &triang, float &version)
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 }
853 
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 }
881 
882 
883 
884 } // end namespace Trk
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
query_example.row
row
Definition: query_example.py:24
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:640
Trk::AlSpaMat::ReadProjected
virtual StatusCode ReadProjected(const std::string &, int &, bool &, float &) override final
Definition: AlSpaMat.cxx:793
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
Trk::AlSpaMat::RemoveCollsRows
virtual int RemoveCollsRows(std::vector< int >) override final
Definition: AlSpaMat.cxx:524
Trk::AlSpaMat::Write
virtual StatusCode Write(const std::string &, bool, bool, double, float) override final
Definition: AlSpaMat.cxx:654
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: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::operator*
AlMat operator*(const AlSymMatBase &) const
Definition: AlSpaMat.cxx:324
Trk::AlSpaMat::SolveWithEigen
int SolveWithEigen(AlVec &RHS)
Definition: AlSpaMat.cxx:401
AlVec.h
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
Trk::AlSpaMat::determinant
virtual double determinant() override final
Definition: AlSpaMat.cxx:396
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:705
Trk::AlSpaMat::operator=
AlSpaMat & operator=(const double &)
Definition: AlSpaMat.cxx:257
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:557
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
Trk::AlVec::size
int size() const
Definition: AlVec.h:109
Trk::AlSpaMat::makeTMatrix
virtual TMatrixDSparse * makeTMatrix() override final
Definition: AlSpaMat.cxx:854
Trk::AlSpaMat::diagonalize
virtual int diagonalize(char jobz, AlVec &w, AlMat &z) override final
Definition: AlSpaMat.cxx:474
beamspotman.n
n
Definition: beamspotman.py:731
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:480
AlMat.h
Trk::AlSpaMat::AlSpaMat
AlSpaMat()
Definition: AlSpaMat.cxx:36
Trk::AlMat
Definition: AlMat.h:27
Trk::AlSpaMat::~AlSpaMat
~AlSpaMat()
Definition: AlSpaMat.cxx:77
Trk::AlSpaMat::reSize
virtual void reSize(long int) override final
Definition: AlSpaMat.cxx:620
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
Trk::AlSpaMat::operator+
AlSpaMat operator+(const AlSpaMat &) const
Definition: AlSpaMat.cxx:266
EventPrimitives.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
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
Trk::AlSpaMat::operator+=
AlSpaMat & operator+=(const AlSpaMat &)
Definition: AlSpaMat.cxx:281
DeMoScan.index
string index
Definition: DeMoScan.py:364
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:377
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
Trk::AlSpaMat::SetPathTxt
virtual void SetPathTxt(const std::string &) override final
Definition: AlSpaMat.cxx:646
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:310
Trk::AlSpaMat
Definition: AlSpaMat.h:27
checkFileSG.ind
list ind
Definition: checkFileSG.py:118
Trk::v
@ v
Definition: ParamDefs.h:78
Trk::AlSpaMat::RemoveModule
virtual void RemoveModule(int) override final
Definition: AlSpaMat.cxx:610
fitman.k
k
Definition: fitman.py:528
Trk::AlSpaMat::Read
virtual StatusCode Read(const std::string &, int &, bool &, float &) override final
Definition: AlSpaMat.cxx:727
Trk::AlSpaMat::operator-
AlSpaMat operator-(const AlSpaMat &) const
Definition: AlSpaMat.cxx:295
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37