FEDRA emulsion software from the OPERA Collaboration
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MATRIX::VtSymMatrix Class Reference

#include <VtSymMatrix.hh>

Inheritance diagram for MATRIX::VtSymMatrix:
Collaboration diagram for MATRIX::VtSymMatrix:

Public Member Functions

— Constructors —
 VtSymMatrix (const unsigned int dim)
 
 VtSymMatrix (const unsigned int dim, double diag)
 initalize diagonal elements More...
 
 VtSymMatrix (const VtSymMatrix &rhs)
 copy constructor More...
 
 VtSymMatrix (const VtSqMatrix &rhs)
 
 VtSymMatrix (const VtMatrix &rhs)
 
virtual ~VtSymMatrix ()
 
— Matrix operations —
bool invert (const bool use_momentum)
 calc inverse using momentum or not More...
 
bool VtDsinv (int dim=0)
 transform to inverse More...
 
const VtSymMatrix dsinv (int dim=0) const
 return inverse More...
 
virtual double det () const
 compute determinant via CERNLIB dsfact() More...
 
double product (const VtVector &rhs, unsigned int dim=0) const
 compute $v^t*A*v$ More...
 
const VtSymMatrix product (const VtMatrix &rhs) const
 compute $B^t*A*B$ More...
 
virtual void place_at (const VtMatrix &rhs, const unsigned int row, const unsigned int col)
 copy a smaller matrix at a certain place More...
 
void copy (const VtSymMatrix &rhs)
 in case matrix dimensions differ More...
 
const VtSymMatrixoperator= (const VtSymMatrix &rhs)
 copy More...
 
const VtSymMatrixoperator+= (const double rhs)
 $\textbf{A} = (a_{\mu\nu} + \alpha)$ More...
 
const VtSymMatrixoperator-= (const double rhs)
 $\textbf{A} = (a_{\mu\nu} - \alpha)$ More...
 
const VtSymMatrixoperator*= (const double rhs)
 $\textbf{A} = (a_{\mu\nu} \cdot\alpha)$ More...
 
const VtSymMatrixoperator/= (const double rhs)
 $\textbf{A} = (a_{\mu\nu} / \alpha)$ More...
 
const VtSymMatrixoperator+= (const VtSymMatrix &rhs)
 $\textbf{A} = (a_{\mu\nu} + b_{\mu\nu})$ More...
 
const VtSymMatrixoperator-= (const VtSymMatrix &rhs)
 $\textbf{A} = (a_{\mu\nu} - b_{\mu\nu})$ More...
 
const VtSymMatrix operator+ (const VtSymMatrix &rhs) const
 
const VtSqMatrix operator+ (const VtSqMatrix &rhs) const
 
const VtMatrix operator+ (const VtMatrix &rhs) const
 
const VtSymMatrix operator- (const VtSymMatrix &rhs) const
 
const VtSymMatrix operator- (void) const
 
const VtSqMatrix operator- (const VtSqMatrix &rhs) const
 
const VtMatrix operator- (const VtMatrix &rhs) const
 
const VtSqMatrix operator* (const VtSymMatrix &rhs) const
 
const VtSqMatrix operator* (const VtSqMatrix &rhs) const
 
const VtMatrix operator* (const VtMatrix &rhs) const
 
const VtVector operator* (const VtVector &rhs) const
 
- Public Member Functions inherited from MATRIX::VtSqMatrix
 VtSqMatrix (const int row)
 
 VtSqMatrix (const VtSqMatrix &rhs)
 
 VtSqMatrix (const VtMatrix &rhs)
 
virtual ~VtSqMatrix ()
 
bool VtDinv (int dim=0)
 transform to inverse More...
 
const VtSqMatrix dinv (int dim=0) const
 return inverse More...
 
virtual const VtMatrix operator* (const VtMatrix &rhs)
 
const VtVector operator* (const VtVector &rhs) const
 
- Public Member Functions inherited from MATRIX::VtMatrix
 VtMatrix (const unsigned int row, const unsigned int col)
 
 VtMatrix (const VtMatrix &rhs)
 
virtual ~VtMatrix ()
 
unsigned int nrow () const
 no of rows $n$ More...
 
