CLHEP 2.4.6.4
C++ Class Library for High Energy Physics
Loading...
Searching...
No Matches
Classes | Public Member Functions | Protected Member Functions | Friends | List of all members
CLHEP::HepMatrix Class Reference

#include <Matrix.h>

+ Inheritance diagram for CLHEP::HepMatrix:

Classes

class  HepMatrix_row
 
class  HepMatrix_row_const
 

Public Member Functions

 HepMatrix ()
 
 HepMatrix (int p, int q)
 
 HepMatrix (int p, int q, int i)
 
 HepMatrix (int p, int q, HepRandom &r)
 
 HepMatrix (const HepMatrix &hm1)
 
 HepMatrix (const HepSymMatrix &)
 
 HepMatrix (const HepDiagMatrix &)
 
 HepMatrix (const HepVector &)
 
virtual ~HepMatrix ()
 
virtual int num_row () const
 
virtual int num_col () const
 
virtual const double & operator() (int row, int col) const
 
virtual double & operator() (int row, int col)
 
HepMatrix & operator*= (double t)
 
HepMatrix & operator/= (double t)
 
HepMatrix & operator+= (const HepMatrix &)
 
HepMatrix & operator+= (const HepSymMatrix &)
 
HepMatrix & operator+= (const HepDiagMatrix &)
 
HepMatrix & operator+= (const HepVector &)
 
HepMatrix & operator-= (const HepMatrix &)
 
HepMatrix & operator-= (const HepSymMatrix &)
 
HepMatrix & operator-= (const HepDiagMatrix &)
 
HepMatrix & operator-= (const HepVector &)
 
HepMatrix & operator= (const HepMatrix &)
 
HepMatrix & operator= (const HepSymMatrix &)
 
HepMatrix & operator= (const HepDiagMatrix &)
 
HepMatrix & operator= (const HepVector &)
 
HepMatrix & operator= (const HepRotation &)
 
HepMatrix operator- () const
 
HepMatrix apply (double(*f)(double, int, int)) const
 
HepMatrix T () const
 
HepMatrix sub (int min_row, int max_row, int min_col, int max_col) const
 
void sub (int row, int col, const HepMatrix &hm1)
 
HepMatrix inverse (int &ierr) const
 
virtual void invert (int &ierr)
 
void invert ()
 
HepMatrix inverse () const
 
double determinant () const
 
double trace () const
 
HepMatrix_row operator[] (int)
 
const HepMatrix_row_const operator[] (int) const
 
- Public Member Functions inherited from CLHEP::HepGenMatrix
virtual ~HepGenMatrix ()
 
virtual int num_row () const =0
 
virtual int num_col () const =0
 
virtual const double & operator() (int row, int col) const =0
 
virtual double & operator() (int row, int col)=0
 
virtual void invert (int &)=0
 
HepGenMatrix_row operator[] (int)
 
const HepGenMatrix_row_const operator[] (int) const
 
virtual bool operator== (const HepGenMatrix &) const
 

Protected Member Functions

virtual int num_size () const
 
virtual void invertHaywood4 (int &ierr)
 
virtual void invertHaywood5 (int &ierr)
 
virtual void invertHaywood6 (int &ierr)
 
- Protected Member Functions inherited from CLHEP::HepGenMatrix
virtual int num_size () const =0
 
void delete_m (int size, double *)
 
double * new_m (int size)
 

Friends

class HepMatrix_row
 
class HepMatrix_row_const
 
class HepVector
 
class HepSymMatrix
 
class HepDiagMatrix
 
void swap (HepMatrix &hm1, HepMatrix &hm2)
 
