00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 #ifndef _SYMMatrix_H_
00106 #define _SYMMatrix_H_
00107
00108 #ifdef GNUPRAGMA
00109 #pragma interface
00110 #endif
00111
00112 #include "CLHEP/Matrix/GenMatrix.h"
00113 #ifdef HEP_USE_RANDOM
00114 class HepRandom;
00115 #endif
00116
00117 #ifdef HEP_NO_INLINE_IN_DECLARATION
00118 #define inline
00119 #endif
00120
00121 class HepMatrix;
00122 class HepDiagMatrix;
00123 class HepVector;
00124
00129 class HepSymMatrix : public HepGenMatrix {
00130 public:
00131 inline HepSymMatrix();
00132
00133
00134
00135 explicit HepSymMatrix(int p);
00136 HepSymMatrix(int p, int);
00137
00138
00139
00140
00141 #ifdef HEP_USE_RANDOM
00142 HepSymMatrix(int p, HepRandom &r);
00143 #endif
00144
00145 HepSymMatrix(const HepSymMatrix &m1);
00146
00147
00148 HepSymMatrix(const HepDiagMatrix &m1);
00149
00150
00151 virtual ~HepSymMatrix();
00152
00153
00154 inline int num_row() const;
00155 inline int num_col() const;
00156
00157
00158 const double & operator()(int row, int col) const;
00159 double & operator()(int row, int col);
00160
00161
00162
00163 const double & fast(int row, int col) const;
00164 double & fast(int row, int col);
00165
00166
00167
00168
00169 void assign(const HepMatrix &m2);
00170
00171
00172 void assign(const HepSymMatrix &m2);
00173
00174
00175 HepSymMatrix & operator*=(double t);
00176
00177
00178 HepSymMatrix & operator/=(double t);
00179
00180
00181 HepSymMatrix & operator+=( const HepSymMatrix &m2);
00182 HepSymMatrix & operator+=( const HepDiagMatrix &m2);
00183 HepSymMatrix & operator-=( const HepSymMatrix &m2);
00184 HepSymMatrix & operator-=( const HepDiagMatrix &m2);
00185
00186
00187 HepSymMatrix & operator=( const HepSymMatrix &m2);
00188 HepSymMatrix & operator=( const HepDiagMatrix &m2);
00189
00190
00191 HepSymMatrix operator- () const;
00192
00193
00194 HepSymMatrix T() const;
00195
00196
00197 HepSymMatrix apply(double (*f)(double, int, int)) const;
00198
00199
00200 HepSymMatrix similarity(const HepMatrix &m1) const;
00201 HepSymMatrix similarity(const HepSymMatrix &m1) const;
00202
00203
00204 HepSymMatrix similarityT(const HepMatrix &m1) const;
00205
00206
00207
00208 double similarity(const HepVector &v) const;
00209
00210
00211 HepSymMatrix sub(int min_row, int max_row) const;
00212
00213 void sub(int row, const HepSymMatrix &m1);
00214
00215 #ifdef HEP_CC_NEED_SUB_WITHOUT_CONST
00216 HepSymMatrix sub(int min_row, int max_row);
00217
00218
00219 #endif
00220
00221 inline HepSymMatrix inverse(int &ifail) const;
00222
00223
00224
00225 void invert(int &ifail);
00226
00227
00228
00229
00230
00231 double determinant() const;
00232
00233
00234 double trace() const;
00235
00236
00237 class HepSymMatrix_row {
00238 public:
00239 inline HepSymMatrix_row(HepSymMatrix&,int);
00240 inline double & operator[](int);
00241 private:
00242 HepSymMatrix& _a;
00243 int _r;
00244 };
00245 class HepSymMatrix_row_const {
00246 public:
00247 inline HepSymMatrix_row_const(const HepSymMatrix&,int);
00248 inline const double & operator[](int) const;
00249 private:
00250 const HepSymMatrix& _a;
00251 int _r;
00252 };
00253
00254
00255 inline HepSymMatrix_row operator[] (int);
00256 inline HepSymMatrix_row_const operator[] (int) const;
00257
00258
00259
00260
00261
00262
00263
00264
00265 void invertCholesky5 (int &ifail);
00266 void invertCholesky6 (int &ifail);
00267
00268
00269
00270 void invertHaywood4 (int & ifail);
00271 void invertHaywood5 (int &ifail);
00272 void invertHaywood6 (int &ifail);
00273 void invertBunchKaufman (int &ifail);
00274
00275 protected:
00276 inline int num_size() const;
00277
00278 private:
00279 friend class HepSymMatrix_row;
00280 friend class HepSymMatrix_row_const;
00281 friend class HepMatrix;
00282 friend class HepDiagMatrix;
00283
00284 friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
00285 friend double condition(const HepSymMatrix &m);
00286 friend void diag_step(HepSymMatrix *t,int begin,int end);
00287 friend void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end);
00288 friend HepMatrix diagonalize(HepSymMatrix *s);
00289 friend HepVector house(const HepSymMatrix &a,int row,int col);
00290 friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col);
00291
00292 friend HepSymMatrix operator+(const HepSymMatrix &m1,
00293 const HepSymMatrix &m2);
00294 friend HepSymMatrix operator-(const HepSymMatrix &m1,
00295 const HepSymMatrix &m2);
00296 friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00297 friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00298 friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00299 friend HepVector operator*(const HepSymMatrix &m1, const HepVector &m2);
00300
00301
00302 friend HepSymMatrix vT_times_v(const HepVector &v);
00303
00304
00305 double *m;
00306 int nrow;
00307 int size;
00308
00309 static double posDefFraction5x5;
00310 static double adjustment5x5;
00311 static const double CHOLESKY_THRESHOLD_5x5;
00312 static const double CHOLESKY_CREEP_5x5;
00313
00314 static double posDefFraction6x6;
00315 static double adjustment6x6;
00316 static const double CHOLESKY_THRESHOLD_6x6;
00317 static const double CHOLESKY_CREEP_6x6;
00318
00319 void invert4 (int & ifail);
00320 void invert5 (int & ifail);
00321 void invert6 (int & ifail);
00322
00323 };
00324
00325
00326
00327
00328
00329
00330 HepStd::ostream& operator<<(HepStd::ostream &s, const HepSymMatrix &q);
00331
00332
00333 HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00334 HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00335 HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00336 HepSymMatrix operator*(double t, const HepSymMatrix &s1);
00337 HepSymMatrix operator*(const HepSymMatrix &s1, double t);
00338
00339
00340
00341 HepSymMatrix operator/(const HepSymMatrix &m1, double t);
00342
00343
00344 HepMatrix operator+(const HepMatrix &m1, const HepSymMatrix &s2);
00345 HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &m2);
00346 HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2);
00347
00348
00349 HepMatrix operator-(const HepMatrix &m1, const HepSymMatrix &s2);
00350 HepMatrix operator-(const HepSymMatrix &m1, const HepMatrix &m2);
00351 HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2);
00352
00353
00354 HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2);
00355
00356
00357 double condition(const HepSymMatrix &m);
00358
00359
00360 void diag_step(HepSymMatrix *t, int begin, int end);
00361 void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end);
00362
00363
00364 HepMatrix diagonalize(HepSymMatrix *s);
00365
00366
00367
00368 HepVector house(const HepSymMatrix &a, int row=1, int col=1);
00369 void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1);
00370
00371
00372 void tridiagonal(HepSymMatrix *a, HepMatrix *hsm);
00373 HepMatrix tridiagonal(HepSymMatrix *a);
00374
00375
00376 #ifdef HEP_NO_INLINE_IN_DECLARATION
00377 #undef inline
00378 #endif
00379
00380 #if defined(HEP_SHORT_NAMES)
00381 typedef HepSymMatrix SymMatrix;
00382 #endif
00383
00384 #ifndef HEP_DEBUG_INLINE
00385 #include "CLHEP/Matrix/SymMatrix.icc"
00386 #endif
00387
00388 #endif