unsigned int ncol () const
 no of columns $m$ More...
 
int size () const
 $m\times n$ More...
 
VtMatrix_row operator[] (int row)
 
VtMatrix_row_const operator[] (int row) const
 
virtual double operator() (unsigned int row, unsigned int col) const
 
virtual double & operator() (const unsigned int row, const unsigned int col)
 
double get (unsigned int row, unsigned int col) const
 
double & get (unsigned int row, unsigned int col)
 
void VtT (void)
 transform into transpose matrix More...
 
const VtMatrix T (void) const
 return transpose More...
 
virtual void place_at (const VtVector &rhs, const unsigned int row, const unsigned int col)
 copy a vector at a certain place More...
 
void copy (const VtMatrix &rhs)
 to be used if matrix dimensions are not equal More...
 
void clear (void)
 set matrix elements to 0 More...
 
const VtMatrixoperator= (const VtMatrix &rhs)
 $\textbf{A} = \textbf{B}$ More...
 
const VtMatrixoperator= (const VtNegMatrix &rhs)
 
const VtMatrixoperator+= (const VtMatrix &rhs)
 $\textbf{A} = (a_{\mu\nu} + b_{\mu\nu})$ More...
 
const VtMatrixoperator-= (const VtMatrix &rhs)
 $\textbf{A} = (a_{\mu\nu} - b_{\mu\nu})$ More...
 
const VtMatrix operator+ (const VtMatrix &rhs) const
 $\textbf{A} + \textbf{B}$ More...
 
const VtMatrix operator+ (const VtNegMatrix &rhs) const
 
const VtMatrix operator- (const VtMatrix &rhs) const
 $\textbf{A} - \textbf{B}$ More...
 
const VtMatrix operator- (const VtNegMatrix &rhs) const
 
const VtNegMatrix operator- (void) const
 $-\textbf{A}$ More...
 
const VtMatrix operator* (const VtMatrix &rhs) const
 $\textbf{A}\cdot\textbf{B} = \sum_{\nu=1}^n a_{\mu\nu}b_{\nu\lambda}$ More...
 
const VtVector operator* (const VtVector &rhs) const
 $\textbf{A}\cdot\vec{v} = (\sum_{\nu=1}^n a_{\mu\nu}v_{\nu})$ More...
 
double * array () const
 return pointer to internal array More...
 

— Expert functions —

virtual void print (std::ostream &os) const
 
void VtAssert (void)
 

Additional Inherited Members

- Protected Attributes inherited from MATRIX::VtMatrix
double * m
 
double * work
 
unsigned int m_nrow
 
unsigned int m_ncol
 

Detailed Description

Class for symmetric matrices

Constructor & Destructor Documentation

◆ VtSymMatrix() [1/5]

MATRIX::VtSymMatrix::VtSymMatrix ( const unsigned int  dim)
inline
53: VtSqMatrix(dim) {}
VtSqMatrix(const int row)
Definition: VtSqMatrix.hh:54

◆ VtSymMatrix() [2/5]

MATRIX::VtSymMatrix::VtSymMatrix ( const unsigned int  dim,
double  diag 
)

initalize diagonal elements

189 : VtSqMatrix(dim) {
190 for(unsigned int i=0; i<nrow(); ++i)
191 operator()(i,i) = diag;
192 }
unsigned int nrow() const
no of rows $n$
Definition: VtMatrix.hh:61

◆ VtSymMatrix() [3/5]

MATRIX::VtSymMatrix::VtSymMatrix ( const VtSymMatrix rhs)

copy constructor

197 : VtSqMatrix(rhs.nrow()) {
198 const unsigned int nrow = rhs.nrow();
199
200 for(unsigned int i=0; i<nrow; ++i) {
201 // operator()(i,i) = rhs(i,i);
202 get(i,i) = rhs.get(i,i);
203 for(unsigned int j=i+1; j<nrow; ++j) {
204 // operator()(i,j) = operator()(j,i) = rhs(i,j);
205 get(i,j) = get(j,i) = rhs.get(i,j);
206 }
207 }
208 }
double get(unsigned int row, unsigned int col) const
Definition: VtMatrix.hh:105

◆ VtSymMatrix() [4/5]

