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::HepSymMatrix Class Reference

#include <SymMatrix.h>

+ Inheritance diagram for CLHEP::HepSymMatrix:

Classes

class  HepSymMatrix_row
 
class  HepSymMatrix_row_const
 

Public Member Functions

 HepSymMatrix ()
 
 HepSymMatrix (int p)
 
 HepSymMatrix (int p, int)
 
 HepSymMatrix (int p, HepRandom &r)
 
 HepSymMatrix (const HepSymMatrix &hm1)
 
 HepSymMatrix (const HepDiagMatrix &hm1)
 
virtual ~HepSymMatrix ()
 
int num_row () const
 
int num_col () const
 
const double & operator() (int row, int col) const
 
double & operator() (int row, int col)
 
const double & fast (int row, int col) const
 
double & fast (int row, int col)
 
void assign (const HepMatrix &hm2)
 
void assign (const HepSymMatrix &hm2)
 
HepSymMatrix & operator*= (double t)
 
HepSymMatrix & operator/= (double t)
 
HepSymMatrix & operator+= (const HepSymMatrix &hm2)
 
HepSymMatrix & operator+= (const HepDiagMatrix &hm2)
 
HepSymMatrix & operator-= (const HepSymMatrix &hm2)
 
HepSymMatrix & operator-= (const HepDiagMatrix &hm2)
 
HepSymMatrix & operator= (const HepSymMatrix &hm2)
 
HepSymMatrix & operator= (const HepDiagMatrix &hm2)
 
HepSymMatrix operator- () const
 
HepSymMatrix T () const
 
HepSymMatrix apply (double(*f)(double, int, int)) const
 
HepSymMatrix similarity (const HepMatrix &hm1) const
 
HepSymMatrix similarity (const HepSymMatrix &hm1) const
 
HepSymMatrix similarityT (const HepMatrix &hm1) const
 
double similarity (const HepVector &v) const
 
HepSymMatrix sub (int min_row, int max_row) const
 
void sub (int row, const HepSymMatrix &hm1)
 
HepSymMatrix sub (int min_row, int max_row)
 
HepSymMatrix inverse (int &ifail) const
 
void invert (int &ifail)
 
void invert ()
 
HepSymMatrix inverse () const
 
double determinant () const
 
double trace () const
 
HepSymMatrix_row operator[] (int)
 
HepSymMatrix_row_const operator[] (int) const
 
void invertCholesky5 (int &ifail)
 
void invertCholesky6 (int &ifail)
 
void invertHaywood4 (int &ifail)
 
void invertHaywood5 (int &ifail)
 
void invertHaywood6 (int &ifail)
 
void invertBunchKaufman (int &ifail)
 
- 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

int num_size () const
 
- 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 HepSymMatrix_row
 
class HepSymMatrix_row_const
 
class HepMatrix
 
class HepDiagMatrix
 
void tridiagonal (HepSymMatrix *a, HepMatrix *hsm)
 
double condition (const HepSymMatrix &m)
 
void diag_step (HepSymMatrix *t, int begin, int end)
 
void diag_step (HepSymMatrix *t, HepMatrix *u, int begin, int end)
 
HepMatrix diagonalize (HepSymMatrix *s)
 
HepVector house (const HepSymMatrix &a, int row, int col)
 
void house_with_update2 (HepSymMatrix *a, HepMatrix *v, int row, int col)
 
HepSymMatrix operator+ (const HepSymMatrix &hm1, const HepSymMatrix &hm2)
 
HepSymMatrix operator- (const HepSymMatrix &hm1, const HepSymMatrix &hm2)
 
HepMatrix operator* (const HepSymMatrix &hm1, const HepSymMatrix &hm2)
 
HepMatrix operator* (const HepSymMatrix &hm1, const HepMatrix &hm2)
 
HepMatrix operator* (const HepMatrix &hm1, const HepSymMatrix &hm2)
 
HepVector operator* (const HepSymMatrix &hm1, const HepVector &hm2)
 
HepSymMatrix vT_times_v (const HepVector &v)
 

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 86 of file SymMatrix.h.

Constructor & Destructor Documentation

◆ HepSymMatrix() [1/6]

CLHEP::HepSymMatrix::HepSymMatrix ( )
inline

◆ HepSymMatrix() [2/6]

CLHEP::HepSymMatrix::HepSymMatrix ( int  p)
explicit

Definition at line 56 of file SymMatrix.cc.

57 : m(p*(p+1)/2), nrow(p)
58{
59 size_ = nrow * (nrow+1) / 2;
60 m.assign(size_,0);
61}

◆ HepSymMatrix() [3/6]

CLHEP::HepSymMatrix::HepSymMatrix ( int  p,
int  init 
)

Definition at line 63 of file SymMatrix.cc.

64 : m(p*(p+1)/2), nrow(p)
65{
66 size_ = nrow * (nrow+1) / 2;
67
68 m.assign(size_,0);
69 switch(init)
70 {
71 case 0:
72 break;
73
74 case 1:
75 {
76 HepMatrix::mIter a;
77 for(int i=0;i<nrow;++i) {
78 a = m.begin() + (i+1)*i/2 + i;
79 *a = 1.0;
80 }
81 break;
82 }
83 default:
84 error("SymMatrix: initialization must be either 0 or 1.");
85 }
86}
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

◆ HepSymMatrix() [4/6]

CLHEP::HepSymMatrix::HepSymMatrix ( int  p,
HepRandom &  r 
)

Definition at line 88 of file SymMatrix.cc.

89 : m(p*(p+1)/2), nrow(p)
90{
91 size_ = nrow * (nrow+1) / 2;
92 HepMatrix::mIter a = m.begin();
93 HepMatrix::mIter b = m.begin() + size_;
94 for(;a<b;a++) *a = r();
95}

◆ HepSymMatrix() [5/6]

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

Definition at line 103 of file SymMatrix.cc.

104 : HepGenMatrix(hm1), m(hm1.size_), nrow(hm1.nrow), size_(hm1.size_)
105{
106 m = hm1.m;
107}

◆ HepSymMatrix() [6/6]

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

Definition at line 109 of file SymMatrix.cc.

110 : m(hm1.nrow*(hm1.nrow+1)/2), nrow(hm1.nrow)
111{
112 size_ = nrow * (nrow+1) / 2;
113
114 int n = num_row();
115 m.assign(size_,0);
116
117 HepMatrix::mIter mrr = m.begin();
118 HepMatrix::mcIter mr = hm1.m.begin();
119 for(int r=1;r<=n;r++) {
120 *mrr = *(mr++);
121 if(r<n) mrr += (r+1);
122 }
123}
CLHEP::HepGenMatrix::mcIter
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
Definition: GenMatrix.h:74
CLHEP::HepSymMatrix::num_row
int num_row() const

◆ ~HepSymMatrix()

CLHEP::HepSymMatrix::~HepSymMatrix ( )
virtual

Definition at line 100 of file SymMatrix.cc.

100 {
101}

Member Function Documentation

◆ apply()

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

Definition at line 696 of file SymMatrix.cc.