HepMatrix operator+ (const HepMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator- (const HepMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator* (const HepMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator* (const HepMatrix &hm1, const HepSymMatrix &hm2)
 
HepMatrix operator* (const HepMatrix &hm1, const HepDiagMatrix &hm2)
 
HepMatrix operator* (const HepSymMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator* (const HepDiagMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator* (const HepVector &hm1, const HepMatrix &hm2)
 
HepVector operator* (const HepMatrix &hm1, const HepVector &hm2)
 
HepMatrix operator* (const HepSymMatrix &hm1, const HepSymMatrix &hm2)
 
HepVector solve (const HepMatrix &, const HepVector &)
 
HepVector qr_solve (HepMatrix *, const HepVector &)
 
HepMatrix qr_solve (HepMatrix *, const HepMatrix &b)
 
void tridiagonal (HepSymMatrix *a, HepMatrix *hsm)
 
void row_house (HepMatrix *, const HepMatrix &, double, int, int, int, int)
 
void row_house (HepMatrix *, const HepVector &, double, int, int)
 
void back_solve (const HepMatrix &R, HepVector *b)
 
void back_solve (const HepMatrix &R, HepMatrix *b)
 
void col_givens (HepMatrix *A, double c, double s, int k1, int k2, int rowmin, int rowmax)
 
void row_givens (HepMatrix *A, double c, double s, int k1, int k2, int colmin, int colmax)
 
void col_house (HepMatrix *, const HepMatrix &, double, int, int, int, int)
 
HepVector house (const HepMatrix &a, int row, int col)
 
void house_with_update (HepMatrix *a, int row, int col)
 
void house_with_update (HepMatrix *a, HepMatrix *v, int row, int col)
 
void house_with_update2 (HepSymMatrix *a, HepMatrix *v, int row, int col)
 

Additional Inherited Members

- Public Types inherited from CLHEP::HepGenMatrix
enum  { size_max = 25 }
 
typedef std::vector< double, Alloc< double, 25 > >::iterator mIter
 
typedef std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
 
- Static Public Member Functions inherited from CLHEP::HepGenMatrix
static void swap (int &, int &)
 
static void swap (std::vector< double, Alloc< double, 25 > > &, std::vector< double, Alloc< double, 25 > > &)
 
static void error (const char *s)
 

Detailed Description

Author

Definition at line 205 of file Matrix.h.

Constructor & Destructor Documentation

◆ HepMatrix() [1/8]

CLHEP::HepMatrix::HepMatrix ( )
inline

◆ HepMatrix() [2/8]

CLHEP::HepMatrix::HepMatrix ( int  p,
int  q 
)

Definition at line 61 of file Matrix.cc.

62 : m(p*q), nrow(p), ncol(q)
63{
64 size_ = nrow * ncol;
65}

◆ HepMatrix() [3/8]

CLHEP::HepMatrix::HepMatrix ( int  p,
int  q,
int  i 
)

Definition at line 67 of file Matrix.cc.

68 : m(p*q), nrow(p), ncol(q)
69{
70 size_ = nrow * ncol;
71
72 if (size_ > 0) {
73 switch(init)
74 {
75 case 0:
76 break;
77
78 case 1:
79 {
80 if ( ncol == nrow ) {
81 mIter a = m.begin();
82 for( int step=0; step < size_; step+=(ncol+1) ) *(a+step) = 1.0;
83 } else {
84 error("Invalid dimension in HepMatrix(int,int,1).");
85 }
86 break;
87 }
88 default:
89 error("Matrix: initialization must be either 0 or 1.");
90 }
91 }
92}
CLHEP::HepGenMatrix::mIter
std::vector< double, Alloc< double, 25 > >::iterator mIter
Definition: GenMatrix.h:73
CLHEP::HepGenMatrix::error
static void error(const char *s)
Definition: GenMatrix.cc:70
init
void init()
Definition: testRandom.cc:29

◆ HepMatrix() [4/8]

CLHEP::HepMatrix::HepMatrix ( int  p,
int  q,
HepRandom &  r 
)

Definition at line 94 of file Matrix.cc.

95 : m(p*q), nrow(p), ncol(q)
96{
97 size_ = nrow * ncol;
98
99 mIter a = m.begin();
100 mIter b = m.end();
101 for(; a<b; a++) *a = r();
102}

◆ HepMatrix() [5/8]

CLHEP::HepMatrix::HepMatrix ( const HepMatrix &  hm1)

Definition at line 109 of file Matrix.cc.

110 : HepGenMatrix(hm1), m(hm1.size_), nrow(hm1.nrow), ncol(hm1.ncol), size_(hm1.size_)
111{
112 m = hm1.m;
113
114}

◆ HepMatrix() [6/8]

CLHEP::HepMatrix::HepMatrix ( const HepSymMatrix &  hm1)

Definition at line 145 of file Matrix.cc.

146 : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
147{
148 size_ = nrow * ncol;
149
150 mcIter sjk = hm1.m.begin();
151 // j >= k
152 for(int j=0; j!=nrow; ++j) {
153 for(int k=0; k<=j; ++k) {
154 m[j*ncol+k] = *sjk;
155 // we could copy the diagonal element twice or check
156 // doing the check may be a tiny bit faster,
157 // so we choose that option for now
158 if(k!=j) m[k*nrow+j] = *sjk;
159 ++sjk;
160 }
161 }
162}
CLHEP::HepGenMatrix::mcIter
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
Definition: GenMatrix.h:74

◆ HepMatrix() [7/8]

CLHEP::HepMatrix::HepMatrix ( const HepDiagMatrix &  hm1)

Definition at line 164 of file Matrix.cc.

165 : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
166{
167 size_ = nrow * ncol;
168
169 int n = num_row();
170 mIter mrr;
171 mcIter mr = hm1.m.begin();
172 for(int r=0;r<n;r++) {
173 mrr = m.begin()+(n+1)*r;
174 *mrr = *(mr++);
175 }
176}
CLHEP::HepMatrix::num_row
virtual int num_row() const
Definition: Matrix.cc:118

◆ HepMatrix() [8/8]

CLHEP::HepMatrix::HepMatrix ( const HepVector &  hm1)

Definition at line 178 of file Matrix.cc.

179 : m(hm1.nrow), nrow(hm1.nrow), ncol(1)
180{
181
182 size_ = nrow;
183 m = hm1.m;
184}

◆ ~HepMatrix()

CLHEP::HepMatrix::~HepMatrix ( )
virtual

Definition at line 106 of file Matrix.cc.

106 {
107}

Member Function Documentation

◆ apply()

HepMatrix CLHEP::HepMatrix::apply ( double(*)(double, int, int)  f) const

Definition at line 474 of file Matrix.cc.

477{
478#else
479{
480 HepMatrix mret(num_row(),num_col());
481#endif
482 mcIter a = m.begin();
483 mIter b = mret.m.begin();
484 for(int ir=1;ir<=num_row();ir++) {
485 for(int ic=1;ic<=num_col();ic++) {
486 *(b++) = (*f)(*(a++), ir, ic);
487 }
488 }
489 return mret;
490}
CLHEP::HepMatrix::num_col
virtual int num_col() const
Definition: Matrix.cc:120
f
void f(void g())
Definition: excDblThrow.cc:38

Referenced by main().

◆ determinant()

double CLHEP::HepMatrix::determinant ( ) const

Definition at line 813 of file Matrix.cc.

813 {
814 static CLHEP_THREAD_LOCAL int max_array = 20;
815 static CLHEP_THREAD_LOCAL int *ir = new int [max_array+1];
816 if(ncol != nrow)
817 error("HepMatrix::determinant: Matrix is not NxN");
818 if (ncol > max_array) {
819 delete [] ir;
820 max_array = nrow;
821 ir = new int [max_array+1];
822 }
823 double det;
824 HepMatrix mt(*this);
825 int i = mt.dfact_matrix(det, ir);
826 if(i==0) return det;
827 return 0;
828}

◆ inverse() [1/2]

HepMatrix CLHEP::HepMatrix::inverse ( ) const
inline

◆ inverse() [2/2]

HepMatrix CLHEP::HepMatrix::inverse ( int &  ierr) const
inline

Referenced by main().

◆ invert() [1/2]

void CLHEP::HepMatrix::invert ( )
inline

◆ invert() [2/2]

void CLHEP::HepMatrix::invert ( int &  ierr)
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 705 of file Matrix.cc.

705 {
706 if(ncol != nrow)
707 error("HepMatrix::invert: Matrix is not NxN");
708
709 static CLHEP_THREAD_LOCAL int max_array = 20;
710 static CLHEP_THREAD_LOCAL int *ir = new int [max_array+1];
711
712 if (ncol > max_array) {
713 delete [] ir;
714 max_array = nrow;
715 ir = new int [max_array+1];
716 }
717 double t1, t2, t3;
718 double det, temp, sd;
719 int ifail;
720 switch(nrow) {
721 case 3:
722 double c11,c12,c13,c21,c22,c23,c31,c32,c33;
723 ifail = 0;
724 c11 = (*(m.begin()+4)) * (*(m.begin()+8)) - (*(m.begin()+5)) * (*(m.begin()+7));
725 c12 = (*(m.begin()+5)) * (*(m.begin()+6)) - (*(m.begin()+3)) * (*(m.begin()+8));
726 c13 = (*(m.begin()+3)) * (*(m.begin()+7)) - (*(m.begin()+4)) * (*(m.begin()+6));
727 c21 = (*(m.begin()+7)) * (*(m.begin()+2)) - (*(m.begin()+8)) * (*(m.begin()+1));
728 c22 = (*(m.begin()+8)) * (*m.begin()) - (*(m.begin()+6)) * (*(m.begin()+2));
729 c23 = (*(m.begin()+6)) * (*(m.begin()+1)) - (*(m.begin()+7)) * (*m.begin());
730 c31 = (*(m.begin()+1)) * (*(m.begin()+5)) - (*(m.begin()+2)) * (*(m.begin()+4));
731 c32 = (*(m.begin()+2)) * (*(m.begin()+3)) - (*m.begin()) * (*(m.begin()+5));
732 c33 = (*m.begin()) * (*(m.begin()+4)) - (*(m.begin()+1)) * (*(m.begin()+3));
733 t1 = fabs(*m.begin());
734 t2 = fabs(*(m.begin()+3));
735 t3 = fabs(*(m.begin()+6));
736 if (t1 >= t2) {
737 if (t3 >= t1) {
738 temp = *(m.begin()+6);
739 det = c23*c12-c22*c13;
740 } else {
741 temp = *(m.begin());
742 det = c22*c33-c23*c32;
743 }
744 } else if (t3 >= t2) {
745 temp = *(m.begin()+6);
746 det = c23*c12-c22*c13;
747 } else {
748 temp = *(m.begin()+3);
749 det = c13*c32-c12*c33;
750 }
751 if (det==0) {
752 ierr = 1;
753 return;
754 }
755 {
756 double s1 = temp/det;
757 mIter hmm = m.begin();
758 *(hmm++) = s1*c11;
759 *(hmm++) = s1*c21;
760 *(hmm++) = s1*c31;
761 *(hmm++) = s1*c12;
762 *(hmm++) = s1*c22;
763 *(hmm++) = s1*c32;
764 *(hmm++) = s1*c13;
765 *(hmm++) = s1*c23;
766 *(hmm) = s1*c33;
767 }
768 break;
769 case 2:
770 ifail = 0;
771 det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
772 if (det==0) {
773 ierr = 1;
774 return;
775 }
776 sd = 1.0/det;
777 temp = sd*(*(m.begin()+3));
778 *(m.begin()+1) *= -sd;
779 *(m.begin()+2) *= -sd;
780 *(m.begin()+3) = sd*(*m.begin());
781 *(m.begin()) = temp;
782 break;
783 case 1:
784 ifail = 0;
785 if ((*(m.begin()))==0) {
786 ierr = 1;
787 return;
788 }
789 *(m.begin()) = 1.0/(*(m.begin()));
790 break;
791 case 4:
792 invertHaywood4(ierr);
793 return;
794 case 5:
795 invertHaywood5(ierr);
796 return;
797 case 6:
798 invertHaywood6(ierr);
799 return;
800 default:
801 ifail = dfact_matrix(det, ir);
802 if(ifail) {
803 ierr = 1;
804 return;
805 }
806 dfinv_matrix(ir);
807 break;
808 }
809 ierr = 0;
810 return;
811}
CLHEP::HepMatrix::invertHaywood4
virtual void invertHaywood4(int &ierr)
Definition: MatrixInvert.cc:112
CLHEP::HepMatrix::invertHaywood6
virtual void invertHaywood6(int &ierr)
Definition: MatrixInvert.cc:442
CLHEP::HepMatrix::invertHaywood5
virtual void invertHaywood5(int &ierr)
Definition: MatrixInvert.cc:216

Referenced by test_inversion().

◆ invertHaywood4()

void CLHEP::HepMatrix::invertHaywood4 ( int &  ierr)
protectedvirtual

Definition at line 112 of file MatrixInvert.cc.

112 {
113
114 ifail = 0;
115
116 // Find all NECESSARY 2x2 dets: (18 of them)
117
118 double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
119 double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
120 double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20]; //
121 double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21]; //
122 double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22]; //
123 double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
124 double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
125 double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
126 double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
127 double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
128 double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
129 double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32]; //
130 double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
131 double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
132 double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
133 double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
134 double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
135 double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
136
137 // Find all NECESSARY 3x3 dets: (16 of them)
138
139 double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02
140 + m[F02]*Det2_12_01;
141 double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03
142 + m[F03]*Det2_12_01; //
143 double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03
144 + m[F03]*Det2_12_02; //
145 double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13
146 + m[F03]*Det2_12_12; //
147 double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02
148 + m[F02]*Det2_13_01;
149 double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
150 + m[F03]*Det2_13_01;
151 double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
152 + m[F03]*Det2_13_02; //
153 double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
154 + m[F03]*Det2_13_12; //
155 double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02
156 + m[F02]*Det2_23_01;
157 double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
158 + m[F03]*Det2_23_01;
159 double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
160 + m[F03]*Det2_23_02;
161 double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
162 + m[F03]*Det2_23_12; //
163 double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02
164 + m[F12]*Det2_23_01;
165 double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03
166 + m[F13]*Det2_23_01;
167 double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03
168 + m[F13]*Det2_23_02;
169 double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13
170 + m[F13]*Det2_23_12;
171
172 // Find the 4x4 det:
173
174 double det = m[F00]*Det3_123_123
175 - m[F01]*Det3_123_023
176 + m[F02]*Det3_123_013
177 - m[F03]*Det3_123_012;
178
179 if ( det == 0 ) {
180#ifdef SINGULAR_DIAGNOSTICS
181 std::cerr << "Kramer's rule inversion of a singular 4x4 matrix: "
182 << *this << "\n";
183#endif
184 ifail = 1;
185 return;
186 }
187
188 double oneOverDet = 1.0/det;
189 double mn1OverDet = - oneOverDet;
190
191 m[F00] = Det3_123_123 * oneOverDet;
192 m[F01] = Det3_023_123 * mn1OverDet;
193 m[F02] = Det3_013_123 * oneOverDet;
194 m[F03] = Det3_012_123 * mn1OverDet;
195
196 m[F10] = Det3_123_023 * mn1OverDet;
197 m[F11] = Det3_023_023 * oneOverDet;
198 m[F12] = Det3_013_023 * mn1OverDet;
199 m[F13] = Det3_012_023 * oneOverDet;
200
201 m[F20] = Det3_123_013 * oneOverDet;
202 m[F21] = Det3_023_013 * mn1OverDet;
203 m[F22] = Det3_013_013 * oneOverDet;
204 m[F23] = Det3_012_013 * mn1OverDet;
205
206 m[F30] = Det3_123_012 * mn1OverDet;
207 m[F31] = Det3_023_012 * oneOverDet;
208 m[F32] = Det3_013_012 * mn1OverDet;
209 m[F33] = Det3_012_012 * oneOverDet;
210
211 return;
212}
F01
#define F01
Definition: MatrixInvert.cc:92
F22
#define F22
Definition: MatrixInvert.cc:103
F00
#define F00
Definition: MatrixInvert.cc:91
F02
#define F02
Definition: MatrixInvert.cc:93
F31
#define F31
Definition: MatrixInvert.cc:107
F32
#define F32
Definition: MatrixInvert.cc:108
F21
#define F21
Definition: MatrixInvert.cc:102
F03
#define F03
Definition: MatrixInvert.cc:94
F30
#define F30
Definition: MatrixInvert.cc:106
F11
#define F11
Definition: MatrixInvert.cc:97
F10
#define F10
Definition: MatrixInvert.cc:96
F33
#define F33
Definition: MatrixInvert.cc:109
F20
#define F20
Definition: MatrixInvert.cc:101
F12
#define F12
Definition: MatrixInvert.cc:98
F23
#define F23
Definition: MatrixInvert.cc:104
F13
#define F13
Definition: MatrixInvert.cc:99

Referenced by invert().

◆ invertHaywood5()

void CLHEP::HepMatrix::invertHaywood5 ( int &  ierr)
protectedvirtual

Definition at line 216 of file MatrixInvert.cc.

216 {
217
218 ifail = 0;
219
220 // Find all NECESSARY 2x2 dets: (30 of them)
221
222 double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
223 double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
224 double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
225 double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
226 double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31]; //
227 double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
228 double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31]; //
229 double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
230 double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32]; //
231 double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33]; //
232 double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
233 double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
234 double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
235 double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
236 double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
237 double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
238 double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
239 double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
240 double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
241 double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43]; //
242 double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
243 double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
244 double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
245 double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
246 double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
247 double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
248 double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
249 double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
250 double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
251 double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
252
253 // Find all NECESSARY 3x3 dets: (40 of them)
254
255 double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02
256 + m[M12]*Det2_23_01;
257 double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03
258 + m[M13]*Det2_23_01;
259 double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04
260 + m[M14]*Det2_23_01; //
261 double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03
262 + m[M13]*Det2_23_02;
263 double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04
264 + m[M14]*Det2_23_02; //
265 double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04
266 + m[M14]*Det2_23_03; //
267 double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13
268 + m[M13]*Det2_23_12;
269 double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14
270 + m[M14]*Det2_23_12; //
271 double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14
272 + m[M14]*Det2_23_13; //
273 double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24
274 + m[M14]*Det2_23_23; //
275 double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02
276 + m[M12]*Det2_24_01;
277 double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03
278 + m[M13]*Det2_24_01;
279 double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04
280 + m[M14]*Det2_24_01;
281 double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03
282 + m[M13]*Det2_24_02;
283 double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04
284 + m[M14]*Det2_24_02;
285 double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04
286 + m[M14]*Det2_24_03; //
287 double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13
288 + m[M13]*Det2_24_12;
289 double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14
290 + m[M14]*Det2_24_12;
291 double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14
292 + m[M14]*Det2_24_13; //
293 double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24
294 + m[M14]*Det2_24_23; //
295 double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02
296 + m[M12]*Det2_34_01;
297 double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03
298 + m[M13]*Det2_34_01;
299 double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04
300 + m[M14]*Det2_34_01;
301 double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03
302 + m[M13]*Det2_34_02;
303 double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04
304 + m[M14]*Det2_34_02;
305 double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04
306 + m[M14]*Det2_34_03;
307 double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13
308 + m[M13]*Det2_34_12;
309 double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14
310 + m[M14]*Det2_34_12;
311 double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14
312 + m[M14]*Det2_34_13;
313 double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24
314 + m[M14]*Det2_34_23; //
315 double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02
316 + m[M22]*Det2_34_01;
317 double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03
318 + m[M23]*Det2_34_01;
319 double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04
320 + m[M24]*Det2_34_01;
321 double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03
322 + m[M23]*Det2_34_02;
323 double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04
324 + m[M24]*Det2_34_02;
325 double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04
326 + m[M24]*Det2_34_03;
327 double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13
328 + m[M23]*Det2_34_12;
329 double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14
330 + m[M24]*Det2_34_12;
331 double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14
332 + m[M24]*Det2_34_13;
333 double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24
334 + m[M24]*Det2_34_23;
335
336 // Find all NECESSARY 4x4 dets: (25 of them)
337
338 double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023
339 + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
340 double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024
341 + m[M02]*Det3_123_014 - m[M04]*Det3_123_012; //
342 double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034
343 + m[M03]*Det3_123_014 - m[M04]*Det3_123_013; //
344 double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034
345 + m[M03]*Det3_123_024 - m[M04]*Det3_123_023; //
346 double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134
347 + m[M03]*Det3_123_124 - m[M04]*Det3_123_123; //
348 double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023
349 + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
350 double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024
351 + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
352 double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034
353 + m[M03]*Det3_124_014 - m[M04]*Det3_124_013; //
354 double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034
355 + m[M03]*Det3_124_024 - m[M04]*Det3_124_023; //
356 double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134
357 + m[M03]*Det3_124_124 - m[M04]*Det3_124_123; //
358 double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023
359 + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
360 double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024
361 + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
362 double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034
363 + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
364 double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034
365 + m[M03]*Det3_134_024 - m[M04]*Det3_134_023; //
366 double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134
367 + m[M03]*Det3_134_124 - m[M04]*Det3_134_123; //
368 double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023
369 + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
370 double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024
371 + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
372 double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034
373 + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
374 double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034
375 + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
376 double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134
377 + m[M03]*Det3_234_124 - m[M04]*Det3_234_123; //
378 double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023
379 + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
380 double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024
381 + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
382 double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034
383 + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
384 double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034
385 + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
386 double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134
387 + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
388
389 // Find the 5x5 det:
390
391 double det = m[M00]*Det4_1234_1234
392 - m[M01]*Det4_1234_0234
393 + m[M02]*Det4_1234_0134
394 - m[M03]*Det4_1234_0124
395 + m[M04]*Det4_1234_0123;
396
397 if ( det == 0 ) {
398#ifdef SINGULAR_DIAGNOSTICS
399 std::cerr << "Kramer's rule inversion of a singular 5x5 matrix: "
400 << *this << "\n";
401#endif
402 ifail = 1;
403 return;
404 }
405
406 double oneOverDet = 1.0/det;
407 double mn1OverDet = - oneOverDet;
408
409 m[M00] = Det4_1234_1234 * oneOverDet;
410 m[M01] = Det4_0234_1234 * mn1OverDet;
411 m[M02] = Det4_0134_1234 * oneOverDet;
412 m[M03] = Det4_0124_1234 * mn1OverDet;
413 m[M04] = Det4_0123_1234 * oneOverDet;
414
415 m[M10] = Det4_1234_0234 * mn1OverDet;
416 m[M11] = Det4_0234_0234 * oneOverDet;
417 m[M12] = Det4_0134_0234 * mn1OverDet;
418 m[M13] = Det4_0124_0234 * oneOverDet;
419 m[M14] = Det4_0123_0234 * mn1OverDet;
420
421 m[M20] = Det4_1234_0134 * oneOverDet;
422 m[M21] = Det4_0234_0134 * mn1OverDet;
423 m[M22] = Det4_0134_0134 * oneOverDet;
424 m[M23] = Det4_0124_0134 * mn1OverDet;
425 m[M24] = Det4_0123_0134 * oneOverDet;
426
427 m[M30] = Det4_1234_0124 * mn1OverDet;
428 m[M31] = Det4_0234_0124 * oneOverDet;
429 m[M32] = Det4_0134_0124 * mn1OverDet;
430 m[M33] = Det4_0124_0124 * oneOverDet;
431 m[M34] = Det4_0123_0124 * mn1OverDet;
432
433 m[M40] = Det4_1234_0123 * oneOverDet;
434 m[M41] = Det4_0234_0123 * mn1OverDet;
435 m[M42] = Det4_0134_0123 * oneOverDet;
436 m[M43] = Det4_0124_0123 * mn1OverDet;
437 m[M44] = Det4_0123_0123 * oneOverDet;
438
439 return;
440}
M00
#define M00
Definition: MatrixInvert.cc:61
M34
#define M34
Definition: MatrixInvert.cc:83
M03
#define M03
Definition: MatrixInvert.cc:64
M10
#define M10
Definition: MatrixInvert.cc:67
M41
#define M41
Definition: MatrixInvert.cc:86
M20
#define M20
Definition: MatrixInvert.cc:73
M42
#define M42
Definition: MatrixInvert.cc:87
M13
#define M13
Definition: MatrixInvert.cc:70
M33
#define M33
Definition: MatrixInvert.cc:82
M04
#define M04
Definition: MatrixInvert.cc:65
M23
#define M23
Definition: MatrixInvert.cc:76
M11
#define M11
Definition: MatrixInvert.cc:68
M44
#define M44
Definition: MatrixInvert.cc:89
M31
#define M31
Definition: MatrixInvert.cc:80
M02
#define M02
Definition: MatrixInvert.cc:63
M21
#define M21
Definition: MatrixInvert.cc:74
M01
#define M01
Definition: MatrixInvert.cc:62
M24
#define M24
Definition: MatrixInvert.cc:77
M12
#define M12
Definition: MatrixInvert.cc:69
M22
#define M22
Definition: MatrixInvert.cc:75
M43
#define M43
Definition: MatrixInvert.cc:88
M14
#define M14
Definition: MatrixInvert.cc:71
M30
#define M30
Definition: MatrixInvert.cc:79
M40
#define M40
Definition: MatrixInvert.cc:85
M32
#define M32
Definition: MatrixInvert.cc:81