MATRIX::VtSymMatrix::VtSymMatrix ( const VtSqMatrix rhs)
inline
59 : VtSqMatrix(rhs) {
60 // Assert();
61 }

◆ VtSymMatrix() [5/5]

MATRIX::VtSymMatrix::VtSymMatrix ( const VtMatrix rhs)
inline
63 : VtSqMatrix(rhs) {
64 // Assert();
65 }

◆ ~VtSymMatrix()

MATRIX::VtSymMatrix::~VtSymMatrix ( )
virtual
210{}

Member Function Documentation

◆ copy()

void MATRIX::VtSymMatrix::copy ( const VtSymMatrix rhs)

in case matrix dimensions differ

649 {
650
651 const unsigned int nrow = (m_nrow < rhs.nrow()) ? m_nrow : rhs.nrow();
652
653 for(unsigned int i=0; i<nrow; ++i) {
654 operator()(i,i) = rhs(i,i);
655 for(unsigned int j=i+1; j<nrow; ++j) {
656 operator()(i,j) = operator()(j,i) = rhs(i,j);
657 }
658 }
659
660 return;
661 }
unsigned int m_nrow
Definition: VtMatrix.hh:173
virtual double operator()(unsigned int row, unsigned int col) const
Definition: VtMatrix.C:109

◆ det()

double MATRIX::VtSymMatrix::det ( void  ) const
virtual

compute determinant via CERNLIB dsfact()

Reimplemented from MATRIX::VtSqMatrix.

329 {
330
331 int ifail, jfail;
332 int n = nrow();
333 double det;
334
335 memcpy( work, m, size() * sizeof(double) );
336
337 Dsfact1(&n, work, &n, &ifail, &det, &jfail);
338 if(jfail != 0) {
339 det = 0;
340#ifdef VtDEBUG
341 cout << "VtSymMatrix::det() error: Determinant computation failed!" << endl;
342#endif
343 }
344
345 return det;
346 }
void Dsfact1(int *idim, double *a, int *n, int *ifail, double *det, int *jfail)
Definition: VtSymMatrix.C:43
double * m
Definition: VtMatrix.hh:171
int size() const
$m\times n$
Definition: VtMatrix.hh:65
double * work
Definition: VtMatrix.hh:172
virtual double det() const
compute determinant via CERNLIB dsfact()
Definition: VtSymMatrix.C:329

◆ dsinv()

const VtSymMatrix MATRIX::VtSymMatrix::dsinv ( int  dim = 0) const

return inverse

275 {
276 if(dim==0) {
277 VtSymMatrix tmp(*this);
278 tmp.VtDsinv(dim);
279 return tmp;
280 } else {
281 VtSymMatrix tmp(dim);
282 tmp.copy(*this);
283 tmp.VtDsinv();
284 return tmp;
285 }
286 // return VtSymMatrix(*this).VtDsinv(dim);
287 }
VtSymMatrix(const unsigned int dim)
Definition: VtSymMatrix.hh:53

◆ invert()

bool MATRIX::VtSymMatrix::invert ( const bool  use_momentum)

calc inverse using momentum or not

265 {
266 if(use_momentum)
267 return VtDsinv(nrow());
268 else
269 return VtDsinv(nrow() - 1);
270 }
bool VtDsinv(int dim=0)
transform to inverse
Definition: VtSymMatrix.C:292

◆ operator*() [1/4]

const VtMatrix MATRIX::VtSymMatrix::operator* ( const VtMatrix rhs) const
632 {
633
634 return VtMatrix::operator*(rhs);
635 }
const VtMatrix operator*(const VtMatrix &rhs) const
$\textbf{A}\cdot\textbf{B} = \sum_{\nu=1}^n a_{\mu\nu}b_{\nu\lambda}$
Definition: VtMatrix.C:283

◆ operator*() [2/4]

const VtSqMatrix MATRIX::VtSymMatrix::operator* ( const VtSqMatrix rhs) const
624 {
625
626 return VtMatrix::operator*(rhs);
627 }

◆ operator*() [3/4]

const VtSqMatrix MATRIX::VtSymMatrix::operator* ( const VtSymMatrix rhs) const
616 {
617
618 return VtMatrix::operator*(rhs);
619 }