700{
701#else
702{
703 HepSymMatrix mret(num_row());
704#endif
705 HepMatrix::mcIter a = m.begin();
706 HepMatrix::mIter b = mret.m.begin();
707 for(int ir=1;ir<=num_row();ir++) {
708 for(int ic=1;ic<=ir;ic++) {
709 *(b++) = (*f)(*(a++), ir, ic);
710 }
711 }
712 return mret;
713}
f
void f(void g())
Definition: excDblThrow.cc:38

Referenced by main().

◆ assign() [1/2]

void CLHEP::HepSymMatrix::assign ( const HepMatrix &  hm2)

Definition at line 715 of file SymMatrix.cc.

716{
717 if(hm1.nrow != nrow)
718 {
719 nrow = hm1.nrow;
720 size_ = nrow * (nrow+1) / 2;
721 m.resize(size_);
722 }
723 HepMatrix::mcIter a = hm1.m.begin();
724 HepMatrix::mIter b = m.begin();
725 for(int r=1;r<=nrow;r++) {
726 HepMatrix::mcIter d = a;
727 for(int c=1;c<=r;c++) {
728 *(b++) = *(d++);
729 }
730 if(r<nrow) a += nrow;
731 }
732}

Referenced by main(), and testRandMultiGauss().

◆ assign() [2/2]

void CLHEP::HepSymMatrix::assign ( const HepSymMatrix &  hm2)

◆ determinant()

double CLHEP::HepSymMatrix::determinant ( ) const

Definition at line 940 of file SymMatrix.cc.

940 {
941 static const int max_array = 20;
942 // ir must point to an array which is ***1 longer than*** nrow
943 static std::vector<int> ir_vec (max_array+1);
944 if (ir_vec.size() <= static_cast<unsigned int>(nrow)) ir_vec.resize(nrow+1);
945 int * ir = &ir_vec[0];
946
947 double det;
948 HepMatrix mt(*this);
949 int i = mt.dfact_matrix(det, ir);
950 if(i==0) return det;
951 return 0.0;
952}
CLHEP::HepSymMatrix::HepMatrix
friend class HepMatrix
Definition: SymMatrix.h:240

Referenced by test_inversion().

◆ fast() [1/2]

double & CLHEP::HepSymMatrix::fast ( int  row,
int  col 
)

◆ fast() [2/2]

const double & CLHEP::HepSymMatrix::fast ( int  row,
int  col 
) const

Referenced by main(), and CLHEP::norm().

◆ inverse() [1/2]

HepSymMatrix CLHEP::HepSymMatrix::inverse ( ) const
inline

◆ inverse() [2/2]

HepSymMatrix CLHEP::HepSymMatrix::inverse ( int &  ifail) const
inline

Referenced by main().

◆ invert() [1/2]

void CLHEP::HepSymMatrix::invert ( )
inline

◆ invert() [2/2]

void CLHEP::HepSymMatrix::invert ( int &  ifail)
virtual

Implements CLHEP::HepGenMatrix.

Definition at line 842 of file SymMatrix.cc.

842 {
843
844 ifail = 0;
845
846 switch(nrow) {
847 case 3:
848 {
849 double det, temp;
850 double t1, t2, t3;
851 double c11,c12,c13,c22,c23,c33;
852 c11 = (*(m.begin()+2)) * (*(m.begin()+5)) - (*(m.begin()+4)) * (*(m.begin()+4));
853 c12 = (*(m.begin()+4)) * (*(m.begin()+3)) - (*(m.begin()+1)) * (*(m.begin()+5));
854 c13 = (*(m.begin()+1)) * (*(m.begin()+4)) - (*(m.begin()+2)) * (*(m.begin()+3));
855 c22 = (*(m.begin()+5)) * (*m.begin()) - (*(m.begin()+3)) * (*(m.begin()+3));
856 c23 = (*(m.begin()+3)) * (*(m.begin()+1)) - (*(m.begin()+4)) * (*m.begin());
857 c33 = (*m.begin()) * (*(m.begin()+2)) - (*(m.begin()+1)) * (*(m.begin()+1));
858 t1 = fabs(*m.begin());
859 t2 = fabs(*(m.begin()+1));
860 t3 = fabs(*(m.begin()+3));
861 if (t1 >= t2) {
862 if (t3 >= t1) {
863 temp = *(m.begin()+3);
864 det = c23*c12-c22*c13;
865 } else {
866 temp = *m.begin();
867 det = c22*c33-c23*c23;
868 }
869 } else if (t3 >= t2) {
870 temp = *(m.begin()+3);
871 det = c23*c12-c22*c13;
872 } else {
873 temp = *(m.begin()+1);
874 det = c13*c23-c12*c33;
875 }
876 if (det==0) {
877 ifail = 1;
878 return;
879 }
880 {
881 double ds = temp/det;
882 HepMatrix::mIter hmm = m.begin();
883 *(hmm++) = ds*c11;
884 *(hmm++) = ds*c12;
885 *(hmm++) = ds*c22;
886 *(hmm++) = ds*c13;
887 *(hmm++) = ds*c23;
888 *(hmm) = ds*c33;
889 }
890 }
891 break;
892 case 2:
893 {
894 double det, temp, ds;
895 det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
896 if (det==0) {
897 ifail = 1;
898 return;
899 }
900 ds = 1.0/det;
901 *(m.begin()+1) *= -ds;
902 temp = ds*(*(m.begin()+2));
903 *(m.begin()+2) = ds*(*m.begin());
904 *m.begin() = temp;
905 break;
906 }
907 case 1:
908 {
909 if ((*m.begin())==0) {
910 ifail = 1;
911 return;
912 }
913 *m.begin() = 1.0/(*m.begin());
914 break;
915 }
916 case 5:
917 {
918 invert5(ifail);
919 return;
920 }
921 case 6:
922 {
923 invert6(ifail);
924 return;
925 }
926 case 4:
927 {
928 invert4(ifail);
929 return;
930 }
931 default:
932 {
933 invertBunchKaufman(ifail);
934 return;
935 }
936 }
937 return; // inversion successful
938}
CLHEP::HepSymMatrix::invertBunchKaufman
void invertBunchKaufman(int &ifail)
Definition: SymMatrix.cc:961

Referenced by test_inversion().

◆ invertBunchKaufman()

void CLHEP::HepSymMatrix::invertBunchKaufman ( int &  ifail)

Definition at line 961 of file SymMatrix.cc.

961 {
962 // Bunch-Kaufman diagonal pivoting method
963 // It is decribed in J.R. Bunch, L. Kaufman (1977).
964 // "Some Stable Methods for Calculating Inertia and Solving Symmetric
965 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
966 // Charles F. van Loan, "Matrix Computations" (the second edition
967 // has a bug.) and implemented in "lapack"
968 // Mario Stanke, 09/97
969
970 int i, j, k, is;
971 int pivrow;
972
973 // Establish the two working-space arrays needed: x and piv are
974 // used as pointers to arrays of doubles and ints respectively, each
975 // of length nrow. We do not want to reallocate each time through
976 // unless the size needs to grow. We do not want to leak memory, even
977 // by having a new without a delete that is only done once.
978
979 static const int max_array = 25;
980#ifdef DISABLE_ALLOC
981 static std::vector<double> xvec (max_array);
982 static std::vector<int> pivv (max_array);
983 typedef std::vector<int>::iterator pivIter;
984#else
985 static std::vector<double,Alloc<double,25> > xvec (max_array);
986 static std::vector<int, Alloc<int, 25> > pivv (max_array);
987 typedef std::vector<int,Alloc<int,25> >::iterator pivIter;
988#endif
989 if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
990 if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
991 // Note - resize shuld do nothing if the size is already larger than nrow,
992 // but on VC++ there are indications that it does so we check.
993 // Note - the data elements in a vector are guaranteed to be contiguous,
994 // so x[i] and piv[i] are optimally fast.
995 mIter x = xvec.begin();
996 // x[i] is used as helper storage, needs to have at least size nrow.
997 pivIter piv = pivv.begin();
998 // piv[i] is used to store details of exchanges
999
1000 double temp1, temp2;
1001 HepMatrix::mIter ip, mjj, iq;
1002 double lambda, sigma;
1003 const double alpha = .6404; // = (1+sqrt(17))/8
1004 const double epsilon = 32*DBL_EPSILON;
1005 // whenever a sum of two doubles is below or equal to epsilon
1006 // it is set to zero.
1007 // this constant could be set to zero but then the algorithm
1008 // doesn't neccessarily detect that a matrix is singular
1009
1010 for (i = 0; i < nrow; ++i) piv[i] = i+1;
1011
1012 ifail = 0;
1013
1014 // compute the factorization P*A*P^T = L * D * L^T
1015 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
1016 // L and D^-1 are stored in A = *this, P is stored in piv[]
1017
1018 for (j=1; j < nrow; j+=is) // main loop over columns
1019 {
1020 mjj = m.begin() + j*(j-1)/2 + j-1;
1021 lambda = 0; // compute lambda = max of A(j+1:n,j)
1022 pivrow = j+1;
1023 //ip = m.begin() + (j+1)*j/2 + j-1;
1024 for (i=j+1; i <= nrow ; ++i) {
1025 // calculate ip to avoid going off end of storage array
1026 ip = m.begin() + (i-1)*i/2 + j-1;
1027 if (fabs(*ip) > lambda) {
1028 lambda = fabs(*ip);
1029 pivrow = i;
1030 }
1031 } // for i
1032 if (lambda == 0 ) {
1033 if (*mjj == 0) {
1034 ifail = 1;
1035 return;
1036 }
1037 is=1;
1038 *mjj = 1./ *mjj;
1039 } else { // lambda == 0
1040 if (fabs(*mjj) >= lambda*alpha) {
1041 is=1;
1042 pivrow=j;
1043 } else { // fabs(*mjj) >= lambda*alpha
1044 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
1045 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
1046 for (k=j; k < pivrow; k++) {
1047 if (fabs(*ip) > sigma) sigma = fabs(*ip);
1048 ip++;
1049 } // for k
1050 if (sigma * fabs(*mjj) >= alpha * lambda * lambda) {
1051 is=1;
1052 pivrow = j;
1053 } else if (fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
1054 >= alpha * sigma) {
1055 is=1;
1056 } else {
1057 is=2;
1058 } // if sigma...
1059 } // fabs(*mjj) >= lambda*alpha
1060 if (pivrow == j) { // no permutation neccessary
1061 piv[j-1] = pivrow;
1062 if (*mjj == 0) {
1063 ifail=1;
1064 return;
1065 }
1066 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
1067
1068 // update A(j+1:n, j+1,n)
1069 for (i=j+1; i <= nrow; i++) {
1070 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
1071 ip = m.begin()+i*(i-1)/2+j;
1072 for (k=j+1; k<=i; k++) {
1073 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
1074 if (fabs(*ip) <= epsilon)
1075 *ip=0;
1076 ip++;
1077 }
1078 } // for i
1079 // update L
1080 //ip = m.begin() + (j+1)*j/2 + j-1;
1081 for (i=j+1; i <= nrow; ++i) {
1082 // calculate ip to avoid going off end of storage array
1083 ip = m.begin() + (i-1)*i/2 + j-1;
1084 *ip *= temp2;
1085 }
1086 } else if (is==1) { // 1x1 pivot
1087 piv[j-1] = pivrow;
1088
1089 // interchange rows and columns j and pivrow in
1090 // submatrix (j:n,j:n)
1091 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1092 for (i=j+1; i < pivrow; i++, ip++) {
1093 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1094 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
1095 *ip = temp1;
1096 } // for i
1097 temp1 = *mjj;
1098 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
1099 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
1100 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
1101 iq = ip + pivrow-j;
1102 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) {
1103 temp1 = *iq;
1104 *iq = *ip;
1105 *ip = temp1;
1106 } // for i
1107
1108 if (*mjj == 0) {
1109 ifail = 1;
1110 return;
1111 } // *mjj == 0
1112 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
1113
1114 // update A(j+1:n, j+1:n)
1115 for (i = j+1; i <= nrow; i++) {
1116 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
1117 ip = m.begin()+i*(i-1)/2+j;
1118 for (k=j+1; k<=i; k++) {
1119 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
1120 if (fabs(*ip) <= epsilon)
1121 *ip=0;
1122 ip++;
1123 } // for k
1124 } // for i
1125 // update L
1126 //ip = m.begin() + (j+1)*j/2 + j-1;
1127 for (i=j+1; i <= nrow; ++i) {
1128 // calculate ip to avoid going off end of storage array
1129 ip = m.begin() + (i-1)*i/2 + j-1;
1130 *ip *= temp2;
1131 }
1132 } else { // is=2, ie use a 2x2 pivot
1133 piv[j-1] = -pivrow;
1134 piv[j] = 0; // that means this is the second row of a 2x2 pivot
1135
1136 if (j+1 != pivrow) {
1137 // interchange rows and columns j+1 and pivrow in
1138 // submatrix (j:n,j:n)
1139 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
1140 for (i=j+2; i < pivrow; i++, ip++) {
1141 temp1 = *(m.begin() + i*(i-1)/2 + j);
1142 *(m.begin() + i*(i-1)/2 + j) = *ip;
1143 *ip = temp1;
1144 } // for i
1145 temp1 = *(mjj + j + 1);
1146 *(mjj + j + 1) =
1147 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1148 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1149 temp1 = *(mjj + j);
1150 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
1151 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
1152 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
1153 iq = ip + pivrow-(j+1);
1154 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) {
1155 temp1 = *iq;
1156 *iq = *ip;
1157 *ip = temp1;
1158 } // for i
1159 } // j+1 != pivrow
1160 // invert D(j:j+1,j:j+1)
1161 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
1162 if (temp2 == 0) {
1163 std::cerr
1164 << "SymMatrix::bunch_invert: error in pivot choice"
1165 << std::endl;
1166 }
1167 temp2 = 1. / temp2;
1168 // this quotient is guaranteed to exist by the choice
1169 // of the pivot
1170 temp1 = *mjj;
1171 *mjj = *(mjj + j + 1) * temp2;
1172 *(mjj + j + 1) = temp1 * temp2;
1173 *(mjj + j) = - *(mjj + j) * temp2;
1174
1175 if (j < nrow-1) { // otherwise do nothing
1176 // update A(j+2:n, j+2:n)
1177 for (i=j+2; i <= nrow ; i++) {
1178 ip = m.begin() + i*(i-1)/2 + j-1;
1179 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
1180 if (fabs(temp1 ) <= epsilon) temp1 = 0;
1181 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
1182 if (fabs(temp2 ) <= epsilon) temp2 = 0;
1183 for (k = j+2; k <= i ; k++) {
1184 ip = m.begin() + i*(i-1)/2 + k-1;
1185 iq = m.begin() + k*(k-1)/2 + j-1;
1186 *ip -= temp1 * *iq + temp2 * *(iq+1);
1187 if (fabs(*ip) <= epsilon)
1188 *ip = 0;
1189 } // for k
1190 } // for i
1191 // update L
1192 for (i=j+2; i <= nrow ; i++) {
1193 ip = m.begin() + i*(i-1)/2 + j-1;
1194 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
1195 if (fabs(temp1) <= epsilon) temp1 = 0;
1196 *(ip+1) = *ip * *(mjj + j)
1197 + *(ip+1) * *(mjj + j + 1);
1198 if (fabs(*(ip+1)) <= epsilon) *(ip+1) = 0;
1199 *ip = temp1;
1200 } // for k
1201 } // j < nrow-1
1202 }
1203 }
1204 } // end of main loop over columns
1205
1206 if (j == nrow) { // the the last pivot is 1x1
1207 mjj = m.begin() + j*(j-1)/2 + j-1;
1208 if (*mjj == 0) {
1209 ifail = 1;
1210 return;
1211 } else { *mjj = 1. / *mjj; }
1212 } // end of last pivot code
1213
1214 // computing the inverse from the factorization
1215
1216 for (j = nrow ; j >= 1 ; j -= is) // loop over columns
1217 {
1218 mjj = m.begin() + j*(j-1)/2 + j-1;
1219 if (piv[j-1] > 0) { // 1x1 pivot, compute column j of inverse
1220 is = 1;
1221 if (j < nrow) {
1222 //ip = m.begin() + (j+1)*j/2 + j-1;
1223 //for (i=0; i < nrow-j; ip += 1+j+i++) x[i] = *ip;
1224 ip = m.begin() + (j+1)*j/2 - 1;
1225 for (i=0; i < nrow-j; ++i) {
1226 ip += i + j;
1227 x[i] = *ip;
1228 }
1229 for (i=j+1; i<=nrow ; i++) {
1230 temp2=0;
1231 ip = m.begin() + i*(i-1)/2 + j;
1232 for (k=0; k <= i-j-1; k++) temp2 += *ip++ * x[k];
1233 // avoid setting ip outside the bounds of the storage array
1234 ip -= 1;
1235 // using the value of k from the previous loop
1236 for ( ; k < nrow-j; ++k) {
1237 ip += j+k;
1238 temp2 += *ip * x[k];
1239 }
1240 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1241 } // for i
1242 temp2 = 0;
1243 //ip = m.begin() + (j+1)*j/2 + j-1;
1244 //for (k=0; k < nrow-j; ip += 1+j+k++)
1245 //temp2 += x[k] * *ip;
1246 ip = m.begin() + (j+1)*j/2 - 1;
1247 for (k=0; k < nrow-j; ++k) {
1248 ip += j+k;
1249 temp2 += x[k] * *ip;
1250 }
1251 *mjj -= temp2;
1252 } // j < nrow
1253 } else { //2x2 pivot, compute columns j and j-1 of the inverse
1254 if (piv[j-1] != 0)
1255 std::cerr << "error in piv" << piv[j-1] << std::endl;
1256 is=2;
1257 if (j < nrow) {
1258 //ip = m.begin() + (j+1)*j/2 + j-1;
1259 //for (i=0; i < nrow-j; ip += 1+j+i++) x[i] = *ip;
1260 ip = m.begin() + (j+1)*j/2 - 1;
1261 for (i=0; i < nrow-j; ++i) {
1262 ip += i + j;
1263 x[i] = *ip;
1264 }
1265 for (i=j+1; i<=nrow ; i++) {
1266 temp2 = 0;
1267 ip = m.begin() + i*(i-1)/2 + j;
1268 for (k=0; k <= i-j-1; k++)
1269 temp2 += *ip++ * x[k];
1270 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1271 temp2 += *ip * x[k];
1272 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1273 } // for i
1274 temp2 = 0;
1275 //ip = m.begin() + (j+1)*j/2 + j-1;
1276 //for (k=0; k < nrow-j; ip += 1+j+k++) temp2 += x[k] * *ip;
1277 ip = m.begin() + (j+1)*j/2 - 1;
1278 for (k=0; k < nrow-j; ++k) {
1279 ip += k + j;
1280 temp2 += x[k] * *ip;
1281 }
1282 *mjj -= temp2;
1283 temp2 = 0;
1284 //ip = m.begin() + (j+1)*j/2 + j-2;
1285 //for (i=j+1; i <= nrow; ip += i++) temp2 += *ip * *(ip+1);
1286 ip = m.begin() + (j+1)*j/2 - 2;
1287 for (i=j+1; i <= nrow; ++i) {
1288 ip += i - 1;
1289 temp2 += *ip * *(ip+1);
1290 }
1291 *(mjj-1) -= temp2;
1292 //ip = m.begin() + (j+1)*j/2 + j-2;
1293 //for (i=0; i < nrow-j; ip += 1+j+i++) x[i] = *ip;
1294 ip = m.begin() + (j+1)*j/2 - 2;
1295 for (i=0; i < nrow-j; ++i) {
1296 ip += i + j;
1297 x[i] = *ip;
1298 }
1299 for (i=j+1; i <= nrow ; i++) {
1300 temp2 = 0;
1301 ip = m.begin() + i*(i-1)/2 + j;
1302 for (k=0; k <= i-j-1; k++)
1303 temp2 += *ip++ * x[k];
1304 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1305 temp2 += *ip * x[k];
1306 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1307 } // for i
1308 temp2 = 0;
1309 //ip = m.begin() + (j+1)*j/2 + j-2;
1310 //for (k=0; k < nrow-j; ip += 1+j+k++)
1311 // temp2 += x[k] * *ip;
1312 ip = m.begin() + (j+1)*j/2 - 2;
1313 for (k=0; k < nrow-j; ++k) {
1314 ip += k + j;
1315 temp2 += x[k] * *ip;
1316 }
1317 *(mjj-j) -= temp2;
1318 } // j < nrow
1319 } // else piv[j-1] > 0
1320
1321 // interchange rows and columns j and piv[j-1]
1322 // or rows and columns j and -piv[j-2]
1323
1324 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1325 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1326 for (i=j+1;i < pivrow; i++, ip++) {
1327 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1328 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1329 *ip = temp1;
1330 } // for i
1331 temp1 = *mjj;
1332 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1333 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1334 if (is==2) {
1335 temp1 = *(mjj-1);
1336 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1337 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1338 } // is==2
1339
1340 // problem right here
1341 if( pivrow < nrow ) {
1342 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1343 // adding parenthesis for VC++
1344 iq = ip + (pivrow-j);
1345 for (i = pivrow+1; i <= nrow; i++) {
1346 temp1 = *iq;
1347 *iq = *ip;
1348 *ip = temp1;
1349 if( i < nrow ) {
1350 ip += i;
1351 iq += i;
1352 }
1353 } // for i
1354 } // pivrow < nrow
1355 } // end of loop over columns (in computing inverse from factorization)
1356
1357 return; // inversion successful
1358
1359}