Referenced by invert().

◆ invertHaywood6()

void CLHEP::HepMatrix::invertHaywood6 ( int &  ierr)
protectedvirtual

Definition at line 442 of file MatrixInvert.cc.

442 {
443
444 ifail = 0;
445
446 // Find all NECESSARY 2x2 dets: (45 of them)
447
448 double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
449 double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
450 double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
451 double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
452 double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40]; //
453 double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
454 double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
455 double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
456 double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41]; //
457 double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
458 double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
459 double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42]; //
460 double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
461 double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43]; //
462 double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44]; //
463 double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
464 double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
465 double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
466 double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
467 double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
468 double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
469 double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
470 double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
471 double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
472 double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
473 double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
474 double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
475 double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
476 double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
477 double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54]; //
478 double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
479 double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
480 double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
481 double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
482 double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
483 double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
484 double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
485 double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
486 double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
487 double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
488 double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
489 double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
490 double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
491 double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
492 double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
493
494 // Find all NECESSARY 3x3 dets: (80 of them)
495
496 double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
497 + m[A22]*Det2_34_01;
498 double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
499 + m[A23]*Det2_34_01;
500 double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
501 + m[A24]*Det2_34_01;
502 double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
503 + m[A25]*Det2_34_01; //
504 double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
505 + m[A23]*Det2_34_02;
506 double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
507 + m[A24]*Det2_34_02;
508 double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05
509 + m[A25]*Det2_34_02; //
510 double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
511 + m[A24]*Det2_34_03;
512 double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05
513 + m[A25]*Det2_34_03; //
514 double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05
515 + m[A25]*Det2_34_04; //
516 double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
517 + m[A23]*Det2_34_12;
518 double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
519 + m[A24]*Det2_34_12;
520 double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15
521 + m[A25]*Det2_34_12; //
522 double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
523 + m[A24]*Det2_34_13;
524 double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15
525 + m[A25]*Det2_34_13; //
526 double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15
527 + m[A25]*Det2_34_14; //
528 double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
529 + m[A24]*Det2_34_23;
530 double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
531 + m[A25]*Det2_34_23; //
532 double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
533 + m[A25]*Det2_34_24; //
534 double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
535 + m[A25]*Det2_34_34; //
536 double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
537 + m[A22]*Det2_35_01;
538 double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
539 + m[A23]*Det2_35_01;
540 double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
541 + m[A24]*Det2_35_01;
542 double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
543 + m[A25]*Det2_35_01;
544 double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
545 + m[A23]*Det2_35_02;
546 double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
547 + m[A24]*Det2_35_02;
548 double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
549 + m[A25]*Det2_35_02;
550 double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
551 + m[A24]*Det2_35_03;
552 double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
553 + m[A25]*Det2_35_03;
554 double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05
555 + m[A25]*Det2_35_04; //
556 double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
557 + m[A23]*Det2_35_12;
558 double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
559 + m[A24]*Det2_35_12;
560 double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
561 + m[A25]*Det2_35_12;
562 double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
563 + m[A24]*Det2_35_13;
564 double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
565 + m[A25]*Det2_35_13;
566 double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15
567 + m[A25]*Det2_35_14; //
568 double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
569 + m[A24]*Det2_35_23;
570 double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
571 + m[A25]*Det2_35_23;
572 double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25
573 + m[A25]*Det2_35_24; //
574 double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35
575 + m[A25]*Det2_35_34; //
576 double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
577 + m[A22]*Det2_45_01;
578 double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
579 + m[A23]*Det2_45_01;
580 double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
581 + m[A24]*Det2_45_01;
582 double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
583 + m[A25]*Det2_45_01;
584 double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
585 + m[A23]*Det2_45_02;
586 double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
587 + m[A24]*Det2_45_02;
588 double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
589 + m[A25]*Det2_45_02;
590 double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
591 + m[A24]*Det2_45_03;
592 double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
593 + m[A25]*Det2_45_03;
594 double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
595 + m[A25]*Det2_45_04;
596 double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
597 + m[A23]*Det2_45_12;
598 double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
599 + m[A24]*Det2_45_12;
600 double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
601 + m[A25]*Det2_45_12;
602 double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
603 + m[A24]*Det2_45_13;
604 double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
605 + m[A25]*Det2_45_13;
606 double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
607 + m[A25]*Det2_45_14;
608 double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
609 + m[A24]*Det2_45_23;
610 double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
611 + m[A25]*Det2_45_23;
612 double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
613 + m[A25]*Det2_45_24;
614 double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35
615 + m[A25]*Det2_45_34; //
616 double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
617 + m[A32]*Det2_45_01;
618 double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
619 + m[A33]*Det2_45_01;
620 double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
621 + m[A34]*Det2_45_01;
622 double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
623 + m[A35]*Det2_45_01;
624 double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
625 + m[A33]*Det2_45_02;
626 double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
627 + m[A34]*Det2_45_02;
628 double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
629 + m[A35]*Det2_45_02;
630 double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
631 + m[A34]*Det2_45_03;
632 double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
633 + m[A35]*Det2_45_03;
634 double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
635 + m[A35]*Det2_45_04;
636 double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
637 + m[A33]*Det2_45_12;
638 double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
639 + m[A34]*Det2_45_12;
640 double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
641 + m[A35]*Det2_45_12;
642 double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
643 + m[A34]*Det2_45_13;
644 double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
645 + m[A35]*Det2_45_13;
646 double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
647 + m[A35]*Det2_45_14;
648 double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
649 + m[A34]*Det2_45_23;
650 double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
651 + m[A35]*Det2_45_23;
652 double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
653 + m[A35]*Det2_45_24;
654 double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
655 + m[A35]*Det2_45_34;
656
657 // Find all NECESSARY 4x4 dets: (75 of them)
658
659 double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
660 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
661 double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
662 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
663 double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025
664 + m[A12]*Det3_234_015 - m[A15]*Det3_234_012; //
665 double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
666 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
667 double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
668 + m[A13]*Det3_234_015 - m[A15]*Det3_234_013; //
669 double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
670 + m[A14]*Det3_234_015 - m[A15]*Det3_234_014; //
671 double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
672 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
673 double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035
674 + m[A13]*Det3_234_025 - m[A15]*Det3_234_023; //
675 double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045
676 + m[A14]*Det3_234_025 - m[A15]*Det3_234_024; //
677 double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045
678 + m[A14]*Det3_234_035 - m[A15]*Det3_234_034; //
679 double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
680 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
681 double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135
682 + m[A13]*Det3_234_125 - m[A15]*Det3_234_123; //
683 double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145
684 + m[A14]*Det3_234_125 - m[A15]*Det3_234_124; //
685 double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145
686 + m[A14]*Det3_234_135 - m[A15]*Det3_234_134; //
687 double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245
688 + m[A14]*Det3_234_235 - m[A15]*Det3_234_234; //
689 double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
690 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
691 double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
692 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
693 double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
694 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
695 double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
696 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
697 double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
698 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
699 double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045
700 + m[A14]*Det3_235_015 - m[A15]*Det3_235_014; //
701 double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
702 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
703 double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
704 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
705 double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045
706 + m[A14]*Det3_235_025 - m[A15]*Det3_235_024; //
707 double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045
708 + m[A14]*Det3_235_035 - m[A15]*Det3_235_034; //
709 double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
710 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
711 double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
712 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
713 double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145
714 + m[A14]*Det3_235_125 - m[A15]*Det3_235_124; //
715 double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145
716 + m[A14]*Det3_235_135 - m[A15]*Det3_235_134; //
717 double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245
718 + m[A14]*Det3_235_235 - m[A15]*Det3_235_234; //
719 double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
720 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
721 double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
722 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
723 double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
724 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
725 double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
726 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
727 double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
728 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
729 double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
730 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
731 double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
732 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
733 double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
734 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
735 double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
736 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
737 double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045
738 + m[A14]*Det3_245_035 - m[A15]*Det3_245_034; //
739 double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
740 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
741 double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
742 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
743 double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
744 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
745 double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145
746 + m[A14]*Det3_245_135 - m[A15]*Det3_245_134; //
747 double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245
748 + m[A14]*Det3_245_235 - m[A15]*Det3_245_234; //
749 double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
750 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
751 double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
752 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
753 double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
754 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
755 double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
756 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
757 double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
758 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
759 double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
760 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
761 double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
762 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
763 double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
764 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
765 double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
766 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
767 double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
768 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
769 double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
770 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
771 double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
772 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
773 double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
774 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
775 double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
776 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
777 double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245
778 + m[A14]*Det3_345_235 - m[A15]*Det3_345_234; //
779 double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
780 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
781 double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
782 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
783 double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
784 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
785 double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
786 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
787 double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
788 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
789 double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
790 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
791 double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
792 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
793 double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
794 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
795 double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
796 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
797 double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
798 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
799 double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
800 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
801 double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
802 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
803 double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
804 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
805 double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
806 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
807 double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
808 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
809
810 // Find all NECESSARY 5x5 dets: (36 of them)
811
812 double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
813 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
814 double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
815 + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123; //
816 double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
817 + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124; //
818 double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
819 + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134; //
820 double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
821 + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234; //
822 double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
823 + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234; //
824 double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
825 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
826 double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
827 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
828 double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245
829 + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124; //
830 double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345
831 + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134; //
832 double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345
833 + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234; //
834 double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345
835 + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234; //
836 double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
837 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
838 double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
839 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
840 double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
841 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
842 double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345
843 + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134; //
844 double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345
845 + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234; //
846 double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345
847 + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234; //
848 double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
849 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
850 double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
851 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
852 double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
853 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
854 double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
855 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
856 double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345
857 + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234; //
858 double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345
859 + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
860 double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
861 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
862 double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
863 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
864 double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
865 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
866 double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
867 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
868 double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
869 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
870 double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345
871 + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234; //
872 double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
873 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
874 double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
875 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
876 double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
877 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
878 double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
879 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
880 double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
881 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
882 double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
883 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
884
885 // Find the determinant
886
887 double det = m[A00]*Det5_12345_12345
888 - m[A01]*Det5_12345_02345
889 + m[A02]*Det5_12345_01345
890 - m[A03]*Det5_12345_01245
891 + m[A04]*Det5_12345_01235
892 - m[A05]*Det5_12345_01234;
893
894 if ( det == 0 ) {
895#ifdef SINGULAR_DIAGNOSTICS
896 std::cerr << "Kramer's rule inversion of a singular 6x6 matrix: "
897 << *this << "\n";
898#endif
899 ifail = 1;
900 return;
901 }
902
903 double oneOverDet = 1.0/det;
904 double mn1OverDet = - oneOverDet;
905
906 m[A00] = Det5_12345_12345*oneOverDet;
907 m[A01] = Det5_02345_12345*mn1OverDet;
908 m[A02] = Det5_01345_12345*oneOverDet;
909 m[A03] = Det5_01245_12345*mn1OverDet;
910 m[A04] = Det5_01235_12345*oneOverDet;
911 m[A05] = Det5_01234_12345*mn1OverDet;
912
913 m[A10] = Det5_12345_02345*mn1OverDet;
914 m[A11] = Det5_02345_02345*oneOverDet;
915 m[A12] = Det5_01345_02345*mn1OverDet;
916 m[A13] = Det5_01245_02345*oneOverDet;
917 m[A14] = Det5_01235_02345*mn1OverDet;
918 m[A15] = Det5_01234_02345*oneOverDet;
919
920 m[A20] = Det5_12345_01345*oneOverDet;
921 m[A21] = Det5_02345_01345*mn1OverDet;
922 m[A22] = Det5_01345_01345*oneOverDet;
923 m[A23] = Det5_01245_01345*mn1OverDet;
924 m[A24] = Det5_01235_01345*oneOverDet;
925 m[A25] = Det5_01234_01345*mn1OverDet;
926
927 m[A30] = Det5_12345_01245*mn1OverDet;
928 m[A31] = Det5_02345_01245*oneOverDet;
929 m[A32] = Det5_01345_01245*mn1OverDet;
930 m[A33] = Det5_01245_01245*oneOverDet;
931 m[A34] = Det5_01235_01245*mn1OverDet;
932 m[A35] = Det5_01234_01245*oneOverDet;
933
934 m[A40] = Det5_12345_01235*oneOverDet;
935 m[A41] = Det5_02345_01235*mn1OverDet;
936 m[A42] = Det5_01345_01235*oneOverDet;
937 m[A43] = Det5_01245_01235*mn1OverDet;
938 m[A44] = Det5_01235_01235*oneOverDet;
939 m[A45] = Det5_01234_01235*mn1OverDet;
940
941 m[A50] = Det5_12345_01234*mn1OverDet;
942 m[A51] = Det5_02345_01234*oneOverDet;
943 m[A52] = Det5_01345_01234*mn1OverDet;
944 m[A53] = Det5_01245_01234*oneOverDet;
945 m[A54] = Det5_01235_01234*mn1OverDet;
946 m[A55] = Det5_01234_01234*oneOverDet;
947
948 return;
949}
A41
#define A41
Definition: MatrixInvert.cc:48
A20
#define A20
Definition: MatrixInvert.cc:33
A01
#define A01
Definition: MatrixInvert.cc:20
A53
#define A53
Definition: MatrixInvert.cc:57
A23
#define A23
Definition: MatrixInvert.cc:36
A45
#define A45
Definition: MatrixInvert.cc:52
A13
#define A13
Definition: MatrixInvert.cc:29
A00
#define A00
Definition: MatrixInvert.cc:19
A54
#define A54
Definition: MatrixInvert.cc:58
A40
#define A40
Definition: MatrixInvert.cc:47
A25
#define A25
Definition: MatrixInvert.cc:38
A02
#define A02
Definition: MatrixInvert.cc:21
A24
#define A24
Definition: MatrixInvert.cc:37
A22
#define A22
Definition: MatrixInvert.cc:35
A04
#define A04
Definition: MatrixInvert.cc:23
A30
#define A30
Definition: MatrixInvert.cc:40
A12
#define A12
Definition: MatrixInvert.cc:28
A55
#define A55
Definition: MatrixInvert.cc:59
A35
#define A35
Definition: MatrixInvert.cc:45
A50
#define A50
Definition: MatrixInvert.cc:54
A03
#define A03
Definition: MatrixInvert.cc:22
A14
#define A14
Definition: MatrixInvert.cc:30
A51
#define A51
Definition: MatrixInvert.cc:55
A21
#define A21
Definition: MatrixInvert.cc:34
A11
#define A11
Definition: MatrixInvert.cc:27
A10
#define A10
Definition: MatrixInvert.cc:26
A44
#define A44
Definition: MatrixInvert.cc:51
A05
#define A05
Definition: MatrixInvert.cc:24
A32
#define A32
Definition: MatrixInvert.cc:42
A31
#define A31
Definition: MatrixInvert.cc:41
A33
#define A33
Definition: MatrixInvert.cc:43
A42
#define A42
Definition: MatrixInvert.cc:49
A52
#define A52
Definition: MatrixInvert.cc:56
A15
#define A15
Definition: MatrixInvert.cc:31
A34
#define A34
Definition: MatrixInvert.cc:44
A43
#define A43
Definition: MatrixInvert.cc:50