◆ operator*() [4/4]

const VtVector MATRIX::VtSymMatrix::operator* ( const VtVector rhs) const
640 {
641
642 return VtMatrix::operator*(rhs);
643 }

◆ operator*=()

const VtSymMatrix & MATRIX::VtSymMatrix::operator*= ( const double  rhs)
virtual

$\textbf{A} = (a_{\mu\nu} \cdot\alpha)$

Reimplemented from MATRIX::VtMatrix.

571 {
572 const unsigned int nrow = ncol();
573
574 for(unsigned int i=0; i<nrow; ++i){
575 for(unsigned int j=i; j<nrow; ++j) {
576 operator()(i,j) = operator()(j,i) *= rhs;
577 }
578 }
579 return *this;
580 }
unsigned int ncol() const
no of columns $m$
Definition: VtMatrix.hh:63

◆ operator+() [1/3]

const VtMatrix MATRIX::VtSymMatrix::operator+ ( const VtMatrix rhs) const
443 {
444
445 return VtMatrix::operator+(rhs);
446 }
const VtMatrix operator+(const VtMatrix &rhs) const
$\textbf{A} + \textbf{B}$
Definition: VtMatrix.C:216

◆ operator+() [2/3]

const VtSqMatrix MATRIX::VtSymMatrix::operator+ ( const VtSqMatrix rhs) const
435 {
436
437 return VtMatrix::operator+(rhs);
438 }

◆ operator+() [3/3]

const VtSymMatrix MATRIX::VtSymMatrix::operator+ ( const VtSymMatrix rhs) const
505 {
506
507#ifndef VtFAST
508 assert(ncol() == rhs.ncol());
509#endif
510
511#ifdef VtDEBUG
512 cout << "VtSymMatrix::operator+ called!" << endl;
513#endif
514
515 VtSymMatrix tmp(*this);
516 return tmp += rhs;
517 }

◆ operator+=() [1/2]

const VtSymMatrix & MATRIX::VtSymMatrix::operator+= ( const double  rhs)
virtual

$\textbf{A} = (a_{\mu\nu} + \alpha)$

Reimplemented from MATRIX::VtMatrix.

472 {
473 const unsigned int nrow = ncol();
474
475 for(unsigned int i=0; i<nrow; ++i){
476 for(unsigned int j=i; j<nrow; ++j) {
477 operator()(i,j) = operator()(j,i) += rhs;
478 }
479 }
480 return *this;
481 }

◆ operator+=() [2/2]

const VtSymMatrix & MATRIX::VtSymMatrix::operator+= ( const VtSymMatrix rhs)

$\textbf{A} = (a_{\mu\nu} + b_{\mu\nu})$

486 {
487
488#ifndef VtFAST
489 assert(ncol() == rhs.ncol());
490#endif
491
492 const unsigned int nrow = rhs.nrow();
493
494 for(unsigned int i=0; i<nrow; ++i)
495 for(unsigned int j=i; j<nrow; ++j)
496 operator()(i,j) = operator()(j,i) += rhs(i,j);
497
498 return *this;
499 }

◆ operator-() [1/4]

const VtMatrix MATRIX::VtSymMatrix::operator- ( const VtMatrix rhs) const
607 {
608
609 return VtMatrix::operator-(rhs);
610 }
const VtNegMatrix operator-(void) const
$-\textbf{A}$
Definition: VtMatrix.C:271

◆ operator-() [2/4]

const VtSqMatrix MATRIX::VtSymMatrix::operator- ( const VtSqMatrix rhs) const
599 {
600
601 return VtMatrix::operator-(rhs);
602 }

◆ operator-() [3/4]

const VtSymMatrix MATRIX::VtSymMatrix::operator- ( const VtSymMatrix rhs) const
554 {
555
556#ifndef VtFAST
557 assert(ncol() == rhs.ncol());
558#endif
559
560#ifdef VtDEBUG
561 cout << "VtSymMatrix::operator- called!" << endl;
562#endif
563
564 VtSymMatrix tmp(*this);
565 return tmp -= rhs;
566 }

◆ operator-() [4/4]