Referenced by invert().

◆ invertCholesky5()

void CLHEP::HepSymMatrix::invertCholesky5 ( int &  ifail)

Definition at line 679 of file SymMatrixInvert.cc.

679 {
680
681// Invert by
682//
683// a) decomposing M = G*G^T with G lower triangular
684// (if M is not positive definite this will fail, leaving this unchanged)
685// b) inverting G to form H
686// c) multiplying H^T * H to get M^-1.
687//
688// If the matrix is pos. def. it is inverted and 1 is returned.
689// If the matrix is not pos. def. it remains unaltered and 0 is returned.
690
691 double h10; // below-diagonal elements of H
692 double h20, h21;
693 double h30, h31, h32;
694 double h40, h41, h42, h43;
695
696 double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
697 // diagonal elements of H
698
699 double g10; // below-diagonal elements of G
700 double g20, g21;
701 double g30, g31, g32;
702 double g40, g41, g42, g43;
703
704 ifail = 1; // We start by assuing we won't succeed...
705
706// Form G -- compute diagonal members of H directly rather than of G
707//-------
708
709// Scale first column by 1/sqrt(A00)
710
711 h00 = m[A00];
712 if (h00 <= 0) return;
713 h00 = 1.0 / sqrt(h00);
714
715 g10 = m[A10] * h00;
716 g20 = m[A20] * h00;
717 g30 = m[A30] * h00;
718 g40 = m[A40] * h00;
719
720// Form G11 (actually, just h11)
721
722 h11 = m[A11] - (g10 * g10);
723 if (h11 <= 0) return;
724 h11 = 1.0 / sqrt(h11);
725
726// Subtract inter-column column dot products from rest of column 1 and
727// scale to get column 1 of G
728
729 g21 = (m[A21] - (g10 * g20)) * h11;
730 g31 = (m[A31] - (g10 * g30)) * h11;
731 g41 = (m[A41] - (g10 * g40)) * h11;
732
733// Form G22 (actually, just h22)
734
735 h22 = m[A22] - (g20 * g20) - (g21 * g21);
736 if (h22 <= 0) return;
737 h22 = 1.0 / sqrt(h22);
738
739// Subtract inter-column column dot products from rest of column 2 and
740// scale to get column 2 of G
741
742 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
743 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
744
745// Form G33 (actually, just h33)
746
747 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
748 if (h33 <= 0) return;
749 h33 = 1.0 / sqrt(h33);
750
751// Subtract inter-column column dot product from A43 and scale to get G43
752
753 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
754
755// Finally form h44 - if this is possible inversion succeeds
756
757 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
758 if (h44 <= 0) return;
759 h44 = 1.0 / sqrt(h44);
760
761// Form H = 1/G -- diagonal members of H are already correct
762//-------------
763
764// The order here is dictated by speed considerations
765
766 h43 = -h33 * g43 * h44;
767 h32 = -h22 * g32 * h33;
768 h42 = -h22 * (g32 * h43 + g42 * h44);
769 h21 = -h11 * g21 * h22;
770 h31 = -h11 * (g21 * h32 + g31 * h33);
771 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
772 h10 = -h00 * g10 * h11;
773 h20 = -h00 * (g10 * h21 + g20 * h22);
774 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
775 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
776
777// Change this to its inverse = H^T*H
778//------------------------------------
779
780 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
781 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
782 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
783 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
784 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
785 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
786 m[A03] = h30 * h33 + h40 * h43;
787 m[A13] = h31 * h33 + h41 * h43;
788 m[A23] = h32 * h33 + h42 * h43;
789 m[A33] = h33 * h33 + h43 * h43;
790 m[A04] = h40 * h44;
791 m[A14] = h41 * h44;
792 m[A24] = h42 * h44;
793 m[A34] = h43 * h44;
794 m[A44] = h44 * h44;
795
796 ifail = 0;
797 return;
798
799}
A41
#define A41
Definition: MatrixInvert.cc:48
A20
#define A20
Definition: MatrixInvert.cc:33
A01
#define A01
Definition: MatrixInvert.cc:20
A23
#define A23
Definition: MatrixInvert.cc:36
A13
#define A13
Definition: MatrixInvert.cc:29
A00
#define A00
Definition: MatrixInvert.cc:19
A40
#define A40
Definition: MatrixInvert.cc:47
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
A03
#define A03
Definition: MatrixInvert.cc:22
A14
#define A14
Definition: MatrixInvert.cc:30
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
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
A34
#define A34
Definition: MatrixInvert.cc:44
A43
#define A43
Definition: MatrixInvert.cc:50