Referenced by invert().

◆ num_col()

int CLHEP::HepMatrix::num_col ( ) const
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 120 of file Matrix.cc.

120{ return ncol;}

Referenced by CLHEP::HepVector::HepVector(), main(), operator()(), operator+=(), CLHEP::HepVector::operator+=(), operator-=(), CLHEP::HepVector::operator-=(), CLHEP::operator<<(), CLHEP::HepVector::operator=(), CLHEP::qr_decomp(), CLHEP::qr_inverse(), sub(), and CLHEP::tridiagonal().

◆ num_row()

int CLHEP::HepMatrix::num_row ( ) const
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 118 of file Matrix.cc.

118{ return nrow;}

Referenced by CLHEP::HepDiagMatrix::assign(), CLHEP::col_house(), HepMatrix(), main(), operator()(), operator+=(), CLHEP::HepVector::operator+=(), operator-=(), CLHEP::HepVector::operator-=(), CLHEP::operator<<(), CLHEP::qr_decomp(), CLHEP::qr_inverse(), CLHEP::row_house(), and sub().

◆ num_size()

int CLHEP::HepMatrix::num_size ( ) const
protectedvirtual

Implements CLHEP::HepGenMatrix.

Definition at line 122 of file Matrix.cc.

122{ return size_;}