const VtSymMatrix MATRIX::VtSymMatrix::operator- ( void  ) const
451 {
452
453#ifdef VtDEBUG
454 cout << "VtSymMatrix::operator- (unary) called!" << endl;
455#endif
456
457 const unsigned int nrow = ncol();
458 VtSymMatrix tmp(nrow);
459
460 for(unsigned int i=0; i<nrow; ++i) {
461 for(unsigned int j=i; j<nrow; ++j) {
462 tmp(i,j) = tmp(j,i) = -get(i,j);
463 }
464 }
465
466 return tmp;
467 }

◆ operator-=() [1/2]

const VtSymMatrix & MATRIX::VtSymMatrix::operator-= ( const double  rhs)
virtual

$\textbf{A} = (a_{\mu\nu} - \alpha)$

Reimplemented from MATRIX::VtMatrix.

522 {
523 const unsigned int nrow = ncol();
524
525 for(unsigned int i=0; i<nrow; ++i){
526 for(unsigned int j=i; j<nrow; ++j) {
527 operator()(i,j) = operator()(j,i) -= rhs;
528 }
529 }
530 return *this;
531 }

◆ operator-=() [2/2]

const VtSymMatrix & MATRIX::VtSymMatrix::operator-= ( const VtSymMatrix rhs)

$\textbf{A} = (a_{\mu\nu} - b_{\mu\nu})$

536 {
537
538#ifndef VtFAST
539 assert(ncol() == rhs.ncol());
540#endif
541
542 const unsigned int nrow = rhs.nrow();
543
544 for(unsigned int i=0; i<nrow; ++i)
545 for(unsigned int j=i; j<nrow; ++j)
546 operator()(i,j) = operator()(j,i) -= rhs(i,j);
547
548 return *this;
549 }

◆ operator/=()

const VtSymMatrix & MATRIX::VtSymMatrix::operator/= ( const double  rhs)
virtual

$\textbf{A} = (a_{\mu\nu} / \alpha)$

Reimplemented from MATRIX::VtMatrix.

585 {
586 const unsigned int nrow = ncol();
587
588 for(unsigned int i=0; i<nrow; ++i){
589 for(unsigned int j=i; j<nrow; ++j) {
590 operator()(i,j) = operator()(j,i) /= rhs;
591 }
592 }
593 return *this;
594 }

◆ operator=()

const VtSymMatrix & MATRIX::VtSymMatrix::operator= ( const VtSymMatrix rhs)

copy

215 {
216 if(this == &rhs) return *this;
217
218#ifndef VtFAST
219 assert(ncol() == rhs.ncol());
220#endif
221
222#ifdef VtDEBUG
223 cout << "use VtSymMatrix::operator=" << endl;
224#endif
225 const unsigned int nrow = ncol();
226
227 for(unsigned int i=0; i<nrow; ++i) {
228 // operator()(i,i) = rhs(i,i);
229 get(i,i) = rhs.get(i,i);
230 for(unsigned int j=i+1; j<nrow; ++j) {
231 // operator()(i,j) = operator()(j,i) = rhs(i,j);
232 get(i,j) = get(j,i) = rhs.get(i,j);
233 }
234 }
235
236 return *this;
237 }

◆ place_at()

void MATRIX::VtSymMatrix::place_at ( const VtMatrix rhs,
const unsigned int  row,
const unsigned int  col 
)
virtual

copy a smaller matrix at a certain place

Reimplemented from MATRIX::VtMatrix.

414 {
415#ifndef VtFAST
416 assert(nrow() >= rhs.nrow()+row && ncol() >= rhs.ncol()+col);
417#endif
418
419 const unsigned int rnrow = rhs.nrow();
420 const unsigned int rncol = rhs.ncol();
421
422 for(unsigned int i=0; i<rnrow; ++i) {
423 const unsigned int nr = i + row;
424 for(unsigned int j=0, nc=col; j<rncol; ++j,++nc) {
425 get(nr,nc) = get(nc,nr) = rhs(i,j);
426 }
427 }
428
429 return;
430 }

◆ print()

void MATRIX::VtSymMatrix::print ( std::ostream &  os) const
virtual

Reimplemented from MATRIX::VtMatrix.

