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 <vector>
00113 #include "CLHEP/Matrix/GenMatrix.h"
00114 #ifdef HEP_USE_RANDOM
00115 class HepRandom;
00116 #endif
00117
00118 #ifdef HEP_NO_INLINE_IN_DECLARATION
00119 #define inline
00120 #endif
00121
00122 class HepMatrix;
00123 class HepDiagMatrix;
00124 class HepVector;
00125
00130 class HepSymMatrix : public HepGenMatrix {
00131 public:
00132 inline HepSymMatrix();
00133
00134
00135
00136 explicit HepSymMatrix(int p);
00137 HepSymMatrix(int p, int);
00138
00139
00140
00141
00142 #ifdef HEP_USE_RANDOM
00143 HepSymMatrix(int p, HepRandom &r);
00144 #endif
00145
00146 HepSymMatrix(const HepSymMatrix &m1);
00147
00148
00149 HepSymMatrix(const HepDiagMatrix &m1);
00150
00151
00152 virtual ~HepSymMatrix();
00153
00154
00155 inline int num_row() const;
00156 inline int num_col() const;
00157
00158
00159 const double & operator()(int row, int col) const;
00160 double & operator()(int row, int col);
00161
00162
00163
00164 const double & fast(int row, int col) const;
00165 double & fast(int row, int col);
00166
00167
00168
00169
00170 void assign(const HepMatrix &m2);
00171
00172
00173 void assign(const HepSymMatrix &m2);
00174
00175
00176 HepSymMatrix & operator*=(double t);
00177
00178
00179 HepSymMatrix & operator/=(double t);
00180
00181
00182 HepSymMatrix & operator+=( const HepSymMatrix &m2);
00183 HepSymMatrix & operator+=( const HepDiagMatrix &m2);
00184 HepSymMatrix & operator-=( const HepSymMatrix &m2);
00185 HepSymMatrix & operator-=( const HepDiagMatrix &m2);
00186
00187
00188 HepSymMatrix & operator=( const HepSymMatrix &m2);
00189 HepSymMatrix & operator=( const HepDiagMatrix &m2);
00190
00191
00192 HepSymMatrix operator- () const;
00193
00194
00195 HepSymMatrix T() const;
00196
00197
00198 HepSymMatrix apply(double (*f)(double, int, int)) const;
00199
00200
00201 HepSymMatrix similarity(const HepMatrix &m1) const;
00202 HepSymMatrix similarity(const HepSymMatrix &m1) const;
00203
00204
00205 HepSymMatrix similarityT(const HepMatrix &m1) const;
00206
00207
00208
00209 double similarity(const HepVector &v) const;
00210
00211
00212 HepSymMatrix sub(int min_row, int max_row) const;
00213
00214 void sub(int row, const HepSymMatrix &m1);
00215
00216 #ifdef HEP_CC_NEED_SUB_WITHOUT_CONST
00217 HepSymMatrix sub(int min_row, int max_row);
00218
00219
00220 #endif
00221
00222 inline HepSymMatrix inverse(int &ifail) const;
00223
00224
00225
00226 void invert(int &ifail);
00227
00228
00229
00230
00231
00232 double determinant() const;
00233
00234
00235 double trace() const;
00236
00237
00238 class HepSymMatrix_row {
00239 public:
00240 inline HepSymMatrix_row(HepSymMatrix&,int);
00241 inline double & operator[](int);
00242 private:
00243 HepSymMatrix& _a;
00244 int _r;
00245 };
00246 class HepSymMatrix_row_const {
00247 public:
00248 inline HepSymMatrix_row_const(const HepSymMatrix&,int);
00249 inline const double & operator[](int) const;
00250 private:
00251 const HepSymMatrix& _a;
00252 int _r;
00253 };
00254
00255
00256 inline HepSymMatrix_row operator[] (int);
00257 inline HepSymMatrix_row_const operator[] (int) const;
00258
00259
00260
00261
00262
00263
00264
00265
00266 void invertCholesky5 (int &ifail);
00267 void invertCholesky6 (int &ifail);
00268
00269
00270
00271 void invertHaywood4 (int & ifail);
00272 void invertHaywood5 (int &ifail);
00273 void invertHaywood6 (int &ifail);
00274 void invertBunchKaufman (int &ifail);
00275
00276 protected:
00277 inline int num_size() const;
00278
00279 private:
00280 friend class HepSymMatrix_row;
00281 friend class HepSymMatrix_row_const;
00282 friend class HepMatrix;
00283 friend class HepDiagMatrix;
00284
00285 friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
00286 friend double condition(const HepSymMatrix &m);
00287 friend void diag_step(HepSymMatrix *t,int begin,int end);
00288 friend void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end);
00289 friend HepMatrix diagonalize(HepSymMatrix *s);
00290 friend HepVector house(const HepSymMatrix &a,int row,int col);
00291 friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col);
00292
00293 friend HepSymMatrix operator+(const HepSymMatrix &m1,
00294 const HepSymMatrix &m2);
00295 friend HepSymMatrix operator-(const HepSymMatrix &m1,
00296 const HepSymMatrix &m2);
00297 friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00298 friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00299 friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00300 friend HepVector operator*(const HepSymMatrix &m1, const HepVector &m2);
00301
00302
00303 friend HepSymMatrix vT_times_v(const HepVector &v);
00304
00305
00306 std::vector<double,Alloc<double,25> > m;
00307 int nrow;
00308 int size;
00309
00310 static double posDefFraction5x5;
00311 static double adjustment5x5;
00312 static const double CHOLESKY_THRESHOLD_5x5;
00313 static const double CHOLESKY_CREEP_5x5;
00314
00315 static double posDefFraction6x6;
00316 static double adjustment6x6;
00317 static const double CHOLESKY_THRESHOLD_6x6;
00318 static const double CHOLESKY_CREEP_6x6;
00319
00320 void invert4 (int & ifail);
00321 void invert5 (int & ifail);
00322 void invert6 (int & ifail);
00323
00324 };
00325
00326
00327
00328
00329
00330
00331 HepStd::ostream& operator<<(HepStd::ostream &s, const HepSymMatrix &q);
00332
00333
00334 HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00335 HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00336 HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00337 HepSymMatrix operator*(double t, const HepSymMatrix &s1);
00338 HepSymMatrix operator*(const HepSymMatrix &s1, double t);
00339
00340
00341
00342 HepSymMatrix operator/(const HepSymMatrix &m1, double t);
00343
00344
00345 HepMatrix operator+(const HepMatrix &m1, const HepSymMatrix &s2);
00346 HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &m2);
00347 HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2);
00348
00349
00350 HepMatrix operator-(const HepMatrix &m1, const HepSymMatrix &s2);
00351 HepMatrix operator-(const HepSymMatrix &m1, const HepMatrix &m2);
00352 HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2);
00353
00354
00355 HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2);
00356
00357
00358 double condition(const HepSymMatrix &m);
00359
00360
00361 void diag_step(HepSymMatrix *t, int begin, int end);
00362 void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end);
00363
00364
00365 HepMatrix diagonalize(HepSymMatrix *s);
00366
00367
00368
00369 HepVector house(const HepSymMatrix &a, int row=1, int col=1);
00370 void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1);
00371
00372
00373 void tridiagonal(HepSymMatrix *a, HepMatrix *hsm);
00374 HepMatrix tridiagonal(HepSymMatrix *a);
00375
00376
00377 #ifdef HEP_NO_INLINE_IN_DECLARATION
00378 #undef inline
00379 #endif
00380
00381 #if defined(HEP_SHORT_NAMES)
00382 typedef HepSymMatrix SymMatrix;
00383 #endif
00384
00385 #ifndef HEP_DEBUG_INLINE
00386 #include "CLHEP/Matrix/SymMatrix.icc"
00387 #endif
00388
00389 #endif