◆ operator()() [1/2]

double & CLHEP::HepMatrix::operator() ( int  row,
int  col 
)
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 126 of file Matrix.cc.

127{
128#ifdef MATRIX_BOUND_CHECK
129 if(row<1 || row>num_row() || col<1 || col>num_col())
130 error("Range error in HepMatrix::operator()");
131#endif
132 return *(m.begin()+(row-1)*ncol+col-1);
133}

◆ operator()() [2/2]

const double & CLHEP::HepMatrix::operator() ( int  row,
int  col 
) const
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 135 of file Matrix.cc.

136{
137#ifdef MATRIX_BOUND_CHECK
138 if(row<1 || row>num_row() || col<1 || col>num_col())
139 error("Range error in HepMatrix::operator()");
140#endif
141 return *(m.begin()+(row-1)*ncol+col-1);
142}

◆ operator*=()

HepMatrix & CLHEP::HepMatrix::operator*= ( double  t)

Definition at line 409 of file Matrix.cc.

410{
411 SIMPLE_UOP(*=)
412 return (*this);
413}
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: DiagMatrix.cc:26

◆ operator+=() [1/4]

HepMatrix & CLHEP::HepMatrix::operator+= ( const HepDiagMatrix &  hm2)

Definition at line 451 of file DiagMatrix.cc.

452{
453 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
454 int n = num_row();
455 mIter mrr = m.begin();
456 HepMatrix::mcIter mr = hm2.m.begin();
457 for(int r=1;r<=n;r++) {
458 *mrr += *(mr++);
459 if(r<n) mrr += (n+1);
460 }
461 return (*this);
462}
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: DiagMatrix.cc:44

◆ operator+=() [2/4]

HepMatrix & CLHEP::HepMatrix::operator+= ( const HepMatrix &  hm2)

Definition at line 389 of file Matrix.cc.

390{
391 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
392 SIMPLE_BOP(+=)
393 return (*this);
394}
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: DiagMatrix.cc:31

◆ operator+=() [3/4]

HepMatrix & CLHEP::HepMatrix::operator+= ( const HepSymMatrix &  hm2)

Definition at line 559 of file SymMatrix.cc.

560{
561 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
562 HepMatrix::mcIter sjk = hm2.m.begin();
563 // j >= k
564 for(int j=0; j!=nrow; ++j) {
565 for(int k=0; k<=j; ++k) {
566 m[j*ncol+k] += *sjk;
567 // make sure this is not a diagonal element
568 if(k!=j) m[k*nrow+j] += *sjk;
569 ++sjk;
570 }
571 }
572 return (*this);
573}

◆ operator+=() [4/4]

HepMatrix & CLHEP::HepMatrix::operator+= ( const HepVector &  hm2)

Definition at line 401 of file Vector.cc.

402{
403 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),1,+=);
404 SIMPLE_BOP(+=)
405 return (*this);
406}

◆ operator-()

HepMatrix CLHEP::HepMatrix::operator- ( ) const

Definition at line 259 of file Matrix.cc.

262{
263#else
264{
265 HepMatrix hm2(nrow, ncol);
266#endif
267 mcIter a=m.begin();
268 mIter b=hm2.m.begin();
269 mcIter e=m.end();
270 for(;a<e; a++, b++) (*b) = -(*a);
271 return hm2;
272}

◆ operator-=() [1/4]

HepMatrix & CLHEP::HepMatrix::operator-= ( const HepDiagMatrix &  hm2)

Definition at line 483 of file DiagMatrix.cc.

484{
485 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),-=);
486 int n = num_row();
487 mIter mrr = m.begin();
488 HepMatrix::mcIter mr = hm2.m.begin();
489 for(int r=1;r<=n;r++) {
490 *mrr -= *(mr++);
491 if(r<n) mrr += (n+1);
492 }
493 return (*this);
494}

◆ operator-=() [2/4]

HepMatrix & CLHEP::HepMatrix::operator-= ( const HepMatrix &  hm2)

Definition at line 396 of file Matrix.cc.

397{
398 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),-=);
399 SIMPLE_BOP(-=)
400 return (*this);
401}

◆ operator-=() [3/4]

HepMatrix & CLHEP::HepMatrix::operator-= ( const HepSymMatrix &  hm2)

Definition at line 582 of file SymMatrix.cc.

583{
584 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),-=);
585 HepMatrix::mcIter sjk = hm2.m.begin();
586 // j >= k
587 for(int j=0; j!=nrow; ++j) {
588 for(int k=0; k<=j; ++k) {
589 m[j*ncol+k] -= *sjk;
590 // make sure this is not a diagonal element
591 if(k!=j) m[k*nrow+j] -= *sjk;
592 ++sjk;
593 }
594 }
595 return (*this);
596}

◆ operator-=() [4/4]

HepMatrix & CLHEP::HepMatrix::operator-= ( const HepVector &  hm2)

Definition at line 422 of file Vector.cc.

423{
424 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),1,-=);
425 SIMPLE_BOP(-=)
426 return (*this);
427}

◆ operator/=()

HepMatrix & CLHEP::HepMatrix::operator/= ( double  t)

Definition at line 403 of file Matrix.cc.

404{
405 SIMPLE_UOP(/=)
406 return (*this);
407}

◆ operator=() [1/5]

HepMatrix & CLHEP::HepMatrix::operator= ( const HepDiagMatrix &  hm1)

Definition at line 527 of file DiagMatrix.cc.

528{
529 if(hm1.nrow*hm1.nrow != size_)
530 {
531 size_ = hm1.nrow * hm1.nrow;
532 m.resize(size_);
533 }
534 nrow = hm1.nrow;
535 ncol = hm1.nrow;
536 int n = nrow;
537 m.assign(size_,0);
538 mIter mrr = m.begin();
539 HepMatrix::mcIter mr = hm1.m.begin();
540 for(int r=1;r<=n;r++) {
541 *mrr = *(mr++);
542 if(r<n) mrr += (n+1);
543 }
544 return (*this);
545}

◆ operator=() [2/5]

HepMatrix & CLHEP::HepMatrix::operator= ( const HepMatrix &  hm1)

Definition at line 415 of file Matrix.cc.

416{
417 if(hm1.nrow*hm1.ncol != size_) //??fixme?? hm1.size != size
418 {
419 size_ = hm1.nrow * hm1.ncol;
420 m.resize(size_); //??fixme?? if (size < hm1.size) m.resize(hm1.size);
421 }
422 nrow = hm1.nrow;
423 ncol = hm1.ncol;
424 m = hm1.m;
425 return (*this);
426}

◆ operator=() [3/5]

HepMatrix & CLHEP::HepMatrix::operator= ( const HepRotation &  hm1)

Definition at line 15 of file MatrixEqRotation.cc.

15 {
16 if(9!=size_) {
17 //delete &m;
18 size_ = 9;
19 m.resize(size_);
20 }
21 nrow = ncol = 3;
22 mIter hmm1;
23 hmm1 = m.begin();
24 *hmm1++ = hm1.xx();
25 *hmm1++ = hm1.xy();
26 *hmm1++ = hm1.xz();
27 *hmm1++ = hm1.yx();
28 *hmm1++ = hm1.yy();
29 *hmm1++ = hm1.yz();
30 *hmm1++ = hm1.zx();
31 *hmm1++ = hm1.zy();
32 *hmm1 = hm1.zz();
33 return (*this);
34}

◆ operator=() [4/5]

HepMatrix & CLHEP::HepMatrix::operator= ( const HepSymMatrix &  hm1)

Definition at line 617 of file SymMatrix.cc.

618{
619 // define size, rows, and columns of *this
620 nrow = ncol = hm1.nrow;
621 if(nrow*ncol != size_)
622 {
623 size_ = nrow*ncol;
624 m.resize(size_);
625 }
626 // begin copy
627 mcIter sjk = hm1.m.begin();
628 // j >= k
629 for(int j=0; j!=nrow; ++j) {
630 for(int k=0; k<=j; ++k) {
631 m[j*ncol+k] = *sjk;
632 // we could copy the diagonal element twice or check
633 // doing the check may be a tiny bit faster,
634 // so we choose that option for now
635 if(k!=j) m[k*nrow+j] = *sjk;
636 ++sjk;
637 }
638 }
639 return (*this);
640}

◆ operator=() [5/5]

HepMatrix & CLHEP::HepMatrix::operator= ( const HepVector &  hm1)

Definition at line 455 of file Vector.cc.

456{
457 if(hm1.nrow != size_)
458 {
459 size_ = hm1.nrow;
460 m.resize(size_);
461 }
462 nrow = hm1.nrow;
463 ncol = 1;
464 m = hm1.m;
465 return (*this);
466}

◆ operator[]() [1/2]

HepMatrix_row CLHEP::HepMatrix::operator[] ( int  )
inline

◆ operator[]() [2/2]

const HepMatrix_row_const CLHEP::HepMatrix::operator[] ( int  ) const
inline

◆ sub() [1/2]

HepMatrix CLHEP::HepMatrix::sub ( int  min_row,
int  max_row,
int  min_col,
int  max_col 
) const

Definition at line 193 of file Matrix.cc.

