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 #ifndef G4ErrorSymMatrix_hh
00037 #define G4ErrorSymMatrix_hh
00038
00039 #include <vector>
00040
00041 class G4ErrorMatrix;
00042
00043 class G4ErrorSymMatrix
00044 {
00045 public:
00046
00047 inline G4ErrorSymMatrix();
00048
00049
00050
00051 explicit G4ErrorSymMatrix(G4int p);
00052 G4ErrorSymMatrix(G4int p, G4int);
00053
00054
00055
00056
00057 G4ErrorSymMatrix(const G4ErrorSymMatrix &m1);
00058
00059
00060
00061
00062 virtual ~G4ErrorSymMatrix();
00063
00064
00065 inline G4int num_row() const;
00066 inline G4int num_col() const;
00067
00068
00069 const G4double & operator()(G4int row, G4int col) const;
00070 G4double & operator()(G4int row, G4int col);
00071
00072
00073
00074 const G4double & fast(G4int row, G4int col) const;
00075 G4double & fast(G4int row, G4int col);
00076
00077
00078
00079
00080 void assign(const G4ErrorMatrix &m2);
00081
00082
00083 void assign(const G4ErrorSymMatrix &m2);
00084
00085
00086 G4ErrorSymMatrix & operator*=(G4double t);
00087
00088
00089 G4ErrorSymMatrix & operator/=(G4double t);
00090
00091
00092 G4ErrorSymMatrix & operator+=( const G4ErrorSymMatrix &m2);
00093 G4ErrorSymMatrix & operator-=( const G4ErrorSymMatrix &m2);
00094
00095
00096 G4ErrorSymMatrix & operator=( const G4ErrorSymMatrix &m2);
00097
00098
00099 G4ErrorSymMatrix operator- () const;
00100
00101
00102 G4ErrorSymMatrix T() const;
00103
00104
00105 G4ErrorSymMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
00106
00107
00108 G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const;
00109 G4ErrorSymMatrix similarity(const G4ErrorSymMatrix &m1) const;
00110
00111
00112 G4ErrorSymMatrix similarityT(const G4ErrorMatrix &m1) const;
00113
00114
00115
00116 G4ErrorSymMatrix sub(G4int min_row, G4int max_row) const;
00117
00118
00119 void sub(G4int row, const G4ErrorSymMatrix &m1);
00120
00121
00122 G4ErrorSymMatrix sub(G4int min_row, G4int max_row);
00123
00124
00125
00126 inline G4ErrorSymMatrix inverse(G4int &ifail) const;
00127
00128
00129
00130 void invert(G4int &ifail);
00131
00132
00133
00134
00135
00136 G4double determinant() const;
00137
00138
00139 G4double trace() const;
00140
00141
00142 class G4ErrorSymMatrix_row
00143 {
00144 public:
00145 inline G4ErrorSymMatrix_row(G4ErrorSymMatrix&,G4int);
00146 inline G4double & operator[](G4int);
00147 private:
00148 G4ErrorSymMatrix& _a;
00149 G4int _r;
00150 };
00151 class G4ErrorSymMatrix_row_const
00152 {
00153 public:
00154 inline G4ErrorSymMatrix_row_const(const G4ErrorSymMatrix&,G4int);
00155 inline const G4double & operator[](G4int) const;
00156 private:
00157 const G4ErrorSymMatrix& _a;
00158 G4int _r;
00159 };
00160
00161
00162 inline G4ErrorSymMatrix_row operator[] (G4int);
00163 inline G4ErrorSymMatrix_row_const operator[] (G4int) const;
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 void invertCholesky5 (G4int &ifail);
00174 void invertCholesky6 (G4int &ifail);
00175
00176
00177
00178
00179 void invertHaywood4 (G4int & ifail);
00180 void invertHaywood5 (G4int &ifail);
00181 void invertHaywood6 (G4int &ifail);
00182 void invertBunchKaufman (G4int &ifail);
00183
00184 protected:
00185
00186 inline G4int num_size() const;
00187
00188 private:
00189
00190 friend class G4ErrorSymMatrix_row;
00191 friend class G4ErrorSymMatrix_row_const;
00192 friend class G4ErrorMatrix;
00193
00194 friend void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
00195 friend G4double condition(const G4ErrorSymMatrix &m);
00196 friend void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
00197 friend void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u,
00198 G4int begin, G4int end);
00199 friend G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
00200 friend void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
00201 G4int row, G4int col);
00202
00203 friend G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1,
00204 const G4ErrorSymMatrix &m2);
00205 friend G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1,
00206 const G4ErrorSymMatrix &m2);
00207 friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00208 const G4ErrorSymMatrix &m2);
00209 friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00210 const G4ErrorMatrix &m2);
00211 friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
00212 const G4ErrorSymMatrix &m2);
00213
00214
00215
00216
00217 std::vector<G4double > m;
00218
00219 G4int nrow;
00220 G4int size;
00221
00222 static G4double posDefFraction5x5;
00223 static G4double adjustment5x5;
00224 static const G4double CHOLESKY_THRESHOLD_5x5;
00225 static const G4double CHOLESKY_CREEP_5x5;
00226
00227 static G4double posDefFraction6x6;
00228 static G4double adjustment6x6;
00229 static const G4double CHOLESKY_THRESHOLD_6x6;
00230 static const G4double CHOLESKY_CREEP_6x6;
00231
00232 void invert4 (G4int & ifail);
00233 void invert5 (G4int & ifail);
00234 void invert6 (G4int & ifail);
00235 };
00236
00237
00238
00239
00240
00241
00242 std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q);
00243
00244
00245 G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
00246 const G4ErrorSymMatrix &m2);
00247 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00248 const G4ErrorMatrix &m2);
00249 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00250 const G4ErrorSymMatrix &m2);
00251 G4ErrorSymMatrix operator*(G4double t, const G4ErrorSymMatrix &s1);
00252 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &s1, G4double t);
00253
00254
00255
00256 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1, G4double t);
00257
00258
00259 G4ErrorMatrix operator+(const G4ErrorMatrix &m1,
00260 const G4ErrorSymMatrix &s2);
00261 G4ErrorMatrix operator+(const G4ErrorSymMatrix &s1,
00262 const G4ErrorMatrix &m2);
00263 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &s1,
00264 const G4ErrorSymMatrix &s2);
00265
00266
00267 G4ErrorMatrix operator-(const G4ErrorMatrix &m1,
00268 const G4ErrorSymMatrix &s2);
00269 G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1,
00270 const G4ErrorMatrix &m2);
00271 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &s1,
00272 const G4ErrorSymMatrix &s2);
00273
00274
00275 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &s1,
00276 const G4ErrorSymMatrix &s2);
00277
00278
00279 G4double condition(const G4ErrorSymMatrix &m);
00280
00281
00282 void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
00283 void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u, G4int begin, G4int end);
00284
00285
00286 G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
00287
00288
00289
00290 void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
00291 G4int row=1, G4int col=1);
00292
00293
00294 void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
00295 G4ErrorMatrix tridiagonal(G4ErrorSymMatrix *a);
00296
00297
00298 #include "G4ErrorSymMatrix.icc"
00299
00300 #endif