◆ invertCholesky6()

void CLHEP::HepSymMatrix::invertCholesky6 ( int &  ifail)

Definition at line 802 of file SymMatrixInvert.cc.

802 {
803
804// Invert by
805//
806// a) decomposing M = G*G^T with G lower triangular
807// (if M is not positive definite this will fail, leaving this unchanged)
808// b) inverting G to form H
809// c) multiplying H^T * H to get M^-1.
810//
811// If the matrix is pos. def. it is inverted and 1 is returned.
812// If the matrix is not pos. def. it remains unaltered and 0 is returned.
813
814 double h10; // below-diagonal elements of H
815 double h20, h21;
816 double h30, h31, h32;
817 double h40, h41, h42, h43;
818 double h50, h51, h52, h53, h54;
819
820 double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G =
821 // diagonal elements of H
822
823 double g10; // below-diagonal elements of G
824 double g20, g21;
825 double g30, g31, g32;
826 double g40, g41, g42, g43;
827 double g50, g51, g52, g53, g54;
828
829 ifail = 1; // We start by assuing we won't succeed...
830
831// Form G -- compute diagonal members of H directly rather than of G
832//-------
833
834// Scale first column by 1/sqrt(A00)
835
836 h00 = m[A00];
837 if (h00 <= 0) return;
838 h00 = 1.0 / sqrt(h00);
839
840 g10 = m[A10] * h00;
841 g20 = m[A20] * h00;
842 g30 = m[A30] * h00;
843 g40 = m[A40] * h00;
844 g50 = m[A50] * h00;
845
846// Form G11 (actually, just h11)
847
848 h11 = m[A11] - (g10 * g10);
849 if (h11 <= 0) return;
850 h11 = 1.0 / sqrt(h11);
851
852// Subtract inter-column column dot products from rest of column 1 and
853// scale to get column 1 of G
854
855 g21 = (m[A21] - (g10 * g20)) * h11;
856 g31 = (m[A31] - (g10 * g30)) * h11;
857 g41 = (m[A41] - (g10 * g40)) * h11;
858 g51 = (m[A51] - (g10 * g50)) * h11;
859
860// Form G22 (actually, just h22)
861
862 h22 = m[A22] - (g20 * g20) - (g21 * g21);
863 if (h22 <= 0) return;
864 h22 = 1.0 / sqrt(h22);
865
866// Subtract inter-column column dot products from rest of column 2 and
867// scale to get column 2 of G
868
869 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
870 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
871 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
872
873// Form G33 (actually, just h33)
874
875 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
876 if (h33 <= 0) return;
877 h33 = 1.0 / sqrt(h33);
878
879// Subtract inter-column column dot products from rest of column 3 and
880// scale to get column 3 of G
881
882 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
883 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
884
885// Form G44 (actually, just h44)
886
887 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
888 if (h44 <= 0) return;
889 h44 = 1.0 / sqrt(h44);
890
891// Subtract inter-column column dot product from M54 and scale to get G54
892
893 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
894
895// Finally form h55 - if this is possible inversion succeeds
896
897 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
898 if (h55 <= 0) return;
899 h55 = 1.0 / sqrt(h55);
900
901// Form H = 1/G -- diagonal members of H are already correct
902//-------------
903
904// The order here is dictated by speed considerations
905
906 h54 = -h44 * g54 * h55;
907 h43 = -h33 * g43 * h44;
908 h53 = -h33 * (g43 * h54 + g53 * h55);
909 h32 = -h22 * g32 * h33;
910 h42 = -h22 * (g32 * h43 + g42 * h44);
911 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
912 h21 = -h11 * g21 * h22;
913 h31 = -h11 * (g21 * h32 + g31 * h33);
914 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
915 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
916 h10 = -h00 * g10 * h11;
917 h20 = -h00 * (g10 * h21 + g20 * h22);
918 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
919 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
920 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
921
922// Change this to its inverse = H^T*H
923//------------------------------------
924
925 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
926 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
927 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
928 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
929 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
930 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
931 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
932 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
933 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
934 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
935 m[A04] = h40 * h44 + h50 * h54;
936 m[A14] = h41 * h44 + h51 * h54;
937 m[A24] = h42 * h44 + h52 * h54;
938 m[A34] = h43 * h44 + h53 * h54;
939 m[A44] = h44 * h44 + h54 * h54;
940 m[A05] = h50 * h55;
941 m[A15] = h51 * h55;
942 m[A25] = h52 * h55;
943 m[A35] = h53 * h55;
944 m[A45] = h54 * h55;
945 m[A55] = h55 * h55;
946
947 ifail = 0;
948 return;
949
950}
A53
#define A53
Definition: MatrixInvert.cc:57
A45
#define A45
Definition: MatrixInvert.cc:52
A54
#define A54
Definition: MatrixInvert.cc:58
A25
#define A25
Definition: MatrixInvert.cc:38
A55
#define A55
Definition: MatrixInvert.cc:59
A35
#define A35
Definition: MatrixInvert.cc:45
A50
#define A50
Definition: MatrixInvert.cc:54
A51
#define A51
Definition: MatrixInvert.cc:55
A05
#define A05
Definition: MatrixInvert.cc:24
A52
#define A52
Definition: MatrixInvert.cc:56
A15
#define A15
Definition: MatrixInvert.cc:31