197{
198#else
199{
200 HepMatrix mret(max_row-min_row+1,max_col-min_col+1);
201#endif
202 if(max_row > num_row() || max_col >num_col())
203 error("HepMatrix::sub: Index out of range");
204 mIter a = mret.m.begin();
205 int nc = num_col();
206 mcIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
207 int rowsize = mret.num_row();
208 for(int irow=1; irow<=rowsize; ++irow) {
209 mcIter brc = b1;
210 for(int icol=0; icol<mret.num_col(); ++icol) {
211 *(a++) = *(brc++);
212 }
213 if(irow<rowsize) b1 += nc;
214 }
215 return mret;
216}

Referenced by main(), and matrix_test().

◆ sub() [2/2]

void CLHEP::HepMatrix::sub ( int  row,
int  col,
const HepMatrix &  hm1 
)

Definition at line 218 of file Matrix.cc.

219{
220 if(row <1 || row+hm1.num_row()-1 > num_row() ||
221 col <1 || col+hm1.num_col()-1 > num_col() )
222 error("HepMatrix::sub: Index out of range");
223 mcIter a = hm1.m.begin();
224 int nc = num_col();
225 mIter b1 = m.begin() + (row - 1) * nc + col - 1;
226 int rowsize = hm1.num_row();
227 for(int irow=1; irow<=rowsize; ++irow) {
228 mIter brc = b1;
229 for(int icol=0; icol<hm1.num_col(); ++icol) {
230 *(brc++) = *(a++);
231 }
232 if(irow<rowsize) b1 += nc;
233 }
234}

◆ T()

HepMatrix CLHEP::HepMatrix::T ( ) const

Definition at line 454 of file Matrix.cc.

457{
458#else
459{
460 HepMatrix mret(ncol,nrow);
461#endif
462 mcIter pme = m.begin();
463 mIter pt = mret.m.begin();
464 for( int nr=0; nr<nrow; ++nr) {
465 for( int nc=0; nc<ncol; ++nc) {
466 pt = mret.m.begin() + nr + nrow*nc;
467 (*pt) = (*pme);
468 ++pme;
469 }
470 }
471 return mret;
472}

Referenced by main(), and matrix_test().

◆ trace()

double CLHEP::HepMatrix::trace ( ) const

Definition at line 830 of file Matrix.cc.

830 {
831 double t = 0.0;
832 for (mcIter d = m.begin(); d < m.end(); d += (ncol+1) )
833 t += *d;
834 return t;
835}

Friends And Related Function Documentation

◆ back_solve [1/2]

void back_solve ( const HepMatrix &  R,
HepMatrix *  b 
)
friend

Definition at line 90 of file MatrixLinear.cc.

91{
92 int n = R.num_col();
93 int nb = b->num_row();
94 int nc = b->num_col();
95 HepMatrix::mIter bbi = b->m.begin() + (nb - 2) * nc;
96 for (int i=1;i<=b->num_col();i++) {
97 (*b)(b->num_row(),i) /= R(b->num_row(),b->num_row());
98 HepMatrix::mcIter Rrr = R.m.begin() + (nb-2) * (n+1);
99 HepMatrix::mIter bri = bbi;
100 for (int r=b->num_row()-1;r>=1;--r) {
101 HepMatrix::mIter bci = bri + nc;
102 HepMatrix::mcIter Rrc = Rrr+1;
103 for (int c=r+1;c<=b->num_row();c++) {
104 (*bri)-= (*(Rrc++)) * (*bci);
105 if(c<b->num_row()) bci += nc;
106 }
107 (*bri)/= (*Rrr);
108 if(r>1) {
109 Rrr -= (n+1);
110 bri -= nc;
111 }
112 }
113 bbi++;
114 }
115}

◆ back_solve [2/2]

void back_solve ( const HepMatrix &  R,
HepVector *  b 
)
friend

Definition at line 63 of file MatrixLinear.cc.

64{
65 (*b)(b->num_row()) /= R(b->num_row(),b->num_row());
66 int n = R.num_col();
67 int nb = b->num_row();
68 HepMatrix::mIter br = b->m.begin() + b->num_row() - 2;
69 HepMatrix::mcIter Rrr = R.m.begin() + (nb-2) * (n+1);
70 for (int r=b->num_row()-1;r>=1;--r) {
71 HepMatrix::mIter bc = br+1;
72 HepMatrix::mcIter Rrc = Rrr+1;
73 for (int c=r+1;c<=b->num_row();c++) {
74 (*br)-=(*(Rrc++))*(*(bc++));
75 }
76 (*br)/=(*Rrr);
77 if(r>1) {
78 br--;
79 Rrr -= n+1;
80 }
81 }
82}

◆ col_givens

void col_givens ( HepMatrix *  A,
double  c,
double  s,
int  k1,
int  k2,
int  rowmin = 1,
int  rowmax = 0 
)
friend

Definition at line 124 of file MatrixLinear.cc.

125 {
126 if (row_max<=0) row_max = A->num_row();
127 int n = A->num_col();
128 HepMatrix::mIter Ajk1 = A->m.begin() + (row_min - 1) * n + k1 - 1;
129 HepMatrix::mIter Ajk2 = A->m.begin() + (row_min - 1) * n + k2 - 1;
130 for (int j=row_min;j<=row_max;j++) {
131 double tau1=(*Ajk1); double tau2=(*Ajk2);
132 (*Ajk1)=c*tau1-ds*tau2;(*Ajk2)=ds*tau1+c*tau2;
133 if(j<row_max) {
134 Ajk1 += n;
135 Ajk2 += n;
136 }
137 }
138}

◆ col_house

void col_house ( HepMatrix *  a,
const HepMatrix &  v,
double  vnormsq,
int  row,
int  col,
int  row_start,
int  col_start 
)
friend

Definition at line 154 of file MatrixLinear.cc.

155 {
156 double beta=-2/vnormsq;
157
158 // This is a fast way of calculating w=beta*A.sub(row,n,col,n)*v
159
160 HepVector w(a->num_col()-col+1,0);
161/* not tested */
162 HepMatrix::mIter wptr = w.m.begin();
163 int n = a->num_col();
164 int nv = v.num_col();
165 HepMatrix::mIter acrb = a->m.begin() + (col-1) * n + (row-1);
166 int c;
167 for (c=col;c<=a->num_col();c++) {
168 HepMatrix::mcIter vp = v.m.begin() + (row_start-1) * nv + (col_start-1);
169 HepMatrix::mcIter acr = acrb;
170 for (int r=row;r<=a->num_row();r++) {
171 (*wptr)+=(*(acr++))*(*vp);
172 vp += nv;
173 }
174 wptr++;
175 if(c<a->num_col()) acrb += n;
176 }
177 w*=beta;
178
179 // Fast way of calculating A.sub=A.sub+w*v.T()
180
181 HepMatrix::mIter arcb = a->m.begin() + (row-1) * n + col-1;
182 wptr = w.m.begin();
183 for (int r=row; r<=a->num_row();r++) {
184 HepMatrix::mIter arc = arcb;
185 HepMatrix::mcIter vp = v.m.begin() + (row_start-1) * nv + col_start;
186 for (c=col;c<=a->num_col();c++) {
187 (*(arc++))+=(*vp)*(*wptr);
188 vp += nv;
189 }
190 wptr++;
191 if(r<a->num_row()) arcb += n;
192 }
193}
CLHEP::HepMatrix::HepVector
friend class HepVector
Definition: Matrix.h:343

◆ HepDiagMatrix

friend class HepDiagMatrix
friend

Definition at line 345 of file Matrix.h.

◆ HepMatrix_row

friend class HepMatrix_row
friend

Definition at line 341 of file Matrix.h.

◆ HepMatrix_row_const

friend class HepMatrix_row_const
friend

Definition at line 342 of file Matrix.h.

◆ HepSymMatrix

friend class HepSymMatrix
friend

Definition at line 344 of file Matrix.h.

◆ HepVector

friend class HepVector
friend

Definition at line 343 of file Matrix.h.

◆ house

HepVector house ( const HepMatrix &  a,
int  row = 1,
int  col = 1 
)
friend

Definition at line 371 of file MatrixLinear.cc.

372{
373 HepVector v(a.num_row()-row+1);
374/* not tested */
375 int n = a.num_col();
376 HepMatrix::mcIter aic = a.m.begin() + (row - 1) * n + (col - 1) ;
377 HepMatrix::mIter vp = v.m.begin();
378 for (int i=row;i<=a.num_row();i++) {
379 (*(vp++))=(*aic);
380 aic += n;
381 }
382 v(1)+=sign(a(row,col))*v.norm();
383 return v;
384}

◆ house_with_update [1/2]

void house_with_update ( HepMatrix *  a,
HepMatrix *  v,
int  row = 1,
int  col = 1 
)
friend

Definition at line 424 of file MatrixLinear.cc.

425{
426 double normsq=0;
427 int nv = v->num_col();
428 int na = a->num_col();
429 HepMatrix::mIter vrc = v->m.begin() + (row-1) * nv + (col-1);
430 HepMatrix::mIter arc = a->m.begin() + (row-1) * na + (col-1);
431 int r;
432 for (r=row;r<=a->num_row();r++) {
433 (*vrc)=(*arc);
434 normsq+=(*vrc)*(*vrc);
435 if(r<a->num_row()) {
436 vrc += nv;
437 arc += na;
438 }
439 }
440 double norm=sqrt(normsq);
441 vrc = v->m.begin() + (row-1) * nv + (col-1);
442 normsq-=(*vrc)*(*vrc);
443 (*vrc)+=sign((*a)(row,col))*norm;
444 normsq+=(*vrc)*(*vrc);
445 (*a)(row,col)=-sign((*a)(row,col))*norm;
446 if (row<a->num_row()) {
447 arc = a->m.begin() + row * na + (col-1);
448 for (r=row+1;r<=a->num_row();r++) {
449 (*arc)=0;
450 if(r<a->num_row()) arc += na;
451 }
452 row_house(a,*v,normsq,row,col+1,row,col);
453 }
454}
CLHEP::HepMatrix::row_house
friend void row_house(HepMatrix *, const HepMatrix &, double, int, int, int, int)
Definition: MatrixLinear.cc:652
CLHEP::norm
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:54

◆ house_with_update [2/2]

void house_with_update ( HepMatrix *  a,
int  row = 1,
int  col = 1 
)
friend

Definition at line 396 of file MatrixLinear.cc.

397{
398 HepVector v(a->num_row()-row+1);
399/* not tested */
400 HepMatrix::mIter vp = v.m.begin();
401 int n = a->num_col();
402 HepMatrix::mIter arc = a->m.begin() + (row-1) * n + (col-1);
403 int r;
404 for (r=row;r<=a->num_row();r++) {
405 (*(vp++))=(*arc);
406 if(r<a->num_row()) arc += n;
407 }
408 double normsq=v.normsq();
409 double norm=sqrt(normsq);
410 normsq-=v(1)*v(1);
411 v(1)+=sign((*a)(row,col))*norm;
412 normsq+=v(1)*v(1);
413 (*a)(row,col)=-sign((*a)(row,col))*norm;
414 if (row<a->num_row()) {
415 arc = a->m.begin() + row * n + (col-1);
416 for (r=row+1;r<=a->num_row();r++) {
417 (*arc)=0;
418 if(r<a->num_row()) arc += n;
419 }
420 row_house(a,v,normsq,row,col+1);
421 }
422}

◆ house_with_update2

void house_with_update2 ( HepSymMatrix *  a,
HepMatrix *  v,
int  row = 1,
int  col = 1 
)
friend

Definition at line 462 of file MatrixLinear.cc.

463{
464 double normsq=0;
465 int nv = v->num_col();
466 HepMatrix::mIter vrc = v->m.begin() + (row-1) * nv + (col-1);
467 HepMatrix::mIter arc = a->m.begin() + (row-1) * row / 2 + (col-1);
468 int r;
469 for (r=row;r<=a->num_row();r++)
470 {
471 (*vrc)=(*arc);
472 normsq+=(*vrc)*(*vrc);
473 if(r<a->num_row()) {
474 arc += r;
475 vrc += nv;
476 }
477 }
478 double norm=sqrt(normsq);
479 vrc = v->m.begin() + (row-1) * nv + (col-1);
480 arc = a->m.begin() + (row-1) * row / 2 + (col-1);
481 (*vrc)+=sign(*arc)*norm;
482 (*arc)=-sign(*arc)*norm;
483 arc += row;
484 for (r=row+1;r<=a->num_row();r++) {
485 (*arc)=0;
486 if(r<a->num_row()) arc += r;
487 }
488}

◆ operator* [1/8]

HepMatrix operator* ( const HepDiagMatrix &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 392 of file DiagMatrix.cc.

395{
396#else
397{
398 HepMatrix mret(hm1.num_row(),hm2.num_col());
399#endif
400 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
401 HepMatrix::mcIter mit1=hm2.m.begin();
402 HepMatrix::mIter mir=mret.m.begin();
403 HepMatrix::mcIter mrr = hm1.m.begin();
404 for(int irow=1;irow<=hm2.num_row();irow++) {
405 for(int icol=1;icol<=hm2.num_col();icol++) {
406 *(mir++) = *(mit1++) * (*mrr);
407 }
408 mrr++;
409 }
410 return mret;
411}
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: DiagMatrix.cc:49

◆ operator* [2/8]

HepMatrix operator* ( const HepMatrix &  hm1,
const HepDiagMatrix &  hm2 
)
friend

Definition at line 372 of file DiagMatrix.cc.

375{
376#else
377 {
378 HepMatrix mret(hm1.num_row(),hm2.num_col());
379#endif
380 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
381 HepMatrix::mcIter mit1=hm1.m.begin();
382 HepMatrix::mIter mir=mret.m.begin();
383 for(int irow=1;irow<=hm1.num_row();irow++) {
384 HepMatrix::mcIter mcc = hm2.m.begin();
385 for(int icol=1;icol<=hm1.num_col();icol++) {
386 *(mir++) = *(mit1++) * (*(mcc++));
387 }
388 }
389 return mret;
390 }

◆ operator* [3/8]

HepMatrix operator* ( const HepMatrix &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 349 of file Matrix.cc.

352{
353#else
354{
355 // initialize matrix to 0.0
356 HepMatrix mret(hm1.nrow,hm2.ncol,0);
357#endif
358 CHK_DIM_1(hm1.ncol,hm2.nrow,*);
359
360 int m1cols = hm1.ncol;
361 int m2cols = hm2.ncol;
362
363 for (int i=0; i<hm1.nrow; i++)
364 {
365 for (int j=0; j<m1cols; j++)
366 {
367 double temp = hm1.m[i*m1cols+j];
368 HepMatrix::mIter pt = mret.m.begin() + i*m2cols;
369
370 // Loop over k (the column index in matrix hm2)
371 HepMatrix::mcIter pb = hm2.m.begin() + m2cols*j;
372 const HepMatrix::mcIter pblast = pb + m2cols;
373 while (pb < pblast)
374 {
375 (*pt) += temp * (*pb);
376 pb++;
377 pt++;
378 }
379 }
380 }
381
382 return mret;
383}

◆ operator* [4/8]

HepMatrix operator* ( const HepMatrix &  hm1,
const HepSymMatrix &  hm2 
)
friend

Definition at line 353 of file SymMatrix.cc.

356{
357#else
358 {
359 HepMatrix mret(hm1.num_row(),hm2.num_col());
360#endif
361 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
362 HepMatrix::mcIter mit1, mit2, sp,snp; //mit2=0
363 double temp;
364 HepMatrix::mIter mir=mret.m.begin();
365 for(mit1=hm1.m.begin();
366 mit1<hm1.m.begin()+hm1.num_row()*hm1.num_col();
367 mit1 = mit2)
368 {
369 snp=hm2.m.begin();
370 for(int step=1;step<=hm2.num_row();++step)
371 {
372 mit2=mit1;
373 sp=snp;
374 snp+=step;
375 temp=0;
376 while(sp<snp)
377 temp+=*(sp++)*(*(mit2++));
378 if( step<hm2.num_row() ) { // only if we aren't on the last row
379 sp+=step-1;
380 for(int stept=step+1;stept<=hm2.num_row();stept++)
381 {
382 temp+=*sp*(*(mit2++));
383 if(stept<hm2.num_row()) sp+=stept;
384 }
385 } // if(step
386 *(mir++)=temp;
387 } // for(step
388 } // for(mit1
389 return mret;
390 }

◆ operator* [5/8]

HepVector operator* ( const HepMatrix &  hm1,
const HepVector &  hm2 
)
friend

Definition at line 354 of file Vector.cc.

357{
358#else
359{
360 HepVector mret(hm1.num_row());
361#endif
362 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
363 HepGenMatrix::mcIter hm1p,hm2p,vp;
364 HepGenMatrix::mIter m3p;
365 double temp;
366 m3p=mret.m.begin();
367 for(hm1p=hm1.m.begin();hm1p<hm1.m.begin()+hm1.num_row()*hm1.num_col();hm1p=hm2p)
368 {
369 temp=0;
370 vp=hm2.m.begin();
371 hm2p=hm1p;
372 while(hm2p<hm1p+hm1.num_col())
373 temp+=(*(hm2p++))*(*(vp++));
374 *(m3p++)=temp;
375 }
376 return mret;
377}

◆ operator* [6/8]

HepMatrix operator* ( const HepSymMatrix &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 392 of file SymMatrix.cc.

395{
396#else
397{
398 HepMatrix mret(hm1.num_row(),hm2.num_col());
399#endif
400 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
401 int step,stept;
402 HepMatrix::mcIter mit1,mit2,sp,snp;
403 double temp;
404 HepMatrix::mIter mir=mret.m.begin();
405 for(step=1,snp=hm1.m.begin();step<=hm1.num_row();snp+=step++)
406 for(mit1=hm2.m.begin();mit1<hm2.m.begin()+hm2.num_col();mit1++)
407 {
408 mit2=mit1;
409 sp=snp;
410 temp=0;
411 while(sp<snp+step)
412 {
413 temp+=*mit2*(*(sp++));
414 if( hm2.num_size()-(mit2-hm2.m.begin())>hm2.num_col() ){
415 mit2+=hm2.num_col();
416 }
417 }
418 if(step<hm1.num_row()) { // only if we aren't on the last row
419 sp+=step-1;
420 for(stept=step+1;stept<=hm1.num_row();stept++)
421 {
422 temp+=*mit2*(*sp);
423 if(stept<hm1.num_row()) {
424 mit2+=hm2.num_col();
425 sp+=stept;
426 }
427 }
428 } // if(step
429 *(mir++)=temp;
430 } // for(mit1
431 return mret;
432}

◆ operator* [7/8]

HepMatrix operator* ( const HepSymMatrix &  hm1,
const HepSymMatrix &  hm2 
)
friend

Definition at line 434 of file SymMatrix.cc.

437{
438#else
439{
440 HepMatrix mret(hm1.num_row(),hm1.num_row());
441#endif
442 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
443 int step1,stept1,step2,stept2;
444 HepMatrix::mcIter snp1,sp1,snp2,sp2;
445 double temp;
446 HepMatrix::mIter mr = mret.m.begin();
447 snp1=hm1.m.begin();
448 for(step1=1;step1<=hm1.num_row();++step1) {
449 snp2=hm2.m.begin();
450 for(step2=1;step2<=hm2.num_row();++step2)
451 {
452 sp1=snp1;
453 sp2=snp2;
454 snp2+=step2;
455 temp=0;
456 if(step1<step2)
457 {
458 while(sp1<snp1+step1) {
459 temp+=(*(sp1++))*(*(sp2++));
460 }
461 sp1+=step1-1;
462 for(stept1=step1+1;stept1!=step2+1;++stept1) {
463 temp+=(*sp1)*(*(sp2++));
464 if(stept1<hm2.num_row()) sp1+=stept1;
465 }
466 if(step2<hm2.num_row()) { // only if we aren't on the last row
467 sp2+=step2-1;
468 for(stept2=step2+1;stept2<=hm2.num_row();stept1++,stept2++) {
469 temp+=(*sp1)*(*sp2);
470 if(stept2<hm2.num_row()) {
471 sp1+=stept1;
472 sp2+=stept2;
473 }
474 } // for(stept2
475 } // if(step2
476 } // step1<step2
477 else
478 {
479 while(sp2<snp2) {
480 temp+=(*(sp1++))*(*(sp2++));
481 }
482 if(step2<hm2.num_row()) { // only if we aren't on the last row
483 sp2+=step2-1;
484 for(stept2=step2+1;stept2!=step1+1;stept2++) {
485 temp+=(*(sp1++))*(*sp2);
486 if(stept2<hm1.num_row()) sp2+=stept2;
487 }
488 if(step1<hm1.num_row()) { // only if we aren't on the last row
489 sp1+=step1-1;
490 for(stept1=step1+1;stept1<=hm1.num_row();stept1++,stept2++) {
491 temp+=(*sp1)*(*sp2);
492 if(stept1<hm1.num_row()) {
493 sp1+=stept1;
494 sp2+=stept2;
495 }
496 } // for(stept1
497 } // if(step1
498 } // if(step2
499 } // else
500 *(mr++)=temp;
501 } // for(step2
502 if(step1<hm1.num_row()) snp1+=step1;
503 } // for(step1
504 return mret;
505}

◆ operator* [8/8]

HepMatrix operator* ( const HepVector &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 379 of file Vector.cc.

382{
383#else
384{
385 HepMatrix mret(hm1.num_row(),hm2.num_col());
386#endif
387 CHK_DIM_1(1,hm2.num_row(),*);
388 HepGenMatrix::mcIter hm1p;
389 HepMatrix::mcIter hm2p;
390 HepMatrix::mIter mrp=mret.m.begin();
391 for(hm1p=hm1.m.begin();hm1p<hm1.m.begin()+hm1.num_row();hm1p++)
392 for(hm2p=hm2.m.begin();hm2p<hm2.m.begin()+hm2.num_col();hm2p++)
393 *(mrp++)=*hm1p*(*hm2p);
394 return mret;
395}

◆ operator+

HepMatrix operator+ ( const HepMatrix &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 276 of file Matrix.cc.

279{
280#else
281{
282 HepMatrix mret(hm1.nrow, hm1.ncol);
283#endif
284 CHK_DIM_2(hm1.num_row(),hm2.num_row(), hm1.num_col(),hm2.num_col(),+);
285 SIMPLE_TOP(+)
286 return mret;
287}
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: DiagMatrix.cc:37

◆ operator-

HepMatrix operator- ( const HepMatrix &  hm1,
const HepMatrix &  hm2 
)
friend

Definition at line 293 of file Matrix.cc.

296{
297#else
298{
299 HepMatrix mret(hm1.num_row(), hm1.num_col());
300#endif
301 CHK_DIM_2(hm1.num_row(),hm2.num_row(),
302 hm1.num_col(),hm2.num_col(),-);
303 SIMPLE_TOP(-)
304 return mret;
305}

◆ qr_solve [1/2]

HepMatrix qr_solve ( HepMatrix *  A,
const HepMatrix &  b 
)
friend

Definition at line 738 of file MatrixLinear.cc.

739{
740 HepMatrix Q=qr_decomp(A);
741 // Quick way to to Q.T*b.
742 HepMatrix b2(Q.num_col(),b.num_col(),0);
743 int nb = b.num_col();
744 int nq = Q.num_col();
745 HepMatrix::mcIter b1i = b.m.begin();
746 HepMatrix::mIter b21i = b2.m.begin();
747 for (int i=1;i<=b.num_col();i++) {
748 HepMatrix::mIter Q1r = Q.m.begin();
749 HepMatrix::mIter b2ri = b21i;
750 for (int r=1;r<=b2.num_row();r++) {
751 HepMatrix::mIter Qcr = Q1r;
752 HepMatrix::mcIter bci = b1i;
753 for (int c=1;c<=b.num_row();c++) {
754 *b2ri += (*Qcr) * (*bci);
755 if(c<b.num_row()) {
756 Qcr += nq;
757 bci += nb;
758 }
759 }
760 Q1r++;
761 if(r<b2.num_row()) b2ri += nb;
762 }
763 b1i++;
764 b21i++;
765 }
766 back_solve(*A,&b2);
767 return b2;
768}
CLHEP::HepMatrix::back_solve
friend void back_solve(const HepMatrix &R, HepVector *b)
Definition: MatrixLinear.cc:63
CLHEP::qr_decomp
void qr_decomp(HepMatrix *A, HepMatrix *hsm)
Definition: MatrixLinear.cc:898

◆ qr_solve [2/2]

HepVector qr_solve ( HepMatrix *  A,
const HepVector &  b 
)
friend

Definition at line 710 of file MatrixLinear.cc.

711{
712 HepMatrix Q=qr_decomp(A);
713 // Quick way to to Q.T*b.
714 HepVector b2(Q.num_col(),0);
715 HepMatrix::mIter b2r = b2.m.begin();
716 HepMatrix::mIter Qr = Q.m.begin();
717 int n = Q.num_col();
718 for (int r=1;r<=b2.num_row();r++) {
719 HepMatrix::mcIter bc = b.m.begin();
720 HepMatrix::mIter Qcr = Qr;
721 for (int c=1;c<=b.num_row();c++) {
722 *b2r += (*Qcr) * (*(bc++));
723 if(c<b.num_row()) Qcr += n;
724 }
725 b2r++;
726 Qr++;
727 }
728 back_solve(*A,&b2);
729 return b2;
730}

◆ row_givens

void row_givens ( HepMatrix *  A,
double  c,
double  s,
int  k1,
int  k2,
int  colmin = 1,
int  colmax = 0 
)
friend

Definition at line 587 of file MatrixLinear.cc.

588 {
589 /* not tested */
590 if (col_max==0) col_max = A->num_col();
591 int n = A->num_col();
592 HepMatrix::mIter Ak1j = A->m.begin() + (k1-1) * n + (col_min-1);
593 HepMatrix::mIter Ak2j = A->m.begin() + (k2-1) * n + (col_min-1);
594 for (int j=col_min;j<=col_max;j++) {
595 double tau1=(*Ak1j); double tau2=(*Ak2j);
596 (*(Ak1j++))=c*tau1-ds*tau2;(*(Ak2j++))=ds*tau1+c*tau2;
597 }
598}

◆ row_house [1/2]

void row_house ( HepMatrix *  a,
const HepMatrix &  v,
double  vnormsq,
int  row,
int  col,
int  row_start,
int  col_start 
)
friend

Definition at line 652 of file MatrixLinear.cc.

653 {
654 double beta=-2/vnormsq;
655
656 // This is a fast way of calculating w=beta*A.sub(row,n,col,n).T()*v
657 HepVector w(a->num_col()-col+1,0);
658 int na = a->num_col();
659 int nv = v.num_col();
660 HepMatrix::mIter wptr = w.m.begin();
661 HepMatrix::mIter arcb = a->m.begin() + (row-1) * na + (col-1);
662 HepMatrix::mcIter vpcb = v.m.begin() + (row_start-1) * nv + (col_start-1);
663 int c;
664 for (c=col;c<=a->num_col();c++) {
665 HepMatrix::mIter arc = arcb;
666 HepMatrix::mcIter vpc = vpcb;
667 for (int r=row;r<=a->num_row();r++) {
668 (*wptr)+=(*arc)*(*vpc);
669 if(r<a->num_row()) {
670 arc += na;
671 vpc += nv;
672 }
673 }
674 wptr++;
675 arcb++;
676 }
677 w*=beta;
678
679 arcb = a->m.begin() + (row-1) * na + (col-1);
680 HepMatrix::mcIter vpc = v.m.begin() + (row_start-1) * nv + (col_start-1);
681 for (int r=row; r<=a->num_row();r++) {
682 HepMatrix::mIter arc = arcb;
683 HepMatrix::mIter wptr2 = w.m.begin();
684 for (c=col;c<=a->num_col();c++) {
685 (*(arc++))+=(*vpc)*(*(wptr2++));
686 }
687 if(r<a->num_row()) {
688 arcb += na;
689 vpc += nv;
690 }
691 }
692}

◆ row_house [2/2]

void row_house ( HepMatrix *  a,
const HepVector &  v,
double  vnormsq,
int  row = 1,
int  col = 1 
)
friend

Definition at line 613 of file MatrixLinear.cc.

614 {
615 double beta=-2/vnormsq;
616
617 // This is a fast way of calculating w=beta*A.sub(row,n,col,n).T()*v
618
619 HepVector w(a->num_col()-col+1,0);
620/* not tested */
621 int na = a->num_col();
622 HepMatrix::mIter wptr = w.m.begin();
623 HepMatrix::mIter arcb = a->m.begin() + (row-1) * na + (col-1);
624 int c;
625 for (c=col;c<=a->num_col();c++) {
626 HepMatrix::mcIter vp = v.m.begin();
627 HepMatrix::mIter arc = arcb;
628 for (int r=row;r<=a->num_row();r++) {
629 (*wptr)+=(*arc)*(*(vp++));
630 if(r<a->num_row()) arc += na;
631 }
632 wptr++;
633 arcb++;
634 }
635 w*=beta;
636
637 // Fast way of calculating A.sub=A.sub+v*w.T()
638
639 arcb = a->m.begin() + (row-1) * na + (col-1);
640 HepMatrix::mcIter vp = v.m.begin();
641 for (int r=row; r<=a->num_row();r++) {
642 HepMatrix::mIter wptr2 = w.m.begin();
643 HepMatrix::mIter arc = arcb;
644 for (c=col;c<=a->num_col();c++) {
645 (*(arc++))+=(*vp)*(*(wptr2++));
646 }
647 vp++;
648 if(r<a->num_row()) arcb += na;
649 }
650}

◆ solve

HepVector solve ( const HepMatrix &  a,
const HepVector &  v 
)
friend

Definition at line 575 of file Vector.cc.

578{
579#else
580{
581 HepVector vret(v);
582#endif
583 static CLHEP_THREAD_LOCAL int max_array = 20;
584 static CLHEP_THREAD_LOCAL int *ir = new int [max_array+1];
585
586 if(a.ncol != a.nrow)
587 HepGenMatrix::error("Matrix::solve Matrix is not NxN");
588 if(a.ncol != v.nrow)
589 HepGenMatrix::error("Matrix::solve Vector has wrong number of rows");
590
591 int n = a.ncol;
592 if (n > max_array) {
593 delete [] ir;
594 max_array = n;
595 ir = new int [max_array+1];
596 }
597 double det;
598 HepMatrix mt(a);
599 int i = mt.dfact_matrix(det, ir);
600 if (i!=0) {
601 for (i=1;i<=n;i++) vret(i) = 0;
602 return vret;
603 }
604 double s21, s22;
605 int nxch = ir[n];
606 if (nxch!=0) {
607 for (int hmm=1;hmm<=nxch;hmm++) {
608 int ij = ir[hmm];
609 i = ij >> 12;
610 int j = ij%4096;
611 double te = vret(i);
612 vret(i) = vret(j);
613 vret(j) = te;
614 }
615 }
616 vret(1) = mt(1,1) * vret(1);
617 if (n!=1) {
618 for (i=2;i<=n;i++) {
619 s21 = -vret(i);
620 for (int j=1;j<i;j++) {
621 s21 += mt(i,j) * vret(j);
622 }
623 vret(i) = -mt(i,i)*s21;
624 }
625 for (i=1;i<n;i++) {
626 int nmi = n-i;
627 s22 = -vret(nmi);
628 for (int j=1;j<=i;j++) {
629 s22 += mt(nmi,n-j+1) * vret(n-j+1);
630 }
631 vret(nmi) = -s22;
632 }
633 }
634 return vret;
635}

◆ swap

void swap ( HepMatrix &  hm1,
HepMatrix &  hm2 
)
friend

◆ tridiagonal

void tridiagonal ( HepSymMatrix *  a,
HepMatrix *  hsm 
)
friend

Definition at line 777 of file MatrixLinear.cc.

778{
779 int nh = hsm->num_col();
780 for (int k=1;k<=a->num_col()-2;k++) {
781
782 // If this row is already in tridiagonal form, we can skip the
783 // transformation.
784
785 double scale=0;
786 HepMatrix::mIter ajk = a->m.begin() + k * (k+5) / 2;
787 int j;
788 for (j=k+2;j<=a->num_row();j++) {
789 scale+=fabs(*ajk);
790 if(j<a->num_row()) ajk += j;
791 }
792 if (scale ==0) {
793 HepMatrix::mIter hsmjkp = hsm->m.begin() + k * (nh+1) - 1;
794 for (j=k+1;j<=hsm->num_row();j++) {
795 *hsmjkp = 0;
796 if(j<hsm->num_row()) hsmjkp += nh;
797 }
798 } else {
799 house_with_update2(a,hsm,k+1,k);
800 double normsq=0;
801 HepMatrix::mIter hsmrptrkp = hsm->m.begin() + k * (nh+1) - 1;
802 int rptr;
803 for (rptr=k+1;rptr<=hsm->num_row();rptr++) {
804 normsq+=(*hsmrptrkp)*(*hsmrptrkp);
805 if(rptr<hsm->num_row()) hsmrptrkp += nh;
806 }
807 HepVector p(a->num_row()-k,0);
808 rptr=k+1;
809 HepMatrix::mIter pr = p.m.begin();
810 int r;
811 for (r=1;r<=p.num_row();r++) {
812 HepMatrix::mIter hsmcptrkp = hsm->m.begin() + k * (nh+1) - 1;
813 int cptr;
814 for (cptr=k+1;cptr<=rptr;cptr++) {
815 (*pr)+=a->fast(rptr,cptr)*(*hsmcptrkp);
816 if(cptr<a->num_col()) hsmcptrkp += nh;
817 }
818 for (;cptr<=a->num_col();cptr++) {
819 (*pr)+=a->fast(cptr,rptr)*(*hsmcptrkp);
820 if(cptr<a->num_col()) hsmcptrkp += nh;
821 }
822 (*pr)*=2/normsq;
823 rptr++;
824 pr++;
825 }
826 double pdotv=0;
827 pr = p.m.begin();
828 hsmrptrkp = hsm->m.begin() + k * (nh+1) - 1;
829 for (r=1;r<=p.num_row();r++) {
830 pdotv+=(*(pr++))*(*hsmrptrkp);
831 if(r<p.num_row()) hsmrptrkp += nh;
832 }
833 pr = p.m.begin();
834 hsmrptrkp = hsm->m.begin() + k * (nh+1) - 1;
835 for (r=1;r<=p.num_row();r++) {
836 (*(pr++))-=pdotv*(*hsmrptrkp)/normsq;
837 if(r<p.num_row()) hsmrptrkp += nh;
838 }
839 rptr=k+1;
840 pr = p.m.begin();
841 hsmrptrkp = hsm->m.begin() + k * (nh+1) - 1;
842 for (r=1;r<=p.num_row();r++) {
843 int cptr=k+1;
844 HepMatrix::mIter mpc = p.m.begin();
845 HepMatrix::mIter hsmcptrkp = hsm->m.begin() + k * (nh+1) - 1;
846 for (int c=1;c<=r;c++) {
847 a->fast(rptr,cptr)-= (*hsmrptrkp)*(*(mpc++))+(*pr)*(*hsmcptrkp);
848 cptr++;
849 if(c<r) hsmcptrkp += nh;
850 }
851 pr++;
852 rptr++;
853 if(r<p.num_row()) hsmrptrkp += nh;
854 }
855 }
856 }
857}
CLHEP::HepMatrix::house_with_update2
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
Definition: MatrixLinear.cc:462
The documentation for this class was generated from the following files: