ATLAS Offline Software
AlMat.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 #include "GaudiKernel/StatusCode.h"
6 
11 
12 #include <algorithm>
13 #include <cstdint>
14 #include <exception>
15 #include <fstream>
16 #include <iomanip>
17 #include <sstream>
18 
19 namespace{
20 char *
21 charAddress(auto &val){
22  return reinterpret_cast<char *> (&val);
23 }
24 
25 }
26 
27 namespace Trk {
28 
30  : m_ncol(0)
31  , m_nrow(0)
32  , m_nele(0)
33  , m_pathbin("./")
34  , m_pathtxt("./")
35  , m_ptr_data(nullptr)
36  , m_transpose(false)
37 {
38 }
39 
40 AlMat::AlMat(int N, int M)
41  : m_ncol(M)
42  , m_nrow(N)
43  , m_nele(N * M)
44  , m_pathbin("./")
45  , m_pathtxt("./")
46  , m_ptr_data(new double[m_nele])
47  , m_transpose(false)
48 {
49 
50  double* p = m_ptr_data + m_nele;
51  while (p > m_ptr_data) *(--p) = 0.;
52 }
53 
55  : m_ncol(m.ncol())
56  , m_nrow(m.nrow())
57  , m_nele(m_nrow * m_ncol)
58  , m_pathbin(m.m_pathbin)
59  , m_pathtxt(m.m_pathtxt)
60  , m_ptr_data(new double[m_nele])
61  , m_transpose(false)
62 {
63 
64  copy(m);
65 }
66 
68 {if( m_ptr_data != nullptr ) delete [] m_ptr_data;}
69 
70 void AlMat::copy(const AlMat& m) {
71  int nr = nrow();
72  int nc = ncol();
73  if( nr != m.nrow() || nc != m.ncol() ) {
74  throw std::range_error( "AlMat::copy: sizes do not match!" );
75  }
76 
77  for( int i=0; i<nr; i++ ) {
78  for( int j=0; j<nc; j++ ) {
79  elemr(i, j) = m.elemc(i, j);
80  }
81  }
82 }
83 
84 void AlMat::copy(const AlSymMatBase& m) {
85  long int n = m.size();
86  if( nrow() != n || ncol() != n ) {
87  throw std::range_error( "AlMat::copy: sizes do not match!" );
88  }
89 
90  for( int i=0; i<n; i++ ) {
91  for( int j=0; j<n; j++ ) {
92  *(m_ptr_data+i*m_ncol+j) = m.elemc(i, j);
93  }
94  }
95 }
96 
97 double& AlMat::elemr(long int i,long int j) {
98 #ifdef _DEBUG
99  if( i<0 ) { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
100  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
101  if( j<0 ) { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
102  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" ); };
103 #endif
104 
105  if( m_transpose ) return *(m_ptr_data+j*m_ncol+i);
106  return *(m_ptr_data+i*m_ncol+j);
107 }
108 
109 double AlMat::elemc(long int i,long int j) const {
110 #ifdef _DEBUG
111  if( i<0 ) { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
112  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
113  if( j<0 ) { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
114  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" ); };
115 #endif
116 
117  if( m_transpose ) return *(m_ptr_data+j*m_ncol+i);
118  return *(m_ptr_data+i*m_ncol+j);
119 }
120 
121 long int AlMat::elem(long int i,long int j) const {
122 #ifdef _DEBUG
123  if( i<0 ) { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
124  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
125  if( j<0 ) { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
126  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" ); };
127 #endif
128 
129  if( m_transpose ) return (j*m_ncol+i);
130  return (i*m_ncol+j);
131 }
132 
133 
134 AlMat& AlMat::operator=(const double& d) {
135 
136  double* p = m_ptr_data + m_nele;
137  while (p > m_ptr_data) *(--p) = d;
138 
139  return *this;
140 }
141 
143  if (this==&m) return *this;
144 
145  if(( m_nrow!=0 || m_ncol!=0) && (m_nrow != m.nrow() || m_ncol != m.ncol() ))
146  {
147  throw std::range_error( "AlMat=AlMat Assignment: size do not match!" );
148  }
149 
150  if ( m_ptr_data != m.m_ptr_data ) {
151  reSize(m.nrow(), m.ncol());
152  copy(m);
153  };
154  return (*this);
155 }
156 
158  if( ( m_nrow!=0 || m_ncol!=0) && (m_nrow != m.size() || m_ncol != m.size() ))
159  {
160  throw std::range_error( "AlMat=AlSymMatBase Assignment: size do not match!" );
161  }
162 
163  if ( m_ptr_data != m.ptrData() ) {
164  reSize(m.size(), m.size());
165  copy(m);
166  }
167  return (*this);
168 }
169 
170 AlMat AlMat::operator+(const AlMat& m) const {
171  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
172  {
173  throw std::range_error( "AlMat: operator+: size do not match!" );
174  }
175 
176  AlMat b( m_nrow, m_ncol );
177 
178  double* p = m_ptr_data + m_nele;
179  double* q = m.m_ptr_data + m_nele;
180  double* r = b.m_ptr_data + m_nele;
181  while (p > m_ptr_data) *(--r) = (*(--p))+(*(--q));
182 
183  return b;
184 }
185 
187 
188  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
189  {
190  throw std::range_error( "AlMat: operator+=: size do not match!" );
191  }
192 
193  double* p = m_ptr_data + m_nele;
194  double* q = m.m_ptr_data + m_nele;
195  while (p > m_ptr_data) *(--p) += *(--q);
196 
197  return *this;
198 }
199 
200 AlMat AlMat::operator-(const AlMat& m) const {
201  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
202  {
203  throw std::range_error( "AlMat: operator-: size do not match!" );
204  }
205 
206  AlMat b( m_nrow, m_ncol );
207 
208  double* p = m_ptr_data + m_nele;
209  double* q = m.m_ptr_data + m_nele;
210  double* r = b.m_ptr_data + m_nele;
211  while (p > m_ptr_data) *(--r) = (*(--p))-(*(--q));
212 
213  return b;
214 }
215 
217  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
218  {
219  throw std::range_error( "AlMat: operator-=: size do not match!");
220  }
221 
222  double* p = m_ptr_data + m_nele;
223  double* q = m.m_ptr_data + m_nele;
224  while (p > m_ptr_data) *(--p) -= *(--q);
225 
226  return *this;
227 }
228 
229 AlMat AlMat::operator*(const AlMat& m) const {
230  if( ncol() != m.nrow() )
231  {
232  throw std::range_error( "AlMat: operator*: size do not match!" );
233  }
234 
235  int k(nrow());
236  int l(m.ncol());
237  int f(ncol());
238 
239  AlMat b( k, l);
240  for( int i=0; i<k; i++ ) {
241  for( int j=0; j<l; j++ ) {
242  for( int n=0; n<f; n++) b.elemr(i, j) += elemc(i, n)*m.elemc(n, j);
243  }
244  }
245  return b;
246 }
247 
249  if( ncol() != m.size())
250  {
251  throw std::range_error( "AlMat: operator*: size do not match!" );
252  }
253 
254  int k(nrow());
255  int l(m.size());
256 
257  AlMat b( k, l);
258  for( int i=0; i<k; i++ ) {
259  for( int j=0; j<l; j++ ) {
260  for( int n=0; n<l; n++) b.elemr(i, j) += elemc(i, n)*m.elemc(n, j);
261  }
262  }
263  return b;
264 }
265 
266 
267 AlVec AlMat::operator*(const AlVec& v) const {
268  if( ncol() != v.size() )
269  {
270  throw std::range_error( "AlMat: operator*: size do not match! " );
271  }
272 
273  int k(nrow());
274  int l(ncol());
275  double* p;
276  double* q;
277 
278  AlVec b(k);
279  for( int i=0; i<k; i++ ) {
280  p = b.ptrData()+i;
281  q = m_ptr_data+i*l;
282  for( int j=0; j<l; j++ ) *p += (*(q+j))*(*(v.ptrData()+j));
283  }
284 
285  return b;
286 }
287 
288 AlMat& AlMat::operator*=(const double& d) {
289  double* p = m_ptr_data+m_nele;
290  while (p > m_ptr_data)
291  *(--p) *= d;
292 
293  return *this;
294 }
295 
296 // transposition
297 AlMat AlMat::T() const {
298  AlMat b(m_ncol,m_nrow);
299  for( int i=0; i<b.m_nrow; i++ ) {
300  for( int j=0; j<b.m_ncol; j++ ) {
301  b[i][j]= *(m_ptr_data+j*m_ncol+i);
302  }
303  }
304  return b;
305 }
306 
307 // transposition
309 
310  m_transpose = true;
311 
312  return *this;
313 }
314 
315 // normal position
317 
318  m_transpose = false;
319 
320  return *this;
321 }
322 
323 
324 //invert sym matrix declared as non-symetric for convenience
325 void AlMat::invertS(int& ierr, double Norm = 1.) {
326  if(m_nrow!=m_ncol) {
327  throw std::range_error( "AlMat invertS: non-square matrix!" );
328  }
329 
330  AlSymMat b(m_nrow);
331  for( int i=0; i<b.size(); i++ ) {
332  for( int j=0; j<=i; j++ ) {
333  b.elemr(i, j) = elemc(i, j);
334  }
335  }
336 
337  b*=Norm;
338  ierr = b.invert();
339  b*=Norm;
340 
341  if (ierr==0) (*this).copy(b);
342 
343 }
344 
345 // reSize
346 void AlMat::reSize(int Nnew, int Mnew) {
347  if ( Nnew != m_nrow || Mnew != m_ncol ) {
348 
349  double* p = m_ptr_data;
350  m_nele = Nnew*Mnew;
351  m_ptr_data = new double[m_nele];
352  int nrow_old = m_nrow;
353  int ncol_old = m_ncol;
354  m_transpose = false;
355 
356  m_nrow = Nnew;
357  m_ncol = Mnew;
358  int k = m_nrow <= nrow_old ? m_nrow : nrow_old;
359  int l = m_ncol <= ncol_old ? m_ncol : ncol_old;
360 
361  for( int i=0; i<k; i++ ) {
362  for( int j=0; j<l; j++ ) {
363  *(m_ptr_data+i*m_ncol+j) = *(p+i*ncol_old+j);
364  }
365  }
366 
367  delete [] p;
368  }
369 }
370 
371 
373  std::ifstream inmat(filename.c_str(), std::ios::binary);
374  if(inmat.fail())
375  return StatusCode::FAILURE;
376 
377  int32_t mNrow = m_nrow;
378  inmat.read(charAddress(mNrow), sizeof (mNrow));
379  int32_t mNcol = m_ncol;
380  inmat.read(charAddress(mNcol), sizeof (mNcol));
381 
382  m_nrow=std::abs(mNrow);
383  m_ncol=std::abs(mNcol);
385  m_transpose = false;
386 
387 
388 
389  double melem=0;
390  for(int i=0; i<m_nrow; i++) {
391  for(int j=0; j<m_ncol; j++) {
392  inmat.read(charAddress(melem), sizeof (melem));
393  *(m_ptr_data+i*m_ncol+j) = melem;
394  }
395  }
396 
397  inmat.close();
398  return StatusCode::SUCCESS;
399 }
400 
401 //______________________________________________________________________________
402 void AlMat::SetPathBin(const std::string &path)
403 {
404  m_pathbin.assign(path);
405 }
406 
407 //______________________________________________________________________________
408 void AlMat::SetPathTxt(const std::string &path)
409 {
410  m_pathtxt.assign(path);
411 }
412 
413 //______________________________________________________________________________
414 StatusCode AlMat::Write(const std::string &filename, bool binary, unsigned int precision)
415 {
416  std::ofstream outmat;
417 
418  if(binary) {
419  outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
420  if(outmat.fail())
421  return StatusCode::FAILURE;
422 
423  int32_t mNrow = m_nrow;
424  outmat.write(charAddress(mNrow), sizeof (mNrow));
425  int32_t mNcol = m_ncol;
426  outmat.write(charAddress(mNcol), sizeof (mNcol));
427  }
428  else {
429  outmat.open((m_pathtxt+filename).c_str());
430  if(outmat.fail())
431  return StatusCode::FAILURE;
432  outmat.setf(std::ios::fixed);
433  if(precision>0)
434  outmat.setf(std::ios::showpoint);
435  outmat.precision(precision);
436  }
437 
438  double melem=0;
439 
440  for( int i=0; i<m_nrow; i++) {
441  for( int j=0; j<m_ncol; j++) {
442  melem = *(m_ptr_data+i*m_ncol+j);
443  if(binary)
444  outmat.write(charAddress(melem), sizeof (melem));
445  else
446  outmat << " " << std::setw(12) << melem;
447  }
448  if(!binary)
449  outmat << std::endl;
450  }
451 
452  outmat.close();
453  return StatusCode::SUCCESS;
454 }
455 
456 //______________________________________________________________________________
457 const std::string AlMat::Print(const int NColsPerSet)
458 {
459  std::ostringstream textmatrix;
460  int FirstCol = 0;
461  int LastCol = NColsPerSet;
462 
463  int NSets = (ncol()/NColsPerSet);
464  if (ncol()%NColsPerSet>0) NSets++;
465  textmatrix << std::endl;
466  textmatrix << " ++++++++++++++ AlMat ++ "
467  << nrow()
468  << " x "
469  << ncol()
470  << " +++++++++++++++++++++++++++++"<< std::endl;
471  for (int set=0; set < NSets; set++) {
472  FirstCol = set*NColsPerSet;
473  LastCol = (set+1)*NColsPerSet;
474  if (LastCol>ncol()) LastCol = ncol();
475  textmatrix << " |row\\col| ";
476  for (int col=FirstCol; col < LastCol; col++) {
477  textmatrix << std::setw(6) << col << " | ";
478  }
479  textmatrix << std::endl;
480  textmatrix << " |-------|";
481  for (int col=FirstCol; col < LastCol; col++) {
482  textmatrix << "--------|";
483  }
484  textmatrix << std::endl;
485  for (int row=0; row<nrow(); row++) {
486  textmatrix << " |" << std::setw(6) << row << " | ";
487  for (int col=FirstCol; col<LastCol; col++) {
488  double melem = *(m_ptr_data+row*ncol()+col);
489  textmatrix << std::setprecision(5) << std:: setw(6) << melem << " | ";
490  }
491  textmatrix << std::endl;
492  }
493  if (set != (NSets-1)) textmatrix << std::endl;
494  }
495  textmatrix << " +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<< std::endl;
496  return textmatrix.str();
497 }
498 
499 } // end namespace Trk
500 
Trk::AlMat::operator=
AlMat & operator=(const double &)
Definition: AlMat.cxx:134
beamspotman.r
def r
Definition: beamspotman.py:674
Trk::AlMat::m_pathbin
std::string m_pathbin
Definition: AlMat.h:112
Trk::AlMat::T
AlMat T() const
Definition: AlMat.cxx:297
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
Trk::AlMat::AlMat
AlMat()
Definition: AlMat.cxx:29
Trk::AlVec
Definition: AlVec.h:23
AlSpaMat.h
Trk::AlSymMatBase
Definition: AlSymMatBase.h:38
hist_file_dump.d
d
Definition: hist_file_dump.py:142
Trk::AlMat::m_transpose
bool m_transpose
Definition: AlMat.h:115
Trk::AlMat::operator*
AlMat operator*(const AlMat &m) const
Definition: AlMat.cxx:229
CSV_InDetExporter.new
new
Definition: CSV_InDetExporter.py:145
Trk::AlMat::reSize
void reSize(int Nnew, int Mnew)
Definition: AlMat.cxx:346
Trk::AlMat::ReadScalaPack
StatusCode ReadScalaPack(const std::string &)
Definition: AlMat.cxx:372
AlVec.h
keylayer_zslicemap.row
row
Definition: keylayer_zslicemap.py:155
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
Trk::AlMat::SetPathTxt
void SetPathTxt(const std::string &)
Definition: AlMat.cxx:408
Trk::AlMat::operator+
AlMat operator+(const AlMat &m) const
Definition: AlMat.cxx:170
Trk::AlMat::m_nrow
int m_nrow
Definition: AlMat.h:110
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
Trk::AlMat::Print
const std::string Print(const int NColsPerSet=10)
Definition: AlMat.cxx:457
AlSymMat.h
Trk::AlMat::ncol
long int ncol() const
Definition: AlMat.h:162
Trk::AlSymMat
Definition: AlSymMat.h:26
Trk::AlMat::operator+=
AlMat & operator+=(const AlMat &m)
Definition: AlMat.cxx:186
Trk::AlMat::invertS
void invertS(int &ierr, double Norm)
Definition: AlMat.cxx:325
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
Trk::AlMat::Transpose
AlMat & Transpose()
Definition: AlMat.cxx:308
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::AlMat::copy
void copy(const AlMat &m)
Definition: AlMat.cxx:70
beamspotman.n
n
Definition: beamspotman.py:729
Trk::AlMat::operator-=
AlMat & operator-=(const AlMat &m)
Definition: AlMat.cxx:216
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::AlMat::SetPathBin
void SetPathBin(const std::string &)
Definition: AlMat.cxx:402
Trk::AlMat::m_ncol
int m_ncol
Definition: AlMat.h:109
JetVoronoiDiagramHelpers::Norm
Point Norm(const Point &a)
Definition: JetVoronoiDiagramHelpers.cxx:79
AlMat.h
Trk::AlMat
Definition: AlMat.h:27
hist_file_dump.f
f
Definition: hist_file_dump.py:140
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Trk::AlMat::nrow
long int nrow() const
Definition: AlMat.h:157
CxxUtils::set
constexpr std::enable_if_t< is_bitmask_v< E >, E & > set(E &lhs, E rhs)
Convenience function to set bits in a class enum bitmask.
Definition: bitmask.h:232
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::AlMat::Write
StatusCode Write(const std::string &, bool, unsigned int precision=6)
Definition: AlMat.cxx:414
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
Trk::AlMat::m_nele
int m_nele
Definition: AlMat.h:111
Trk::AlMat::elemc
double elemc(long int, long int) const
Definition: AlMat.cxx:109
Trk::AlMat::~AlMat
virtual ~AlMat()
Definition: AlMat.cxx:67
Trk::AlMat::m_ptr_data
double * m_ptr_data
Definition: AlMat.h:114
Trk::AlMat::elemr
double & elemr(long int, long int)
Definition: AlMat.cxx:97
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
Trk::AlMat::Normal
AlMat & Normal()
Definition: AlMat.cxx:316
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:23
extractSporadic.q
list q
Definition: extractSporadic.py:97
Trk::AlMat::operator-
AlMat operator-(const AlMat &m) const
Definition: AlMat.cxx:200
Trk::AlMat::operator*=
AlMat & operator*=(const double &d)
Definition: AlMat.cxx:288
Trk::AlMat::elem
virtual long int elem(long int, long int) const
Definition: AlMat.cxx:121
plotBeamSpotMon.nc
int nc
Definition: plotBeamSpotMon.py:82
Trk::v
@ v
Definition: ParamDefs.h:78
Trk::AlMat::m_pathtxt
std::string m_pathtxt
Definition: AlMat.h:113
fitman.k
k
Definition: fitman.py:528
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106