◆ invertHaywood4()

void CLHEP::HepSymMatrix::invertHaywood4 ( int &  ifail)

Definition at line 1034 of file SymMatrixInvert.cc.

1034 {
1035 invert4(ifail); // For the 4x4 case, the method we use for invert is already
1036 // the Haywood method.
1037}

◆ invertHaywood5()

void CLHEP::HepSymMatrix::invertHaywood5 ( int &  ifail)

Definition at line 120 of file SymMatrixInvert.cc.

120 {
121
122 ifail = 0;
123
124 // Find all NECESSARY 2x2 dets: (25 of them)
125
126 double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
127 double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
128 double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
129 double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
130 double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
131 double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
132 double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
133 double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
134 double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
135 double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
136 double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
137 double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
138 double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
139 double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
140 double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
141 double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
142 double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
143 double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
144 double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
145 double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
146 double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
147 double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
148 double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
149 double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
150 double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
151
152 // Find all NECESSARY 3x3 dets: (30 of them)
153
154 double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
155 + m[A12]*Det2_23_01;
156 double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
157 + m[A13]*Det2_23_01;
158 double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
159 + m[A13]*Det2_23_02;
160 double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
161 + m[A13]*Det2_23_12;
162 double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02
163 + m[A12]*Det2_24_01;
164 double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03
165 + m[A13]*Det2_24_01;
166 double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04
167 + m[A14]*Det2_24_01;
168 double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03
169 + m[A13]*Det2_24_02;
170 double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04
171 + m[A14]*Det2_24_02;
172 double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13
173 + m[A13]*Det2_24_12;
174 double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14
175 + m[A14]*Det2_24_12;
176 double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02
177 + m[A12]*Det2_34_01;
178 double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03
179 + m[A13]*Det2_34_01;
180 double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04
181 + m[A14]*Det2_34_01;
182 double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03
183 + m[A13]*Det2_34_02;
184 double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04
185 + m[A14]*Det2_34_02;
186 double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04
187 + m[A14]*Det2_34_03;
188 double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13
189 + m[A13]*Det2_34_12;
190 double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14
191 + m[A14]*Det2_34_12;
192 double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14
193 + m[A14]*Det2_34_13;
194 double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
195 + m[A22]*Det2_34_01;
196 double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
197 + m[A23]*Det2_34_01;
198 double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
199 + m[A24]*Det2_34_01;
200 double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
201 + m[A23]*Det2_34_02;
202 double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
203 + m[A24]*Det2_34_02;
204 double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
205 + m[A24]*Det2_34_03;
206 double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
207 + m[A23]*Det2_34_12;
208 double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
209 + m[A24]*Det2_34_12;
210 double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
211 + m[A24]*Det2_34_13;
212 double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
213 + m[A24]*Det2_34_23;
214
215 // Find all NECESSARY 4x4 dets: (15 of them)
216
217 double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023
218 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
219 double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023
220 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
221 double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024
222 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
223 double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023
224 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
225 double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024
226 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
227 double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034
228 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
229 double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023
230 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
231 double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024
232 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
233 double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034
234 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
235 double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034
236 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
237 double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
238 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
239 double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
240 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
241 double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
242 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
243 double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
244 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
245 double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
246 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
247
248 // Find the 5x5 det:
249
250 double det = m[A00]*Det4_1234_1234
251 - m[A01]*Det4_1234_0234
252 + m[A02]*Det4_1234_0134
253 - m[A03]*Det4_1234_0124
254 + m[A04]*Det4_1234_0123;
255
256 if ( det == 0 ) {
257#ifdef SINGULAR_DIAGNOSTICS
258 std::cerr << "Kramer's rule inversion of a singular 5x5 matrix: "
259 << *this << "\n";
260#endif
261 ifail = 1;
262 return;
263 }
264
265 double oneOverDet = 1.0/det;
266 double mn1OverDet = - oneOverDet;
267
268 m[A00] = Det4_1234_1234 * oneOverDet;
269 m[A01] = Det4_1234_0234 * mn1OverDet;
270 m[A02] = Det4_1234_0134 * oneOverDet;
271 m[A03] = Det4_1234_0124 * mn1OverDet;
272 m[A04] = Det4_1234_0123 * oneOverDet;
273
274 m[A11] = Det4_0234_0234 * oneOverDet;
275 m[A12] = Det4_0234_0134 * mn1OverDet;
276 m[A13] = Det4_0234_0124 * oneOverDet;
277 m[A14] = Det4_0234_0123 * mn1OverDet;
278
279 m[A22] = Det4_0134_0134 * oneOverDet;
280 m[A23] = Det4_0134_0124 * mn1OverDet;
281 m[A24] = Det4_0134_0123 * oneOverDet;
282
283 m[A33] = Det4_0124_0124 * oneOverDet;
284 m[A34] = Det4_0124_0123 * mn1OverDet;
285
286 m[A44] = Det4_0123_0123 * oneOverDet;
287
288 return;
289}

◆ invertHaywood6()

void CLHEP::HepSymMatrix::invertHaywood6 ( int &  ifail)

Definition at line 291 of file SymMatrixInvert.cc.