666 {
667 os.setf(ios::scientific,ios::floatfield);
668 VtMatrix::print(os);
669 }
virtual void print(std::ostream &os) const
Definition: VtMatrix.C:358

◆ product() [1/2]

const VtSymMatrix MATRIX::VtSymMatrix::product ( const VtMatrix rhs) const

compute $B^t*A*B$

384 {
385
386#ifndef VtFAST
387 assert(ncol() == rhs.nrow());
388#endif
389
390 const unsigned int rncol = rhs.ncol();
391
392 VtSymMatrix tmp(rncol);
393
394 for(unsigned int i=0; i<rncol; ++i) {
395 for(unsigned int j=i; j<rncol; ++j) {
396 for(unsigned int k=0; k<m_ncol; ++k) {
397 double sum = 0;
398 for(unsigned int l=0; l<m_ncol; ++l)
399 sum += get(k,l) * rhs(l,j);
400 tmp(i,j) += rhs(k,i) * sum;
401 }
402 tmp(j,i) = tmp(i,j);
403 }
404 }
405
406 return tmp;
407 }
unsigned int m_ncol
Definition: VtMatrix.hh:174

◆ product() [2/2]

double MATRIX::VtSymMatrix::product ( const VtVector rhs,
unsigned int  dim = 0 
) const

compute $v^t*A*v$

351 {
352
353#ifndef VtFAST
354 assert(ncol() == rhs.nrow());
355#endif
356
357 if(dim == 0) dim = ncol();
358#ifndef VtFAST
359 else
360 assert(dim <= ncol());
361#endif
362
363 double result = 0.;
364
365 // contribution from diagonal elements
366 for(unsigned int i=0; i<dim; ++i)
367 result += get(i,i) * rhs[i]*rhs[i];
368
369 // make use of fact that Matrix is symmetric
370 for(unsigned int i=0; i<(dim-1); ++i) {
371 double sum = 0.;
372 for(unsigned int j=i+1; j<dim; ++j) {
373 sum += get(i,j) * rhs[j];
374 }
375 result += 2. * rhs[i] * sum;
376 }
377
378 return result;
379 }

◆ VtAssert()

void MATRIX::VtSymMatrix::VtAssert ( void  )
private
242 {
243 for(unsigned int i=0; i<nrow(); ++i)
244 for(unsigned int j=i+1; j<ncol(); ++j) {
245 cout.setf(ios::scientific,ios::floatfield);
246 double ij = operator()(i,j);
247 double ji = operator()(j,i);
248 double dd = 0.;
249 if(fabs(ij) < 1e-20)
250 dd = ij - ji;
251 else
252 dd = (ij - ji)/ij;
253// cout << "i: " << i << " j: " << j
254// << "(i,j): " << ij << " "
255// << "(j,i): " << ji << " "
256// << "diff: " << fabs(dd) << endl;
257 assert( dd < 1.e-6 );
258 }
259 return;
260 }
Expr< UnaryOp< Fabs< T >, Expr< A, T, D >, T >, T, D > fabs(const Expr< A, T, D > &rhs)
Definition: UnaryOperators.hh:96

◆ VtDsinv()

bool MATRIX::VtSymMatrix::VtDsinv ( int  dim = 0)

transform to inverse

292 {
293 bool valid = false;
294
295 if(dim == 0) dim = nrow();
296#ifndef VtFAST
297 else {
298 assert(dim <= static_cast<int>(nrow()));
299 }
300#endif
301
302 const size_t msize = size() * sizeof(double);
303
304 memcpy( work, m, msize );
305
306 int ifail = 0;
307 int isize = nrow();
308 Dsinv1(&dim, work, &isize, &ifail);
309
310 if(ifail != -1) {
311 memcpy( m, work, msize );
312 valid = true;
313 }
314#ifdef VtDEBUG
315 else {
316 cout << "VtSymMatrix::dsinv() error: Matrix inversion failed!!!"
317 << " (nrow=" << nrow() << ",ncol=" << ncol()
318 << ",size=" << size()
319 << ",ifail=" << ifail << ")" << endl;
320 }
321#endif
322
323 return valid;
324 }
void Dsinv1(int *idim, double *a, int *n, int *ifail)
Definition: VtSymMatrix.C:92

The documentation for this class was generated from the following files: