ATLAS Offline Software
AlSymMat.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 
10 #include "TrkAlgebraUtils/AlMat.h"
11 
12 #include <algorithm>
13 #include <cfloat> //for DBL_EPSILON
14 #include <cmath>
15 #include <cstdint>
16 #include <fstream>
17 #include <iomanip>
18 #include <iostream>
19 #include <stdexcept>
20 
21 #include <TMatrixDSparse.h>
22 
23 
24 extern "C" {
25  void dsptrf_(char*, const int*, double*, int*, int* );
26 }
27 
28 extern "C" {
29  void dsptri_(char*, const int*, double*, int*, double*, int* );
30 }
31 
32 extern "C" {
33  void dspev_(char*, char*, const int*, double*,
34  double*, double*, const int*, double*, int* );
35 }
36 
37 namespace Trk {
38 
39 //______________________________________________________________________________
41  : m_ptr_data(nullptr)
42  , m_pathbin("./")
43  , m_pathtxt("./")
44 {
45  m_matrix_type = 1;
46 }
47 
48 
49 //______________________________________________________________________________
51  : AlSymMatBase(N)
52  , m_ptr_data(new double[m_nele])
53  , m_pathbin("./")
54  , m_pathtxt("./")
55 {
56  m_matrix_type = 1;
57 
58  double* p = m_ptr_data + m_nele;
59  while (p > m_ptr_data) *(--p) = 0.;
60 
61 }
62 
63 //______________________________________________________________________________
65  : AlSymMatBase(m)
66  , m_ptr_data(new double[m_nele])
67  , m_pathbin(m.m_pathbin)
68  , m_pathtxt(m.m_pathtxt)
69 {
70  m_matrix_type = 1;
71  m_size = m.size();
72  m_nele = m.m_nele;
73 
74  copy(m);
75 }
76 
77 //______________________________________________________________________________
79  : AlSymMatBase (m.size()),
80  m_ptr_data (new double[m_nele]),
81  m_pathbin (m.pathBin()),
82  m_pathtxt (m.pathTxt())
83 {
85 }
86 
87 //______________________________________________________________________________
89 {
90  m_matrix_type = 1;
91  m_size = m.size();
92  m_nele = m_size*(m_size+1)/2;
93  m_pathbin = m.pathBin();
94  m_pathtxt = m.pathTxt();
95 
96  m_ptr_data = new double[m_nele];
97  copy(m);
98  return *this;
99 }
100 
101 //______________________________________________________________________________
103 {
104  if( m_ptr_data != nullptr ) delete [] m_ptr_data;
105  //ptr_map.clear();
106 }
107 
108 //______________________________________________________________________________
110 {
111  // this one implements the fast copy alg.
112  if( size() != m.size()) {
113  throw std::range_error( "AlSymMat::copy: size do not match!" );
114  }
115 
116  double * p = m_ptr_data + m_nele;
117  double* r = m.m_ptr_data + m_nele;
118  while (p > m_ptr_data) *(--p) = *(--r);
119 }
120 
121 //______________________________________________________________________________
123 {
124  //
125  if( size() != m.size()) {
126  throw std::range_error( "AlSymMat::copy: size do not match!" );
127  }
128 
129  (*this) = 0.;
130  //Convert Storage System
131  long int i, j;
132  for (const datamap::value_type& p : *m.ptrMap()) {
133  m.elem(p.first, i, j);
134  AlSymMat::elemr(i,j) = p.second;
135  }
136 
137  }
138 
139 //______________________________________________________________________________
140 void AlSymMat::copy(const AlMat& m)
141 {
142  // copy the lower triangle only!
143  int si = size();
144  if( si != m.nrow() || si != m.ncol() ) {
145  throw std::range_error("AlSymMat::copy: sizes do not match!" );
146  }
147 
148  for( int i=0; i<si; i++ )
149  for( int j=0; j<=i; j++ )
150  elemr(i, j) = m.elemc(i, j);
151 }
152 
153 //______________________________________________________________________________
155 {
156  if (this==&m)
157  return *this;
158 
159  if ( m_ptr_data != m.m_ptr_data ) {
160  reSize(m.size());
161  copy(m);
162  }
163  return *this;
164 }
165 
166 //______________________________________________________________________________
168 {
169  if( m.nrow() != m.ncol() ) {
170  throw std::range_error( "AlSymMat::operator=: a squared matrix is required!" );
171  }
172 
173  reSize(m.nrow());
174  copy(m);
175 
176  return *this;
177 }
178 
179 //______________________________________________________________________________
181 {
182  double * p = m_ptr_data + m_nele;
183  while (p > m_ptr_data) *(--p) = d;
184 
185  return *this;
186 }
187 
188 //______________________________________________________________________________
190 {
191  if( size() != m.size()) {
192  throw std::range_error( "AlSymMat::operator+: size do not match!" );
193  }
194 
195  AlSymMat b(size());
196  double * p = m_ptr_data + m_nele;
197  double * q = m.m_ptr_data + m_nele;
198  double * r = b.m_ptr_data + m_nele;
199  while (p > m_ptr_data) *(--r) = (*(--p))+(*(--q));
200 
201  return b;
202 }
203 
204 //______________________________________________________________________________
206 {
207  if( size() != m.size()){
208  throw std::range_error( "AlSymMat::operator+=: size do not match!" );
209  }
210 
211  double * p = m_ptr_data + m_nele;
212  double * q = m.m_ptr_data + m_nele;
213  while (p > m_ptr_data) *(--p) += *(--q);
214 
215  return *this;
216 }
217 
218 //______________________________________________________________________________
220 {
221  if( size() != m.size()) {
222  throw std::range_error( "AlSymMat::operator-: size do not match!" );
223  }
224 
225  AlSymMat b(size());
226  double * p = m_ptr_data + m_nele;
227  double * q = m.m_ptr_data + m_nele;
228  double * r = b.m_ptr_data + m_nele;
229  while (p > m_ptr_data) *(--r) = (*(--p))-(*(--q));
230 
231  return b;
232 }
233 
234 //______________________________________________________________________________
236 {
237  if( size() != m.size()) {
238  throw std::range_error( "AlSymMat::operator-=: size do not match!" );
239  }
240 
241  double * p = m_ptr_data + m_nele;
242  double * q = m.m_ptr_data + m_nele;
243  while (p > m_ptr_data) *(--p) -= *(--q);
244 
245  return *this;
246 }
247 
248 //______________________________________________________________________________
250 {
251  if( size() != m.size() ) {
252  throw std::range_error( "AlSymMat::operator*: size do not match!" );
253  }
254 
255  AlMat b( size(), size());
256  double x;
257  long int k, l, n, ii, jj;
258  long int siz(size());
259 
260  // optimised version of symmetric matrix product:
261  for( long int i=0; i<siz; i++ )
262  for( long int j=0; j<siz; j++ ) {
263  x = 0.;
264  l = i < j ? i : j;
265  n = i < j ? j : i;
266  ii = (i+1)*i/2;
267  jj = (j+1)*j/2;
268  for( k=0; k<l; k++)
269  x += (*(m_ptr_data+ii+k))*(*(m.m_ptr_data+jj+k));
270  if( i<j ) {
271  for( k=l; k<n; k++)
272  x += (*(m_ptr_data+(k+1)*k/2+i))*(*(m.m_ptr_data+jj+k));
273  }
274  else {
275  for( k=l; k<n; k++)
276  x += (*(m_ptr_data+ii+k))*(*(m.m_ptr_data+(k+1)*k/2+j));
277  }
278  for( k=n; k<siz; k++)
279  x += (*(m_ptr_data+(k+1)*k/2+i))*(*(m.m_ptr_data+(k+1)*k/2+j));
280  b.elemr(i,j) = x;
281  }
282 
283  return b;
284 }
285 
286 //______________________________________________________________________________
288 {
289  if( size() != m.nrow() ) {
290  throw std::range_error("AlSymMat::operator*: size do not match!" );
291  }
292 
293  AlMat b( size(), m.ncol());
294 
295  double x;
296  long int k, ii;
297  long int siz(size()), nco(m.ncol());
298 
299  // optimised version of matrix product:
300  for( long int i=0; i<siz; i++ )
301  for( long int j=0; j<nco; j++ ) {
302  x = 0.;
303  ii = (i+1)*i/2;
304  for( k=0; k<i; k++)
305  x += (*(m_ptr_data+ii+k))*m.elemc(k,j);
306  for( k=i; k<siz; k++)
307  x += (*(m_ptr_data+(k+1)*k/2+i))*m.elemc(k,j);
308  b.elemr(i,j) = x;
309  }
310 
311  return b;
312 }
313 
314 //______________________________________________________________________________
316 {
317  if( size() != v.size() ) {
318  throw std::range_error( "AlSymMat::operator*: size do not match! " );
319  }
320 
321  AlVec b(size());
322 
323  long int ii;
324 
325  for( long int i=0; i<size(); i++ ) {
326  ii = (i+1)*i/2;
327  for( long int j=0; j<i; j++ )
328  b[i] += (*(m_ptr_data+ii+j))*v[j];
329  for( long int j=i; j<size(); j++ )
330  b[i] += (*(m_ptr_data+(j+1)*j/2+i))*v[j];
331  }
332 
333  return b;
334 }
335 
336 //______________________________________________________________________________
338 {
339  double * p = m_ptr_data + m_nele;
340  while (p > m_ptr_data) *(--p) *= d;
341 
342  return *this;
343 }
344 
345 //______________________________________________________________________________
346 AlSymMat AlSymMat::operator*(const double& d) const
347 {
348  AlSymMat a(size());
349  double * p = m_ptr_data + m_nele;
350  double * q = a.m_ptr_data + m_nele;
351  while (p > m_ptr_data) *(--q) = (*(--p))*d;
352 
353  return a;
354 }
355 
356 //______________________________________________________________________________
358 {
359  // add a protection to detect singular matrices:
360  // double det = (*this).determinant();
361  // double detcut = 1.E-20;
362 
363  // if (fabs(det)>detcut) {
364 
365  int ierr = 0;
366  char uplo = 'U';
367  int N = size();
368  int * ipiv = new int[N];
369  double * work = new double[N];
370  double * mx = m_ptr_data;
371 
372  dsptrf_(&uplo, &N, mx, ipiv, &ierr);
373 
374  dsptri_(&uplo, &N, mx, ipiv, work, &ierr);
375 
376  delete [] ipiv;
377  delete [] work;
378 
379  // } else {
380  // ierr=999;
381  // }
382  return ierr;
383 }
384 
385 //______________________________________________________________________________
387 {
388  double deter = 1.;
389 
390  // get determinant using LU decomposition + Crout algorithm
391  AlMat A(size(),size());
392  A = (*this);
393 
394  double AMAX,DUM,SUM;
395  //int CODE;
396  double TINY = 1.E-20;
397  double D;
398  int N=size();
399 
400  D=1.; //CODE=0;
401 
402  for(int I=0;I<N;I++) {
403  AMAX=0.;
404  for(int J=0;J<N;J++)
405  if (fabs(A[I][J])>AMAX)
406  AMAX=fabs(A[I][J]);
407  if(AMAX<TINY) {
408 // CODE = 1; // matrix is singular
409  return 0;
410  }
411  }
412 
413  for(int J=0;J<N;J++) {
414  for(int I=0;I<N;I++) {
415  SUM = A[I][J];
416  int KMAX = (I<J)?I:J;
417  for(int K=0;K<KMAX;K++)
418  SUM = SUM - A[I][K]*A[K][J];
419  A[I][J]=SUM;
420  }
421 
422  // find pivot and exchange if necessary
423  int IMAX=J;
424  for (int I=J+1;I<N;I++) {
425  if(fabs(A[I][J])>fabs(A[IMAX][J])) IMAX=I;
426  }
427 
428  if(IMAX!=J) {
429  for(int K=0;K<N;K++) {
430  DUM=A[IMAX][K];
431  A[IMAX][K]=A[J][K];
432  A[J][K]=DUM;
433  }
434  D=-D;
435  }
436 
437  if(J<N && (fabs(A[J][J])>TINY)) {
438  DUM = 1./ A[J][J];
439  for(int I=J+1;I<N;I++) A[I][J] = A[I][J]*DUM;
440  }
441  }
442 
443  deter=deter*D;
444  for(int I=0;I<N;I++)
445  deter = deter*A[I][I];
446 
447  return deter;
448 }
449 
450 //______________________________________________________________________________
451 int AlSymMat::diagonalize(char jobz, AlVec &w, AlMat &z)
452 {
453  int info = 0;
454  char uplo = 'U';
455  int N = size();
456  double * work = new double[3*N];
457 
458  double * ap = m_ptr_data;
459 
460  dspev_(&jobz, &uplo, &N, ap, w.ptrData(), z.ptrData(), &N, work, &info);
461 
462  delete [] work;
463 
464  return info;
465 }
466 
467 //______________________________________________________________________________
468 // index here is module number within matrix (size=6*max_index)
470 {
471  const int shift=6;
472 
473  if((shift*index)<size() && index>=0) {
474  for(int row=(shift*index); row<(size()-shift); row+=shift) {
475  for(int k=0; k<shift; k++) {
476  for(int col=0; col<row+k+1; col++) {
477  // cout << "(" << row+k << "," << col << ") -> " << *(m_ptr_data+elem(row+k,col)) << endl;
478 
479  if(col<shift*index)
480  *(m_ptr_data+elem(row+k,col))=*(m_ptr_data+elem(row+k+shift,col));
481  else
482  *(m_ptr_data+elem(row+k,col))=*(m_ptr_data+elem(row+k+shift,col+shift));
483  }
484  }
485  }
486  }
487 }
488 
489 
490 //______________________________________________________________________________
491 int AlSymMat::RemoveCollsRows(std::vector<int> indices)
492 {
493  int n = indices.size();
494  if (n==0) {
495  return m_size;
496  }
497  if (n>m_size) {
498  throw std::range_error("Vector of indices larger than matrix size.");
499  }
500 
501  // first sort the list of indices descending
502  // maybe not the fastest way to do that but it works
503  for (int i=0;i<n-1;i++)
504  for (int j=i+1; j<n; j++)
505  if (indices[j] > indices[i]) {
506  int itmp = indices[i];
507  indices[i] = indices[j];
508  indices[j] = itmp;
509  }
510 
511  // remove rows and columns starting from largest indices
512  for (int i=0;i<n;i++) {
513  int index = indices[i];
514  if (index > m_size-1) {
515  throw std::out_of_range("AlSymMat::RemoveCollsRows: Index goes beyond matrix.");
516  }
517 
518  for (int j=index; j<m_size-1; j++)
519  for (int col=0; col<=j; col++)
520  if(col<index)
521  *(m_ptr_data+elem(j,col)) = *(m_ptr_data+elem(j+1,col));
522  else
523  *(m_ptr_data+elem(j,col)) = *(m_ptr_data+elem(j+1,col+1));
524 
525  m_size--;
526  }
527 
528  return m_size;
529 }
530 
531 //______________________________________________________________________________
532 // index here is alignment parameter to remove
533 void AlSymMat::RemoveAlignPar(int index, int control)
534 {
535  // Dynamic shiftRow
536  int shiftRow = 1;
537  int counterRow = 0;
538 
539  // Dynamic shiftCol
540  int shiftCol;
541  int counterCol = 0;
542 
543  index = index-control;
544 
545  for(int row=index; row<(size()-shiftRow); row++) {
546 
547  shiftCol = 0;
548  counterCol = 0;
549 
550  for(int col=0; col<row+1; col++) {
551 
552  if (col==index) {
553  shiftCol = 1;
554  counterCol=0;
555  }
556 
557  *(m_ptr_data+elem(row,col))=*(m_ptr_data+elem(row+shiftRow,col+shiftCol));
558 
559  counterCol++;
560  if (counterCol==5-control) {
561  counterCol=0;
562  shiftCol++;
563  }
564  }
565 
566  counterRow++;
567  if (counterRow==5-control) {
568  counterRow=0;
569  shiftRow++;
570  }
571  }
572 }
573 
574 //______________________________________________________________________________
575 void AlSymMat::SetPathBin(const std::string &path)
576 {
577  m_pathbin.assign(path);
578 }
579 
580 //______________________________________________________________________________
581 void AlSymMat::SetPathTxt(const std::string &path)
582 {
583  m_pathtxt.assign(path);
584 }
585 
586 //______________________________________________________________________________
587 StatusCode AlSymMat::Write(const std::string &filename, bool binary,
588  bool square, double, float version)
589 {
590  std::ofstream outmat;
591 
592  int32_t msizz = square ? size() : (-1)*size();
593 
594  if(binary) {
595  outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
596  if(outmat.fail())
597  return StatusCode::FAILURE;
598  // if triangular, change size sign to distinguish
599  // from square format
600  outmat.write((char*)&msizz, sizeof (msizz));
601  outmat.write((char*)&version, sizeof (version));
602  }
603  else {
604  outmat.open((m_pathtxt+filename).c_str());
605  if(outmat.fail())
606  return StatusCode::FAILURE;
607  outmat.setf(std::ios::fixed);
608  outmat.setf(std::ios::showpoint);
609  outmat.precision(3);
610  outmat << "msizz: " << std::setw(6) << msizz << std::endl;
611  outmat << "AlSymMat version: " << std::setw(6) << version << std::endl;
612  }
613 
614  double melem=0;
615 
616  if (!square) {
617  for( int i=0; i<size(); i++) {
618  for( int j=0; j<=i; j++) {
619  melem = *(m_ptr_data+(i+1)*i/2+j);
620  if(binary)
621  outmat.write((char*)&(melem), sizeof (melem));
622  else
623  outmat << std::setw(14) << melem;
624  }
625  if(!binary)
626  outmat << std::endl;
627  }
628  }
629  else {
630  for( int i=0; i<size(); i++) {
631  for( int j=0; j<size(); j++) {
632  if(i>=j) melem = *(m_ptr_data+(i+1)*i/2+j);
633  else melem = *(m_ptr_data+(j+1)*j/2+i);
634  if(binary)
635  outmat.write((char*)&(melem), sizeof (melem));
636  else
637  outmat << std::setw(14) << melem;
638  }
639  if(!binary)
640  outmat << std::endl;
641  }
642  }
643  outmat.close();
644  return StatusCode::SUCCESS;
645 }
646 
647 //______________________________________________________________________________
648 StatusCode AlSymMat::CheckMatVersion(const std::string& filename, bool &StdUnits)
649 {
650  std::ifstream inmat((filename).c_str(), std::ios::binary);
651  if(inmat.fail())
652  return StatusCode::FAILURE;
653 
654  int32_t msiz=0;
655  inmat.read((char*)&msiz, sizeof (msiz));
656 
657  float version=0.0;
658  inmat.read((char*)&version, sizeof (version));
659 
660  StdUnits = version>=2.0;
661 
662  inmat.close();
663 
664  // std::cout << "AlSymMat::StdUnits: " << StdUnits << std::endl;
665 
666  return StatusCode::SUCCESS;
667 }
668 
669 //______________________________________________________________________________
670 StatusCode AlSymMat::Read(const std::string &filename, int &dofs,
671  bool &triang, float &version)
672 {
673  bool stdUnits = true;
674  if (StatusCode::SUCCESS != CheckMatVersion(m_pathbin+filename, stdUnits))
675  return StatusCode::FAILURE;
676 
677  // std::cout << "AlSymMat::StdUnits: " << stdUnits << std::endl;
678 
679  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
680  if(inmat.fail())
681  return StatusCode::FAILURE;
682 
683  int32_t msiz=0;
684  inmat.read((char*)&msiz, sizeof (msiz));
685  dofs = abs(msiz);
686  m_size = abs(msiz);
687 
688  if (stdUnits)
689  inmat.read((char*)&version, sizeof (version));
690 
691  double melem=0;
692 
693  if(msiz>0) { // square format
694  triang=false;
695  for(int i=0; i<msiz; i++) {
696  for(int j=0; j<msiz; j++) {
697  inmat.read((char*)&melem, sizeof (melem));
698  if( i>=j )
699  *(m_ptr_data+(i+1)*i/2+j) = melem;
700  }
701  }
702  }
703  else { // triangular format
704  triang=true;
705  msiz = (-1)*msiz;
706  // std::cout << "msiz=" << msiz << std::endl;
707  for( int i=0; i<msiz; i++) {
708  for( int j=0; j<=i; j++) {
709  inmat.read((char*)&melem, sizeof (melem));
710  *(m_ptr_data+(i+1)*i/2+j) = melem;
711  }
712  }
713  }
714 
715  inmat.close();
716  return StatusCode::SUCCESS;
717 }
718 
719 //______________________________________________________________________________
720 StatusCode AlSymMat::ReadProjected(const std::string &filename, int &dofs,
721  bool &triang, float &version)
722 {
723  std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
724  if(inmat.fail())
725  return StatusCode::FAILURE;
726 
727  int32_t msiz=0;
728  inmat.read((char*)&msiz, sizeof (msiz));
729  dofs = abs(msiz);
730  m_size = abs(msiz);
731 
732  inmat.read((char*)&version, sizeof (version));
733 
734  double melem=0;
735 
736  if(msiz>0) { // square format
737  triang=false;
738  for(int i=0; i<msiz; i++) {
739  for(int j=0; j<msiz; j++) {
740  inmat.read((char*)&melem, sizeof (melem));
741  if( i>=j )
742  *(m_ptr_data+(i+1)*i/2+j) = melem;
743  }
744  }
745  }
746  else { // triangular format
747  triang=true;
748  msiz = (-1)*msiz;
749  // std::cout << "msiz=" << msiz << std::endl;
750  for( int i=0; i<msiz; i++) {
751  for( int j=0; j<=i; j++) {
752  inmat.read((char*)&melem, sizeof (melem));
753  *(m_ptr_data+(i+1)*i/2+j) = melem;
754  }
755  }
756  }
757 
758  inmat.close();
759  return StatusCode::SUCCESS;
760 }
761 
762 //______________________________________________________________________________
763 void AlSymMat::reSize(long int Nnew)
764 {
765  if ( Nnew != size() ) {
766  double * p = m_ptr_data;
767  long int ii;
768 
769  m_nele = Nnew*(Nnew+1)/2;
770  m_ptr_data = new double[m_nele];
771  long int Size_old = size();
772  m_size = Nnew;
773  long int k = size() <= Size_old ? size() : Size_old;
774  for( long int i=0; i<k; i++ ) {
775  ii = (i+1)*i/2;
776  for( long int j=0; j<=i; j++ )
777  *(m_ptr_data+ii+j) = *(p+ii+j);
778  }
779  delete [] p;
780  }
781 }
782 
783 //______________________________________________________________________________
784 double& AlSymMat::elemr(long int i,long int j)
785 {
786 #ifdef _DEBUG
787  if( i<0 ) {
788  throw std::underflow_error( "AlSymMat::elemr: Index 1 < zero!" );
789  }
790  if( i>=size() ) {
791  throw std::overflow_error( "AlSymMat::elemr: Index 1 too large!" );
792  }
793  if( j<0 ) {
794  throw std::underflow_error( "AlSymMat::elemr: Index 2 < zero!" );
795  }
796  if( j>=size() ) {
797  throw std::overflow_error( "AlSymMat::elemr: Index 2 too large!" );
798  }
799 #endif
800  if( j<=i )
801  return *(m_ptr_data+(i+1)*i/2+j);
802 
803  return *(m_ptr_data+(j+1)*j/2+i);
804 }
805 
806 //______________________________________________________________________________
807 double AlSymMat::elemc(long int i,long int j) const
808 {
809 #ifdef _DEBUG
810  if( i<0 ) {
811  throw std::underflow_error( "AlSymMat::elemc: Index 1 < zero!" );
812  }
813  if( i>=size() ) {
814  throw std::overflow_error( "AlSymMat::elemc: Index 1 too large!" );
815  }
816  if( j<0 ) {
817  throw std::underflow_error( "AlSymMat::elemc: Index 2 < zero!" );
818  }
819  if( j>=size() ) {
820  throw std::overflow_error( "AlSymMat::elemc: Index 2 too large!" );
821  }
822 #endif
823  if( j<=i )
824  return *(m_ptr_data+(i+1)*i/2+j);
825 
826  return *(m_ptr_data+(j+1)*j/2+i);
827 }
828 
829 
831 
832  int si = size();
833 
834  int nnz(0);
835 
836  for( int i=0; i<si; i++ )
837  for( int j=0; j<=i; j++ )
838  if( *(m_ptr_data+(i+1)*i/2+j) != 0 )
839  ++nnz;
840 
841 
842  int *irow = new int[nnz];
843  int *icol = new int[nnz];
844  double *val = new double[nnz];
845  int n(0);
846 
847  for (int i=0;i<si;i++) {
848  for (int j=0;j<=i;j++) {
849  if ( *(m_ptr_data+(i+1)*i/2+j) != 0.) {
850  irow[n] = i;
851  icol[n] = j;
852  val[n] = *(m_ptr_data+(i+1)*i/2+j);
853  ++n;
854  }
855  }
856  }
857 
858 
859  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,si-1,0,si-1);
860  myTMatrix->SetMatrixArray(nnz,irow,icol,val);
861 
862  delete [] irow;
863  delete [] icol;
864  delete [] val;
865 
866 
867  return myTMatrix;
868 
869 }
870 
871 
872 } // end namespace Trk
grepfile.info
info
Definition: grepfile.py:38
Trk::AlSymMat::RemoveAlignPar
virtual void RemoveAlignPar(int, int) override final
Definition: AlSymMat.cxx:533
query_example.row
row
Definition: query_example.py:24
beamspotman.r
def r
Definition: beamspotman.py:676
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
Trk::AlSymMat::operator+
virtual AlSymMat operator+(const AlSymMat &) const
Definition: AlSymMat.cxx:189
Trk::AlVec
Definition: AlVec.h:23
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::indices
std::pair< long int, long int > indices
Definition: AlSymMatBase.h:24
Trk::AlSymMat::Write
virtual StatusCode Write(const std::string &, bool, bool, double, float) override final
Definition: AlSymMat.cxx:587
CSV_InDetExporter.new
new
Definition: CSV_InDetExporter.py:145
Trk::AlSymMat::operator+=
virtual AlSymMat & operator+=(const AlSymMat &)
Definition: AlSymMat.cxx:205
AlVec.h
Trk::AlSymMat::m_pathbin
std::string m_pathbin
Definition: AlSymMat.h:89
Trk::AlSymMat::elem
long int elem(long int, long int) const
Definition: AlSymMat.h:96
fitman.mx
mx
Definition: fitman.py:520
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
AlSymMat.h
Trk::AlSymMat::invert
virtual int invert() override final
Definition: AlSymMat.cxx:357
Trk::AlSymMat::SetPathBin
virtual void SetPathBin(const std::string &) override final
Definition: AlSymMat.cxx:575
python.AtlRunQueryParser.ap
ap
Definition: AtlRunQueryParser.py:826
Trk::AlSymMat::operator*=
virtual AlSymMat & operator*=(const double &)
Definition: AlSymMat.cxx:337
atlasStyleMacro.icol
int icol
Definition: atlasStyleMacro.py:13
Trk::AlSymMat
Definition: AlSymMat.h:26
TMatrixDSparse
class TMatrixTSparse< double > TMatrixDSparse
Definition: AlSpaMat.h:14
Trk::AlSymMatBase::size
long int size() const
Definition: AlSymMatBase.h:145
Trk::AlSymMat::SetPathTxt
virtual void SetPathTxt(const std::string &) override final
Definition: AlSymMat.cxx:581
Trk::AlSymMat::operator=
AlSymMat & operator=(const AlSpaMat &)
Definition: AlSymMat.cxx:88
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
lumiFormat.i
int i
Definition: lumiFormat.py:85
dsptrf_
void dsptrf_(char *, const int *, double *, int *, int *)
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::AlSymMat::RemoveModule
virtual void RemoveModule(int) override final
Definition: AlSymMat.cxx:469
Trk::AlSymMat::Read
virtual StatusCode Read(const std::string &, int &, bool &, float &) override final
Definition: AlSymMat.cxx:670
AlMat.h
Trk::AlMat
Definition: AlMat.h:27
Trk::AlSymMatBase::m_nele
long int m_nele
Definition: AlSymMatBase.h:104
Trk::AlSymMat::elemc
virtual double elemc(long int, long int) const override final
Definition: AlSymMat.cxx:807
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Trk::AlSymMat::m_pathtxt
std::string m_pathtxt
Definition: AlSymMat.h:90
Trk::AlSymMat::makeTMatrix
virtual TMatrixDSparse * makeTMatrix() override final
Definition: AlSymMat.cxx:830
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::AlSymMat::copy
void copy(const AlSymMat &m)
Definition: AlSymMat.cxx:109
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
query_example.col
col
Definition: query_example.py:7
dspev_
void dspev_(char *, char *, const int *, double *, double *, double *, const int *, double *, int *)
Trk::AlSymMat::ReadProjected
virtual StatusCode ReadProjected(const std::string &, int &, bool &, float &) override final
Definition: AlSymMat.cxx:720
get_generator_info.version
version
Definition: get_generator_info.py:33
Trk::AlSymMat::diagonalize
virtual int diagonalize(char jobz, AlVec &w, AlMat &z) override final
Definition: AlSymMat.cxx:451
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::AlSymMat::AlSymMat
AlSymMat()
Definition: AlSymMat.cxx:40
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
dsptri_
void dsptri_(char *, const int *, double *, int *, double *, int *)
Trk::AlSymMat::RemoveCollsRows
virtual int RemoveCollsRows(std::vector< int >) override final
Definition: AlSymMat.cxx:491
Trk::AlSymMat::reSize
virtual void reSize(long int Nnew) override final
Definition: AlSymMat.cxx:763
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
extractSporadic.q
list q
Definition: extractSporadic.py:98
Trk::AlSymMat::operator*
virtual AlMat operator*(const AlSymMat &) const
Definition: AlSymMat.cxx:249
Trk::AlSymMatBase::m_matrix_type
int m_matrix_type
Definition: AlSymMatBase.h:100
I
#define I(x, y, z)
Definition: MD5.cxx:116
Trk::AlSymMat::operator-
virtual AlSymMat operator-(const AlSymMat &) const
Definition: AlSymMat.cxx:219
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::x
@ x
Definition: ParamDefs.h:55
Trk::AlSymMatBase::m_size
long int m_size
Definition: AlSymMatBase.h:103
Trk::AlSymMat::m_ptr_data
double * m_ptr_data
Definition: AlSymMat.h:87
Trk::AlSymMat::~AlSymMat
~AlSymMat()
Definition: AlSymMat.cxx:102
Trk::AlSymMat::determinant
virtual double determinant() override final
Definition: AlSymMat.cxx:386
Trk::AlSymMat::CheckMatVersion
virtual StatusCode CheckMatVersion(const std::string &, bool &) override final
Definition: AlSymMat.cxx:648
Trk::AlSpaMat
Definition: AlSpaMat.h:27
Trk::v
@ v
Definition: ParamDefs.h:78
Trk::AlSymMat::elemr
virtual double & elemr(long int, long int) override final
Definition: AlSymMat.cxx:784
fitman.k
k
Definition: fitman.py:528
Trk::AlSymMat::operator-=
virtual AlSymMat & operator-=(const AlSymMat &)
Definition: AlSymMat.cxx:235