291 {
292
293 ifail = 0;
294
295 // Find all NECESSARY 2x2 dets: (39 of them)
296
297 double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
298 double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
299 double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
300 double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
301 double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
302 double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
303 double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
304 double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
305 double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
306 double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
307 double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
308 double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
309 double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
310 double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
311 double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
312 double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
313 double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
314 double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
315 double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
316 double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
317 double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
318 double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
319 double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
320 double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
321 double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
322 double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
323 double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
324 double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
325 double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
326 double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
327 double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
328 double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
329 double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
330 double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
331 double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
332 double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
333 double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
334 double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
335 double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
336
337 // Find all NECESSARY 3x3 dets: (65 of them)
338
339 double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
340 + m[A22]*Det2_34_01;
341 double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
342 + m[A23]*Det2_34_01;
343 double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
344 + m[A24]*Det2_34_01;
345 double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
346 + m[A23]*Det2_34_02;
347 double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
348 + m[A24]*Det2_34_02;
349 double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
350 + m[A24]*Det2_34_03;
351 double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
352 + m[A23]*Det2_34_12;
353 double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
354 + m[A24]*Det2_34_12;
355 double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
356 + m[A24]*Det2_34_13;
357 double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
358 + m[A24]*Det2_34_23;
359 double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
360 + m[A22]*Det2_35_01;
361 double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
362 + m[A23]*Det2_35_01;
363 double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
364 + m[A24]*Det2_35_01;
365 double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
366 + m[A25]*Det2_35_01;
367 double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
368 + m[A23]*Det2_35_02;
369 double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
370 + m[A24]*Det2_35_02;
371 double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
372 + m[A25]*Det2_35_02;
373 double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
374 + m[A24]*Det2_35_03;
375 double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
376 + m[A25]*Det2_35_03;
377 double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
378 + m[A23]*Det2_35_12;
379 double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
380 + m[A24]*Det2_35_12;
381 double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
382 + m[A25]*Det2_35_12;
383 double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
384 + m[A24]*Det2_35_13;
385 double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
386 + m[A25]*Det2_35_13;
387 double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
388 + m[A24]*Det2_35_23;
389 double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
390 + m[A25]*Det2_35_23;
391 double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
392 + m[A22]*Det2_45_01;
393 double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
394 + m[A23]*Det2_45_01;
395 double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
396 + m[A24]*Det2_45_01;
397 double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
398 + m[A25]*Det2_45_01;
399 double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
400 + m[A23]*Det2_45_02;
401 double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
402 + m[A24]*Det2_45_02;
403 double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
404 + m[A25]*Det2_45_02;
405 double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
406 + m[A24]*Det2_45_03;
407 double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
408 + m[A25]*Det2_45_03;
409 double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
410 + m[A25]*Det2_45_04;
411 double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
412 + m[A23]*Det2_45_12;
413 double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
414 + m[A24]*Det2_45_12;
415 double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
416 + m[A25]*Det2_45_12;
417 double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
418 + m[A24]*Det2_45_13;
419 double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
420 + m[A25]*Det2_45_13;
421 double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
422 + m[A25]*Det2_45_14;
423 double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
424 + m[A24]*Det2_45_23;
425 double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
426 + m[A25]*Det2_45_23;
427 double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
428 + m[A25]*Det2_45_24;
429 double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
430 + m[A32]*Det2_45_01;
431 double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
432 + m[A33]*Det2_45_01;
433 double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
434 + m[A34]*Det2_45_01;
435 double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
436 + m[A35]*Det2_45_01;
437 double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
438 + m[A33]*Det2_45_02;
439 double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
440 + m[A34]*Det2_45_02;
441 double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
442 + m[A35]*Det2_45_02;
443 double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
444 + m[A34]*Det2_45_03;
445 double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
446 + m[A35]*Det2_45_03;
447 double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
448 + m[A35]*Det2_45_04;
449 double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
450 + m[A33]*Det2_45_12;
451 double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
452 + m[A34]*Det2_45_12;
453 double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
454 + m[A35]*Det2_45_12;
455 double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
456 + m[A34]*Det2_45_13;
457 double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
458 + m[A35]*Det2_45_13;
459 double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
460 + m[A35]*Det2_45_14;
461 double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
462 + m[A34]*Det2_45_23;
463 double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
464 + m[A35]*Det2_45_23;
465 double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
466 + m[A35]*Det2_45_24;
467 double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
468 + m[A35]*Det2_45_34;
469
470 // Find all NECESSARY 4x4 dets: (55 of them)
471
472 double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
473 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
474 double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
475 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
476 double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
477 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
478 double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
479 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
480 double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
481 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
482 double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
483 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
484 double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
485 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
486 double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
487 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
488 double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
489 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
490 double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
491 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
492 double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
493 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
494 double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
495 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
496 double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
497 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
498 double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
499 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
500 double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
501 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
502 double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
503 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
504 double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
505 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
506 double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
507 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
508 double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
509 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
510 double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
511 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
512 double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
513 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
514 double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
515 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
516 double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
517 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
518 double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
519 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
520 double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
521 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
522 double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
523 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
524 double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
525 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
526 double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
527 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
528 double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
529 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
530 double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
531 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
532 double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
533 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
534 double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
535 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
536 double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
537 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
538 double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
539 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
540 double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
541 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
542 double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
543 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
544 double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
545 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
546 double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
547 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
548 double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
549 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
550 double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
551 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
552 double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
553 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
554 double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
555 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
556 double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
557 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
558 double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
559 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
560 double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
561 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
562 double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
563 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
564 double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
565 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
566 double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
567 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
568 double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
569 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
570 double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
571 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
572 double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
573 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
574 double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
575 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
576 double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
577 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
578 double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
579 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
580 double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
581 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
582
583 // Find all NECESSARY 5x5 dets: (19 of them)
584
585 double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
586 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
587 double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
588 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
589 double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
590 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
591 double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
592 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
593 double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
594 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
595 double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
596 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
597 double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
598 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
599 double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
600 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
601 double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
602 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
603 double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
604 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
605 double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
606 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
607 double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
608 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
609 double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
610 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
611 double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
612 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
613 double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
614 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
615 double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
616 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
617 double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
618 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
619 double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
620 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
621 double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
622 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
623 double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
624 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
625 double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
626 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
627
628 // Find the determinant
629
630 double det = m[A00]*Det5_12345_12345
631 - m[A01]*Det5_12345_02345
632 + m[A02]*Det5_12345_01345
633 - m[A03]*Det5_12345_01245
634 + m[A04]*Det5_12345_01235
635 - m[A05]*Det5_12345_01234;
636
637 if ( det == 0 ) {
638#ifdef SINGULAR_DIAGNOSTICS
639 std::cerr << "Kramer's rule inversion of a singular 6x6 matrix: "
640 << *this << "\n";
641#endif
642 ifail = 1;
643 return;
644 }
645
646 double oneOverDet = 1.0/det;
647 double mn1OverDet = - oneOverDet;
648
649 m[A00] = Det5_12345_12345*oneOverDet;
650 m[A01] = Det5_12345_02345*mn1OverDet;
651 m[A02] = Det5_12345_01345*oneOverDet;
652 m[A03] = Det5_12345_01245*mn1OverDet;
653 m[A04] = Det5_12345_01235*oneOverDet;
654 m[A05] = Det5_12345_01234*mn1OverDet;
655
656 m[A11] = Det5_02345_02345*oneOverDet;
657 m[A12] = Det5_02345_01345*mn1OverDet;
658 m[A13] = Det5_02345_01245*oneOverDet;
659 m[A14] = Det5_02345_01235*mn1OverDet;
660 m[A15] = Det5_02345_01234*oneOverDet;
661
662 m[A22] = Det5_01345_01345*oneOverDet;
663 m[A23] = Det5_01345_01245*mn1OverDet;
664 m[A24] = Det5_01345_01235*oneOverDet;
665 m[A25] = Det5_01345_01234*mn1OverDet;
666
667 m[A33] = Det5_01245_01245*oneOverDet;
668 m[A34] = Det5_01245_01235*mn1OverDet;
669 m[A35] = Det5_01245_01234*oneOverDet;
670
671 m[A44] = Det5_01235_01235*oneOverDet;
672 m[A45] = Det5_01235_01234*mn1OverDet;
673
674 m[A55] = Det5_01234_01234*oneOverDet;
675
676 return;
677}

◆ num_col()

int CLHEP::HepSymMatrix::num_col ( ) const
inlinevirtual

Implements CLHEP::HepGenMatrix.

Referenced by main(), operator+=(), CLHEP::HepMatrix::operator+=(), operator-=(), CLHEP::HepMatrix::operator-=(), CLHEP::operator<<(), and CLHEP::tridiagonal().

◆ num_row()

int CLHEP::HepSymMatrix::num_row ( ) const
inlinevirtual

Implements CLHEP::HepGenMatrix.

Referenced by CLHEP::HepDiagMatrix::assign(), CLHEP::RandMultiGauss::fire(), CLHEP::RandMultiGauss::fireArray(), HepSymMatrix(), main(), CLHEP::norm(), operator+=(), CLHEP::HepMatrix::operator+=(), operator-=(), CLHEP::HepMatrix::operator-=(), CLHEP::operator<<(), CLHEP::RandMultiGauss::RandMultiGauss(), sub(), and CLHEP::tridiagonal().

◆ num_size()

int CLHEP::HepSymMatrix::num_size ( ) const
inlineprotectedvirtual

Implements CLHEP::HepGenMatrix.

◆ operator()() [1/2]

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

Implements CLHEP::HepGenMatrix.

◆ operator()() [2/2]

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

Implements CLHEP::HepGenMatrix.

◆ operator*=()

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

Definition at line 611 of file SymMatrix.cc.

612{
613 SIMPLE_UOP(*=)
614 return (*this);
615}
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: DiagMatrix.cc:26

◆ operator+=() [1/2]

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

Definition at line 464 of file DiagMatrix.cc.

465{
466 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
467 HepMatrix::mIter a=m.begin();
468 HepMatrix::mcIter b=hm2.m.begin();
469 for(int i=1;i<=num_row();i++) {
470 *a += *(b++);
471 if(i<num_row()) a += (i+1);
472 }
473 return (*this);
474}
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: DiagMatrix.cc:44
CLHEP::HepSymMatrix::num_col
int num_col() const

◆ operator+=() [2/2]

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

Definition at line 575 of file SymMatrix.cc.

576{
577 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
578 SIMPLE_BOP(+=)
579 return (*this);
580}
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: DiagMatrix.cc:31

◆ operator-()

HepSymMatrix CLHEP::HepSymMatrix::operator- ( ) const

Definition at line 211 of file SymMatrix.cc.

214{
215#else
216{
217 HepSymMatrix hm2(nrow);
218#endif
219 HepMatrix::mcIter a=m.begin();
220 HepMatrix::mIter b=hm2.m.begin();
221 HepMatrix::mcIter e=m.begin()+num_size();
222 for(;a<e; a++, b++) (*b) = -(*a);
223 return hm2;
224}
CLHEP::HepSymMatrix::num_size
int num_size() const

◆ operator-=() [1/2]

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

Definition at line 496 of file DiagMatrix.cc.

497{
498 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
499 HepMatrix::mIter a=m.begin();
500 HepMatrix::mcIter b=hm2.m.begin();
501 for(int i=1;i<=num_row();i++) {
502 *a -= *(b++);
503 if(i<num_row()) a += (i+1);
504 }
505 return (*this);
506}

◆ operator-=() [2/2]

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

Definition at line 598 of file SymMatrix.cc.

599{
600 CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),-=);
601 SIMPLE_BOP(-=)
602 return (*this);
603}

◆ operator/=()

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

Definition at line 605 of file SymMatrix.cc.

606{
607 SIMPLE_UOP(/=)
608 return (*this);
609}

◆ operator=() [1/2]

HepSymMatrix & CLHEP::HepSymMatrix::operator= ( const HepDiagMatrix &  hm2)

Definition at line 654 of file SymMatrix.cc.

655{
656 if(hm1.nrow != nrow)
657 {
658 nrow = hm1.nrow;
659 size_ = nrow * (nrow+1) / 2;
660 m.resize(size_);
661 }
662
663 m.assign(size_,0);
664 HepMatrix::mIter mrr = m.begin();
665 HepMatrix::mcIter mr = hm1.m.begin();
666 for(int r=1; r<=nrow; r++) {
667 *mrr = *(mr++);
668 if(r<nrow) mrr += (r+1);
669 }
670 return (*this);
671}

◆ operator=() [2/2]

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

Definition at line 642 of file SymMatrix.cc.

643{
644 if(hm1.nrow != nrow)
645 {
646 nrow = hm1.nrow;
647 size_ = hm1.size_;
648 m.resize(size_);
649 }
650 m = hm1.m;
651 return (*this);
652}

◆ operator[]() [1/2]

HepSymMatrix_row CLHEP::HepSymMatrix::operator[] ( int  )
inline

◆ operator[]() [2/2]

HepSymMatrix_row_const CLHEP::HepSymMatrix::operator[] ( int  ) const
inline

◆ similarity() [1/3]

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

Definition at line 734 of file SymMatrix.cc.

737{
738#else
739{
740 HepSymMatrix mret(hm1.num_row());
741#endif
742 HepMatrix temp = hm1*(*this);
743// If hm1*(*this) has correct dimensions, then so will the hm1.T multiplication.
744// So there is no need to check dimensions again.
745 int n = hm1.num_col();
746 HepMatrix::mIter mr = mret.m.begin();
747 HepMatrix::mIter tempr1 = temp.m.begin();
748 for(int r=1;r<=mret.num_row();r++) {
749 HepMatrix::mcIter hm1c1 = hm1.m.begin();
750 for(int c=1;c<=r;c++) {
751 double tmp = 0.0;
752 HepMatrix::mIter tempri = tempr1;
753 HepMatrix::mcIter hm1ci = hm1c1;
754 for(int i=1;i<=hm1.num_col();i++) {
755 tmp+=(*(tempri++))*(*(hm1ci++));
756 }
757 *(mr++) = tmp;
758 hm1c1 += n;
759 }
760 tempr1 += n;
761 }
762 return mret;
763}

Referenced by main(), and symmatrix_test().

◆ similarity() [2/3]

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

Definition at line 765 of file SymMatrix.cc.

768{
769#else
770{
771 HepSymMatrix mret(hm1.num_row());
772#endif
773 HepMatrix temp = hm1*(*this);
774 int n = hm1.num_col();
775 HepMatrix::mIter mr = mret.m.begin();
776 HepMatrix::mIter tempr1 = temp.m.begin();
777 for(int r=1;r<=mret.num_row();r++) {
778 HepMatrix::mcIter hm1c1 = hm1.m.begin();
779 int c;
780 for(c=1;c<=r;c++) {
781 double tmp = 0.0;
782 HepMatrix::mIter tempri = tempr1;
783 HepMatrix::mcIter hm1ci = hm1c1;
784 int i;
785 for(i=1;i<c;i++) {
786 tmp+=(*(tempri++))*(*(hm1ci++));
787 }
788 for(i=c;i<=hm1.num_col();i++) {
789 tmp+=(*(tempri++))*(*(hm1ci));
790 if(i<hm1.num_col()) hm1ci += i;
791 }
792 *(mr++) = tmp;
793 hm1c1 += c;
794 }
795 tempr1 += n;
796 }
797 return mret;
798}

◆ similarity() [3/3]

double CLHEP::HepSymMatrix::similarity ( const HepVector &  v) const

Definition at line 800 of file SymMatrix.cc.

801 {
802 double mret = 0.0;
803 HepVector temp = (*this) *hm1;
804// If hm1*(*this) has correct dimensions, then so will the hm1.T multiplication.
805// So there is no need to check dimensions again.
806 HepMatrix::mIter a=temp.m.begin();
807 HepMatrix::mcIter b=hm1.m.begin();
808 HepMatrix::mIter e=a+hm1.num_row();
809 for(;a<e;) mret += (*(a++)) * (*(b++));
810 return mret;
811}

◆ similarityT()

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

Definition at line 813 of file SymMatrix.cc.

816{
817#else
818{
819 HepSymMatrix mret(hm1.num_col());
820#endif
821 HepMatrix temp = (*this)*hm1;
822 int n = hm1.num_col();
823 HepMatrix::mIter mrc = mret.m.begin();
824 HepMatrix::mIter temp1r = temp.m.begin();
825 for(int r=1;r<=mret.num_row();r++) {
826 HepMatrix::mcIter m11c = hm1.m.begin();
827 for(int c=1;c<=r;c++) {
828 double tmp = 0.0;
829 for(int i=1;i<=hm1.num_row();i++) {
830 HepMatrix::mIter tempir = temp1r + n*(i-1);
831 HepMatrix::mcIter hm1ic = m11c + n*(i-1);
832 tmp+=(*(tempir))*(*(hm1ic));
833 }
834 *(mrc++) = tmp;
835 m11c++;
836 }
837 temp1r++;
838 }
839 return mret;
840}

Referenced by main(), symmatrix_test(), and testRandMultiGauss().

◆ sub() [1/3]

HepSymMatrix CLHEP::HepSymMatrix::sub ( int  min_row,
int  max_row 
)

Definition at line 154 of file SymMatrix.cc.

155{
156 HepSymMatrix mret(max_row-min_row+1);
157 if(max_row > num_row())
158 error("HepSymMatrix::sub: Index out of range");
159 HepMatrix::mIter a = mret.m.begin();
160 HepMatrix::mIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
161 int rowsize=mret.num_row();
162 for(int irow=1; irow<=rowsize; irow++) {
163 HepMatrix::mIter b = b1;
164 for(int icol=0; icol<irow; ++icol) {
165 *(a++) = *(b++);
166 }
167 if(irow<rowsize) b1 += irow+min_row-1;
168 }
169 return mret;
170}

◆ sub() [2/3]

HepSymMatrix CLHEP::HepSymMatrix::sub ( int  min_row,
int  max_row 
) const

Definition at line 131 of file SymMatrix.cc.

134{
135#else
136{
137 HepSymMatrix mret(max_row-min_row+1);
138#endif
139 if(max_row > num_row())
140 error("HepSymMatrix::sub: Index out of range");
141 HepMatrix::mIter a = mret.m.begin();
142 HepMatrix::mcIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
143 int rowsize=mret.num_row();
144 for(int irow=1; irow<=rowsize; irow++) {
145 HepMatrix::mcIter b = b1;
146 for(int icol=0; icol<irow; ++icol) {
147 *(a++) = *(b++);
148 }
149 if(irow<rowsize) b1 += irow+min_row-1;
150 }
151 return mret;
152}

Referenced by main(), matrix_test2(), and symmatrix_test().

◆ sub() [3/3]

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

Definition at line 172 of file SymMatrix.cc.

173{
174 if(row <1 || row+hm1.num_row()-1 > num_row() )
175 error("HepSymMatrix::sub: Index out of range");
176 HepMatrix::mcIter a = hm1.m.begin();
177 HepMatrix::mIter b1 = m.begin() + (row+2)*(row-1)/2;
178 int rowsize=hm1.num_row();
179 for(int irow=1; irow<=rowsize; ++irow) {
180 HepMatrix::mIter b = b1;
181 for(int icol=0; icol<irow; ++icol) {
182 *(b++) = *(a++);
183 }
184 if(irow<rowsize) b1 += irow+row-1;
185 }
186}

◆ T()

HepSymMatrix CLHEP::HepSymMatrix::T ( ) const

Referenced by main().

◆ trace()

double CLHEP::HepSymMatrix::trace ( ) const

Definition at line 954 of file SymMatrix.cc.

954 {
955 double t = 0.0;
956 for (int i=0; i<nrow; i++)
957 t += *(m.begin() + (i+3)*i/2);
958 return t;
959}

Friends And Related Function Documentation

◆ condition

double condition ( const HepSymMatrix &  m)
friend

Definition at line 201 of file MatrixLinear.cc.

202{
203 HepSymMatrix mcopy=hm;
204 diagonalize(&mcopy);
205 double max,min;
206 max=min=fabs(mcopy(1,1));
207
208 int n = mcopy.num_row();
209 HepMatrix::mIter mii = mcopy.m.begin() + 2;
210 for (int i=2; i<=n; i++) {
211 if (max<fabs(*mii)) max=fabs(*mii);
212 if (min>fabs(*mii)) min=fabs(*mii);
213 if(i<n) mii += i+1;
214 }
215 return max/min;
216}
CLHEP::HepSymMatrix::diagonalize
friend HepMatrix diagonalize(HepSymMatrix *s)
Definition: MatrixLinear.cc:318

◆ diag_step [1/2]

void diag_step ( HepSymMatrix *  t,
HepMatrix *  u,
int  begin,
int  end 
)
friend

Definition at line 270 of file MatrixLinear.cc.

271{
272 double d=(t->fast(end-1,end-1)-t->fast(end,end))/2;
273 double mu=t->fast(end,end)-t->fast(end,end-1)*t->fast(end,end-1)/
274 (d+sign(d)*sqrt(d*d+ t->fast(end,end-1)*t->fast(end,end-1)));
275 double x=t->fast(begin,begin)-mu;
276 double z=t->fast(begin+1,begin);
277 HepMatrix::mIter tkk = t->m.begin() + (begin+2)*(begin-1)/2;
278 HepMatrix::mIter tkp1k = tkk + begin;
279 HepMatrix::mIter tkp2k = tkk + 2 * begin + 1;
280 for (int k=begin;k<=end-1;k++) {
281 double c,ds;
282 givens(x,z,&c,&ds);
283 col_givens(u,c,ds,k,k+1);
284
285 // This is the result of G.T*t*G, making use of the special structure
286 // of t and G. Note that since t is symmetric, only the lower half
287 // needs to be updated. Equations were gotten from maple.
288
289 if (k!=begin) {
290 *(tkk-1)= (*(tkk-1))*c-(*(tkp1k-1))*ds;
291 *(tkp1k-1)=0;
292 }
293 double ap=(*tkk);
294 double bp=(*tkp1k);
295 double aq=(*(tkp1k+1));
296 (*tkk)=ap*c*c-2*c*bp*ds+aq*ds*ds;
297 (*tkp1k)=c*ap*ds+bp*c*c-bp*ds*ds-ds*aq*c;
298 (*(tkp1k+1))=ap*ds*ds+2*c*bp*ds+aq*c*c;
299 if (k<end-1) {
300 double bq=(*(tkp2k+1));
301 (*tkp2k)=-bq*ds;
302 (*(tkp2k+1))=bq*c;
303 x=(*tkp1k);
304 z=(*(tkp2k));
305 tkk += k+1;
306 tkp1k += k+2;
307 }
308 if(k<end-2) tkp2k += k+3;
309 }
310}
CLHEP::givens
void givens(double a, double b, double *c, double *s)
Definition: MatrixLinear.cc:882
CLHEP::col_givens
void col_givens(HepMatrix *A, double c, double s, int k1, int k2, int row_min=1, int row_max=0)
Definition: MatrixLinear.cc:124

◆ diag_step [2/2]

void diag_step ( HepSymMatrix *  t,
int  begin,
int  end 
)
friend

Definition at line 227 of file MatrixLinear.cc.

228{
229 double d=(t->fast(end-1,end-1)-t->fast(end,end))/2;
230 double mu=t->fast(end,end)-t->fast(end,end-1)*t->fast(end,end-1)/
231 (d+sign(d)*sqrt(d*d+ t->fast(end,end-1)*t->fast(end,end-1)));
232 double x=t->fast(begin,begin)-mu;
233 double z=t->fast(begin+1,begin);
234 HepMatrix::mIter tkk = t->m.begin() + (begin+2)*(begin-1)/2;
235 HepMatrix::mIter tkp1k = tkk + begin;
236 HepMatrix::mIter tkp2k = tkk + 2 * begin + 1;
237 for (int k=begin;k<=end-1;k++) {
238 double c,ds;
239 givens(x,z,&c,&ds);
240
241 // This is the result of G.T*t*G, making use of the special structure
242 // of t and G. Note that since t is symmetric, only the lower half
243 // needs to be updated. Equations were gotten from maple.
244
245 if (k!=begin)
246 {
247 *(tkk-1)= *(tkk-1)*c-(*(tkp1k-1))*ds;
248 *(tkp1k-1)=0;
249 }
250 double ap=(*tkk);
251 double bp=(*tkp1k);
252 double aq=(*tkp1k+1);
253 (*tkk)=ap*c*c-2*c*bp*ds+aq*ds*ds;
254 (*tkp1k)=c*ap*ds+bp*c*c-bp*ds*ds-ds*aq*c;
255 (*(tkp1k+1))=ap*ds*ds+2*c*bp*ds+aq*c*c;
256 if (k<end-1)
257 {
258 double bq=(*(tkp2k+1));
259 (*tkp2k)=-bq*ds;
260 (*(tkp2k+1))=bq*c;
261 x=(*tkp1k);
262 z=(*tkp2k);
263 tkk += k+1;
264 tkp1k += k+2;
265 }
266 if(k<end-2) tkp2k += k+3;
267 }
268}

◆ diagonalize

HepMatrix diagonalize ( HepSymMatrix *  s)
friend

Definition at line 318 of file MatrixLinear.cc.

319{
320 const double tolerance = 1e-12;
321 HepMatrix u= tridiagonal(hms);
322 int begin=1;
323 int end=hms->num_row();
324 while(begin!=end)
325 {
326 HepMatrix::mIter sii = hms->m.begin() + (begin+2)*(begin-1)/2;
327 HepMatrix::mIter sip1i = sii + begin;
328 for (int i=begin;i<=end-1;i++) {
329 if (fabs(*sip1i)<=
330 tolerance*(fabs(*sii)+fabs(*(sip1i+1)))) {
331 (*sip1i)=0;
332 }
333 if(i<end-1) {
334 sii += i+1;
335 sip1i += i+2;
336 }
337 }
338 while(begin<end && hms->fast(begin+1,begin) ==0) begin++;
339 while(end>begin && hms->fast(end,end-1) ==0) end--;
340 if (begin!=end)
341 diag_step(hms,&u,begin,end);
342 }
343 return u;
344}
CLHEP::HepSymMatrix::tridiagonal
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
Definition: MatrixLinear.cc:777
CLHEP::HepSymMatrix::diag_step
friend void diag_step(HepSymMatrix *t, int begin, int end)
Definition: MatrixLinear.cc:227
CLHEP::HepSymMatrix::fast
const double & fast(int row, int col) const

◆ HepDiagMatrix

friend class HepDiagMatrix
friend

Definition at line 241 of file SymMatrix.h.

◆ HepMatrix

friend class HepMatrix
friend

Definition at line 240 of file SymMatrix.h.

◆ HepSymMatrix_row

friend class HepSymMatrix_row
friend

Definition at line 238 of file SymMatrix.h.

◆ HepSymMatrix_row_const

friend class HepSymMatrix_row_const
friend

Definition at line 239 of file SymMatrix.h.

◆ house

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

Definition at line 353 of file MatrixLinear.cc.

354{
355 HepVector v(a.num_row()-row+1);
356/* not tested */
357 HepMatrix::mIter vp = v.m.begin();
358 HepMatrix::mcIter aci = a.m.begin() + col * (col - 1) / 2 + row - 1;
359 int i;
360 for (i=row;i<=col;i++) {
361 (*(vp++))=(*(aci++));
362 }
363 for (;i<=a.num_row();i++) {
364 (*(vp++))=(*aci);
365 aci += i;
366 }
367 v(1)+=sign(a(row,col))*v.norm();
368 return v;
369}

◆ 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}
CLHEP::norm
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:54

◆ operator* [1/4]

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 }
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: DiagMatrix.cc:49

◆ operator* [2/4]

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* [3/4]

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* [4/4]

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

Definition at line 507 of file SymMatrix.cc.

510{
511#else
512{
513 HepVector mret(hm1.num_row());
514#endif
515 CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
516 HepMatrix::mcIter sp,snp,vpt;
517 double temp;
518 int step,stept;
519 HepMatrix::mIter vrp=mret.m.begin();
520 for(step=1,snp=hm1.m.begin();step<=hm1.num_row();++step)
521 {
522 sp=snp;
523 vpt=hm2.m.begin();
524 snp+=step;
525 temp=0;
526 while(sp<snp)
527 temp+=*(sp++)*(*(vpt++));
528 if(step<hm1.num_row()) sp+=step-1;
529 for(stept=step+1;stept<=hm1.num_row();stept++)
530 {
531 temp+=*sp*(*(vpt++));
532 if(stept<hm1.num_row()) sp+=stept;
533 }
534 *(vrp++)=temp;
535 } // for(step
536 return mret;
537}

◆ operator+

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

Definition at line 253 of file SymMatrix.cc.

256{
257#else
258{
259 HepSymMatrix mret(hm1.nrow);
260#endif
261 CHK_DIM_1(hm1.nrow, hm2.nrow,+);
262 SIMPLE_TOP(+)
263 return mret;
264}
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: DiagMatrix.cc:37

◆ operator-

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

Definition at line 297 of file SymMatrix.cc.

300{
301#else
302{
303 HepSymMatrix mret(hm1.num_row());
304#endif
305 CHK_DIM_1(hm1.num_row(),hm2.num_row(),-);
306 SIMPLE_TOP(-)
307 return mret;
308}

◆ 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::HepSymMatrix::house_with_update2
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
Definition: MatrixLinear.cc:462

◆ vT_times_v

HepSymMatrix vT_times_v ( const HepVector &  v)
friend

Definition at line 539 of file SymMatrix.cc.

542{
543#else
544{
545 HepSymMatrix mret(v.num_row());
546#endif
547 HepMatrix::mIter mr=mret.m.begin();
548 HepMatrix::mcIter vt1,vt2;
549 for(vt1=v.m.begin();vt1<v.m.begin()+v.num_row();vt1++)
550 for(vt2=v.m.begin();vt2<=vt1;vt2++)
551 *(mr++)=(*vt1)*(*vt2);
552 return mret;
553}
The documentation for this class was generated from the following files: