| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4ErrorSymMatrix.cc,v 1.3 2007/06/21 15:04:10 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02 $
26 // 28 //
27 // ------------------------------------------- 29 // ------------------------------------------------------------
28 // GEANT 4 class implementation file 30 // GEANT 4 class implementation file
29 // ------------------------------------------- 31 // ------------------------------------------------------------
30 32
31 #include "globals.hh" 33 #include "globals.hh"
32 #include <iostream> 34 #include <iostream>
33 #include <cmath> 35 #include <cmath>
34 36
35 #include "G4ErrorSymMatrix.hh" 37 #include "G4ErrorSymMatrix.hh"
36 #include "G4ErrorMatrix.hh" 38 #include "G4ErrorMatrix.hh"
37 39
38 // Simple operation for all elements 40 // Simple operation for all elements
39 41
40 #define SIMPLE_UOP(OPER) << 42 #define SIMPLE_UOP(OPER) \
41 G4ErrorMatrixIter a = m.begin(); << 43 G4ErrorMatrixIter a=m.begin(); \
42 G4ErrorMatrixIter e = m.begin() + num_size() << 44 G4ErrorMatrixIter e=m.begin()+num_size(); \
43 for(; a < e; a++) << 45 for(;a<e; a++) (*a) OPER t;
44 (*a) OPER t; << 46
45 << 47 #define SIMPLE_BOP(OPER) \
46 #define SIMPLE_BOP(OPER) << 48 G4ErrorMatrixIter a=m.begin(); \
47 G4ErrorMatrixIter a = m.begin(); << 49 G4ErrorMatrixConstIter b=m2.m.begin(); \
48 G4ErrorMatrixConstIter b = mat2.m.begin(); << 50 G4ErrorMatrixConstIter e=m.begin()+num_size(); \
49 G4ErrorMatrixConstIter e = m.begin() + num_s << 51 for(;a<e; a++, b++) (*a) OPER (*b);
50 for(; a < e; a++, b++) << 52
51 (*a) OPER(*b); << 53 #define SIMPLE_TOP(OPER) \
52 << 54 G4ErrorMatrixConstIter a=m1.m.begin(); \
53 #define SIMPLE_TOP(OPER) << 55 G4ErrorMatrixConstIter b=m2.m.begin(); \
54 G4ErrorMatrixConstIter a = mat1.m.begin(); << 56 G4ErrorMatrixIter t=mret.m.begin(); \
55 G4ErrorMatrixConstIter b = mat2.m.begin(); << 57 G4ErrorMatrixConstIter e=m1.m.begin()+m1.num_size(); \
56 G4ErrorMatrixIter t = mret.m.begin(); << 58 for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
57 G4ErrorMatrixConstIter e = mat1.m.begin() + << 59
58 for(; a < e; a++, b++, t++) << 60 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
59 (*t) = (*a) OPER(*b); << 61 if (r1!=r2 || c1!=c2) { \
60 << 62 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
61 #define CHK_DIM_2(r1, r2, c1, c2, fun) << 63 }
62 if(r1 != r2 || c1 != c2) << 64
63 { << 65 #define CHK_DIM_1(c1,r2,fun) \
64 G4ErrorMatrix::error("Range error in Matri << 66 if (c1!=r2) { \
65 } << 67 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
66 << 68 }
67 #define CHK_DIM_1(c1, r2, fun) <<
68 if(c1 != r2) <<
69 { <<
70 G4ErrorMatrix::error("Range error in Matri <<
71 } <<
72 69
73 // Constructors. (Default constructors are inl 70 // Constructors. (Default constructors are inlined and in .icc file)
74 71
75 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p) 72 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
76 : m(p * (p + 1) / 2) << 73 : m(p*(p+1)/2), nrow(p)
77 , nrow(p) <<
78 { 74 {
79 size = nrow * (nrow + 1) / 2; << 75 size = nrow * (nrow+1) / 2;
80 m.assign(size, 0); << 76 m.assign(size,0);
81 } 77 }
82 78
83 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4 79 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4int init)
84 : m(p * (p + 1) / 2) << 80 : m(p*(p+1)/2), nrow(p)
85 , nrow(p) <<
86 { 81 {
87 size = nrow * (nrow + 1) / 2; << 82 size = nrow * (nrow+1) / 2;
88 83
89 m.assign(size, 0); << 84 m.assign(size,0);
90 switch(init) << 85 switch(init)
91 { << 86 {
92 case 0: << 87 case 0:
93 break; 88 break;
94 << 89
95 case 1: << 90 case 1:
96 { <<
97 G4ErrorMatrixIter a = m.begin(); <<
98 for(G4int i = 1; i <= nrow; i++) <<
99 { 91 {
100 *a = 1.0; << 92 G4ErrorMatrixIter a = m.begin();
101 a += (i + 1); << 93 for(G4int i=1;i<=nrow;i++)
>> 94 {
>> 95 *a = 1.0;
>> 96 a += (i+1);
>> 97 }
>> 98 break;
102 } 99 }
103 break; << 100 default:
104 } <<
105 default: <<
106 G4ErrorMatrix::error("G4ErrorSymMatrix: 101 G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
107 } << 102 }
108 } 103 }
109 104
110 // 105 //
111 // Destructor 106 // Destructor
112 // 107 //
113 108
114 G4ErrorSymMatrix::~G4ErrorSymMatrix() {} << 109 G4ErrorSymMatrix::~G4ErrorSymMatrix()
>> 110 {
>> 111 }
115 112
116 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4Err << 113 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4ErrorSymMatrix &m1)
117 : m(mat1.size) << 114 : m(m1.size), nrow(m1.nrow), size(m1.size)
118 , nrow(mat1.nrow) <<
119 , size(mat1.size) <<
120 { 115 {
121 m = mat1.m; << 116 m = m1.m;
122 } 117 }
123 118
124 // 119 //
125 // 120 //
126 // Sub matrix 121 // Sub matrix
127 // 122 //
128 // 123 //
129 124
130 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m 125 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) const
131 { 126 {
132 G4ErrorSymMatrix mret(max_row - min_row + 1) << 127 G4ErrorSymMatrix mret(max_row-min_row+1);
133 if(max_row > num_row()) 128 if(max_row > num_row())
134 { << 129 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
135 G4ErrorMatrix::error("G4ErrorSymMatrix::su << 130 G4ErrorMatrixIter a = mret.m.begin();
136 } << 131 G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
137 G4ErrorMatrixIter a = mret.m.begin(); << 132 for(G4int irow=1; irow<=mret.num_row(); irow++)
138 G4ErrorMatrixConstIter b1 = m.begin() + (min <<
139 for(G4int irow = 1; irow <= mret.num_row(); <<
140 { 133 {
141 G4ErrorMatrixConstIter b = b1; 134 G4ErrorMatrixConstIter b = b1;
142 for(G4int icol = 1; icol <= irow; icol++) << 135 for(G4int icol=1; icol<=irow; icol++)
143 { 136 {
144 *(a++) = *(b++); 137 *(a++) = *(b++);
145 } 138 }
146 b1 += irow + min_row - 1; << 139 b1 += irow+min_row-1;
147 } 140 }
148 return mret; 141 return mret;
149 } 142 }
150 143
151 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m << 144 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row)
152 { 145 {
153 G4ErrorSymMatrix mret(max_row - min_row + 1) << 146 G4ErrorSymMatrix mret(max_row-min_row+1);
154 if(max_row > num_row()) 147 if(max_row > num_row())
155 { << 148 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
156 G4ErrorMatrix::error("G4ErrorSymMatrix::su << 149 G4ErrorMatrixIter a = mret.m.begin();
157 } << 150 G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
158 G4ErrorMatrixIter a = mret.m.begin(); << 151 for(G4int irow=1; irow<=mret.num_row(); irow++)
159 G4ErrorMatrixIter b1 = m.begin() + (min_row <<
160 for(G4int irow = 1; irow <= mret.num_row(); <<
161 { 152 {
162 G4ErrorMatrixIter b = b1; 153 G4ErrorMatrixIter b = b1;
163 for(G4int icol = 1; icol <= irow; icol++) << 154 for(G4int icol=1; icol<=irow; icol++)
164 { 155 {
165 *(a++) = *(b++); 156 *(a++) = *(b++);
166 } 157 }
167 b1 += irow + min_row - 1; << 158 b1 += irow+min_row-1;
168 } 159 }
169 return mret; 160 return mret;
170 } 161 }
171 162
172 void G4ErrorSymMatrix::sub(G4int row, const G4 << 163 void G4ErrorSymMatrix::sub(G4int row,const G4ErrorSymMatrix &m1)
173 { 164 {
174 if(row < 1 || row + mat1.num_row() - 1 > num << 165 if(row <1 || row+m1.num_row()-1 > num_row() )
175 { << 166 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
176 G4ErrorMatrix::error("G4ErrorSymMatrix::su << 167 G4ErrorMatrixConstIter a = m1.m.begin();
177 } << 168 G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
178 G4ErrorMatrixConstIter a = mat1.m.begin(); << 169 for(G4int irow=1; irow<=m1.num_row(); irow++)
179 G4ErrorMatrixIter b1 = m.begin() + (row <<
180 for(G4int irow = 1; irow <= mat1.num_row(); <<
181 { 170 {
182 G4ErrorMatrixIter b = b1; 171 G4ErrorMatrixIter b = b1;
183 for(G4int icol = 1; icol <= irow; icol++) << 172 for(G4int icol=1; icol<=irow; icol++)
184 { 173 {
185 *(b++) = *(a++); 174 *(b++) = *(a++);
186 } 175 }
187 b1 += irow + row - 1; << 176 b1 += irow+row-1;
188 } 177 }
189 } 178 }
190 179
191 // 180 //
192 // Direct sum of two matricies 181 // Direct sum of two matricies
193 // 182 //
194 183
195 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix& << 184 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &m1,
196 const G4ErrorSymMatrix& << 185 const G4ErrorSymMatrix &m2)
197 { 186 {
198 G4ErrorSymMatrix mret(mat1.num_row() + mat2. << 187 G4ErrorSymMatrix mret(m1.num_row() + m2.num_row(), 0);
199 mret.sub(1, mat1); << 188 mret.sub(1,m1);
200 mret.sub(mat1.num_row() + 1, mat2); << 189 mret.sub(m1.num_row()+1,m2);
201 return mret; 190 return mret;
202 } 191 }
203 192
204 /* ------------------------------------------- 193 /* -----------------------------------------------------------------------
205 This section contains support routines for 194 This section contains support routines for matrix.h. This section contains
206 The two argument functions +,-. They call t 195 The two argument functions +,-. They call the copy constructor and +=,-=.
207 ------------------------------------------- 196 ----------------------------------------------------------------------- */
208 197
209 G4ErrorSymMatrix G4ErrorSymMatrix::operator-() << 198 G4ErrorSymMatrix G4ErrorSymMatrix::operator- () const
210 { 199 {
211 G4ErrorSymMatrix mat2(nrow); << 200 G4ErrorSymMatrix m2(nrow);
212 G4ErrorMatrixConstIter a = m.begin(); << 201 G4ErrorMatrixConstIter a=m.begin();
213 G4ErrorMatrixIter b = mat2.m.begin(); << 202 G4ErrorMatrixIter b=m2.m.begin();
214 G4ErrorMatrixConstIter e = m.begin() + num_s << 203 G4ErrorMatrixConstIter e=m.begin()+num_size();
215 for(; a < e; a++, b++) << 204 for(;a<e; a++, b++) { (*b) = -(*a); }
216 { << 205 return m2;
217 (*b) = -(*a); <<
218 } <<
219 return mat2; <<
220 } 206 }
221 207
222 G4ErrorMatrix operator+(const G4ErrorMatrix& m << 208 G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
223 { 209 {
224 G4ErrorMatrix mret(mat1); << 210 G4ErrorMatrix mret(m1);
225 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 211 CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+);
226 mret += mat2; << 212 mret += m2;
227 return mret; 213 return mret;
228 } 214 }
229 215
230 G4ErrorMatrix operator+(const G4ErrorSymMatrix << 216 G4ErrorMatrix operator+(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
231 { 217 {
232 G4ErrorMatrix mret(mat2); << 218 G4ErrorMatrix mret(m2);
233 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 219 CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),+);
234 mret += mat1; << 220 mret += m1;
235 return mret; 221 return mret;
236 } 222 }
237 223
238 G4ErrorSymMatrix operator+(const G4ErrorSymMat << 224 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1,
239 const G4ErrorSymMat << 225 const G4ErrorSymMatrix &m2)
240 { 226 {
241 G4ErrorSymMatrix mret(mat1.nrow); << 227 G4ErrorSymMatrix mret(m1.nrow);
242 CHK_DIM_1(mat1.nrow, mat2.nrow, +); << 228 CHK_DIM_1(m1.nrow, m2.nrow,+);
243 SIMPLE_TOP(+) 229 SIMPLE_TOP(+)
244 return mret; 230 return mret;
245 } 231 }
246 232
247 // 233 //
248 // operator - 234 // operator -
249 // 235 //
250 236
251 G4ErrorMatrix operator-(const G4ErrorMatrix& m << 237 G4ErrorMatrix operator-(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
252 { 238 {
253 G4ErrorMatrix mret(mat1); << 239 G4ErrorMatrix mret(m1);
254 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 240 CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
255 mret -= mat2; << 241 mret -= m2;
256 return mret; 242 return mret;
257 } 243 }
258 244
259 G4ErrorMatrix operator-(const G4ErrorSymMatrix << 245 G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
260 { 246 {
261 G4ErrorMatrix mret(mat1); << 247 G4ErrorMatrix mret(m1);
262 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 248 CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
263 mret -= mat2; << 249 mret -= m2;
264 return mret; 250 return mret;
265 } 251 }
266 252
267 G4ErrorSymMatrix operator-(const G4ErrorSymMat << 253 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1,
268 const G4ErrorSymMat << 254 const G4ErrorSymMatrix &m2)
269 { 255 {
270 G4ErrorSymMatrix mret(mat1.num_row()); << 256 G4ErrorSymMatrix mret(m1.num_row());
271 CHK_DIM_1(mat1.num_row(), mat2.num_row(), -) << 257 CHK_DIM_1(m1.num_row(),m2.num_row(),-);
272 SIMPLE_TOP(-) 258 SIMPLE_TOP(-)
273 return mret; 259 return mret;
274 } 260 }
275 261
276 /* ------------------------------------------- 262 /* -----------------------------------------------------------------------
277 This section contains support routines for 263 This section contains support routines for matrix.h. This file contains
278 The two argument functions *,/. They call c 264 The two argument functions *,/. They call copy constructor and then /=,*=.
279 ------------------------------------------- 265 ----------------------------------------------------------------------- */
280 266
281 G4ErrorSymMatrix operator/(const G4ErrorSymMat << 267 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1,G4double t)
282 { 268 {
283 G4ErrorSymMatrix mret(mat1); << 269 G4ErrorSymMatrix mret(m1);
284 mret /= t; 270 mret /= t;
285 return mret; 271 return mret;
286 } 272 }
287 273
288 G4ErrorSymMatrix operator*(const G4ErrorSymMat << 274 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &m1,G4double t)
289 { 275 {
290 G4ErrorSymMatrix mret(mat1); << 276 G4ErrorSymMatrix mret(m1);
291 mret *= t; 277 mret *= t;
292 return mret; 278 return mret;
293 } 279 }
294 280
295 G4ErrorSymMatrix operator*(G4double t, const G << 281 G4ErrorSymMatrix operator*(G4double t,const G4ErrorSymMatrix &m1)
296 { 282 {
297 G4ErrorSymMatrix mret(mat1); << 283 G4ErrorSymMatrix mret(m1);
298 mret *= t; 284 mret *= t;
299 return mret; 285 return mret;
300 } 286 }
301 287
302 G4ErrorMatrix operator*(const G4ErrorMatrix& m << 288 G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
303 { 289 {
304 G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 290 G4ErrorMatrix mret(m1.num_row(),m2.num_col());
305 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 291 CHK_DIM_1(m1.num_col(),m2.num_row(),*);
306 G4ErrorMatrixConstIter mit1, mit2, sp, snp; << 292 G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
307 G4double temp; 293 G4double temp;
308 G4ErrorMatrixIter mir = mret.m.begin(); << 294 G4ErrorMatrixIter mir=mret.m.begin();
309 for(mit1 = mat1.m.begin(); << 295 G4int step,stept;
310 mit1 < mat1.m.begin() + mat1.num_row() * << 296 for(mit1=m1.m.begin();
311 { << 297 mit1<m1.m.begin()+m1.num_row()*m1.num_col(); mit1 = mit2)
312 snp = mat2.m.begin(); << 298 {
313 for(int step = 1; step <= mat2.num_row(); << 299 for(step=1,snp=m2.m.begin();step<=m2.num_row();)
314 { << 300 {
315 mit2 = mit1; << 301 mit2=mit1;
316 sp = snp; << 302 sp=snp;
317 snp += step; << 303 snp+=step;
318 temp = 0; << 304 temp=0;
319 while(sp < snp) // Loop checking, 06.08 << 305 while(sp<snp)
320 { << 306 {
321 temp += *(sp++) * (*(mit2++)); << 307 temp+=*(sp++)*(*(mit2++));
322 } << 308 }
323 if(step < mat2.num_row()) << 309 sp+=step-1;
324 { // only if we aren't on the last row << 310 for(stept=++step;stept<=m2.num_row();stept++)
325 sp += step - 1; << 311 {
326 for(int stept = step + 1; stept <= mat << 312 temp+=*sp*(*(mit2++));
327 { << 313 sp+=stept;
328 temp += *sp * (*(mit2++)); << 314 }
329 if(stept < mat2.num_row()) << 315 *(mir++)=temp;
330 sp += stept; << 316 }
331 } << 317 }
332 } // if(step <<
333 *(mir++) = temp; <<
334 } // for(step <<
335 } // for(mit1 <<
336 return mret; 318 return mret;
337 } 319 }
338 320
339 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 321 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
340 { 322 {
341 G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 323 G4ErrorMatrix mret(m1.num_row(),m2.num_col());
342 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 324 CHK_DIM_1(m1.num_col(),m2.num_row(),*);
343 G4int step, stept; << 325 G4int step,stept;
344 G4ErrorMatrixConstIter mit1, mit2, sp, snp; << 326 G4ErrorMatrixConstIter mit1,mit2,sp,snp;
345 G4double temp; 327 G4double temp;
346 G4ErrorMatrixIter mir = mret.m.begin(); << 328 G4ErrorMatrixIter mir=mret.m.begin();
347 for(step = 1, snp = mat1.m.begin(); step <= << 329 for(step=1,snp=m1.m.begin();step<=m1.num_row();snp+=step++)
348 { 330 {
349 for(mit1 = mat2.m.begin(); mit1 < mat2.m.b << 331 for(mit1=m2.m.begin();mit1<m2.m.begin()+m2.num_col();mit1++)
350 { 332 {
351 mit2 = mit1; << 333 mit2=mit1;
352 sp = snp; << 334 sp=snp;
353 temp = 0; << 335 temp=0;
354 while(sp < snp + step) // Loop checking << 336 while(sp<snp+step)
355 { 337 {
356 temp += *mit2 * (*(sp++)); << 338 temp+=*mit2*(*(sp++));
357 mit2 += mat2.num_col(); << 339 mit2+=m2.num_col();
358 } 340 }
359 sp += step - 1; << 341 sp+=step-1;
360 for(stept = step + 1; stept <= mat1.num_ << 342 for(stept=step+1;stept<=m1.num_row();stept++)
361 { 343 {
362 temp += *mit2 * (*sp); << 344 temp+=*mit2*(*sp);
363 mit2 += mat2.num_col(); << 345 mit2+=m2.num_col();
364 sp += stept; << 346 sp+=stept;
365 } 347 }
366 *(mir++) = temp; << 348 *(mir++)=temp;
367 } 349 }
368 } 350 }
369 return mret; 351 return mret;
370 } 352 }
371 353
372 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 354 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
373 const G4ErrorSymMatrix <<
374 { 355 {
375 G4ErrorMatrix mret(mat1.num_row(), mat1.num_ << 356 G4ErrorMatrix mret(m1.num_row(),m1.num_row());
376 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 357 CHK_DIM_1(m1.num_col(),m2.num_row(),*);
377 G4int step1, stept1, step2, stept2; << 358 G4int step1,stept1,step2,stept2;
378 G4ErrorMatrixConstIter snp1, sp1, snp2, sp2; << 359 G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
379 G4double temp; 360 G4double temp;
380 G4ErrorMatrixIter mr = mret.m.begin(); 361 G4ErrorMatrixIter mr = mret.m.begin();
381 for(step1 = 1, snp1 = mat1.m.begin(); step1 << 362 for(step1=1,snp1=m1.m.begin();step1<=m1.num_row();snp1+=step1++)
382 snp1 += step1++) <<
383 { 363 {
384 for(step2 = 1, snp2 = mat2.m.begin(); step << 364 for(step2=1,snp2=m2.m.begin();step2<=m2.num_row();)
385 { 365 {
386 sp1 = snp1; << 366 sp1=snp1;
387 sp2 = snp2; << 367 sp2=snp2;
388 snp2 += step2; << 368 snp2+=step2;
389 temp = 0; << 369 temp=0;
390 if(step1 < step2) << 370 if(step1<step2)
391 { << 371 {
392 while(sp1 < snp1 + step1) // Loop che << 372 while(sp1<snp1+step1)
393 { << 373 { temp+=(*(sp1++))*(*(sp2++)); }
394 temp += (*(sp1++)) * (*(sp2++)); << 374 sp1+=step1-1;
395 } << 375 for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
396 sp1 += step1 - 1; << 376 { temp+=(*sp1)*(*(sp2++)); }
397 for(stept1 = step1 + 1; stept1 != step << 377 sp2+=step2-1;
398 { << 378 for(stept2=++step2;stept2<=m2.num_row();sp1+=stept1++,sp2+=stept2++)
399 temp += (*sp1) * (*(sp2++)); << 379 { temp+=(*sp1)*(*sp2); }
400 } <<
401 sp2 += step2 - 1; <<
402 for(stept2 = ++step2; stept2 <= mat2.n <<
403 sp1 += stept1++, sp2 += stept2++) <<
404 { <<
405 temp += (*sp1) * (*sp2); <<
406 } <<
407 } 380 }
408 else 381 else
409 { 382 {
410 while(sp2 < snp2) // Loop checking, 0 << 383 while(sp2<snp2)
411 { << 384 { temp+=(*(sp1++))*(*(sp2++)); }
412 temp += (*(sp1++)) * (*(sp2++)); << 385 sp2+=step2-1;
413 } << 386 for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
414 sp2 += step2 - 1; << 387 { temp+=(*(sp1++))*(*sp2); }
415 for(stept2 = ++step2; stept2 != step1 << 388 sp1+=step1-1;
416 { << 389 for(stept1=step1+1;stept1<=m1.num_row();sp1+=stept1++,sp2+=stept2++)
417 temp += (*(sp1++)) * (*sp2); << 390 { temp+=(*sp1)*(*sp2); }
418 } <<
419 sp1 += step1 - 1; <<
420 for(stept1 = step1 + 1; stept1 <= mat1 <<
421 sp1 += stept1++, sp2 += stept2++) <<
422 { <<
423 temp += (*sp1) * (*sp2); <<
424 } <<
425 } 391 }
426 *(mr++) = temp; << 392 *(mr++)=temp;
427 } 393 }
428 } 394 }
429 return mret; 395 return mret;
430 } 396 }
431 397
432 /* ------------------------------------------- 398 /* -----------------------------------------------------------------------
433 This section contains the assignment and in 399 This section contains the assignment and inplace operators =,+=,-=,*=,/=.
434 ------------------------------------------- 400 ----------------------------------------------------------------------- */
435 401
436 G4ErrorMatrix& G4ErrorMatrix::operator+=(const << 402 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorSymMatrix &m2)
437 { 403 {
438 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 404 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
439 G4int n = num_col(); << 405 G4int n = num_col();
440 G4ErrorMatrixConstIter sjk = mat2.m.begin(); << 406 G4ErrorMatrixConstIter sjk = m2.m.begin();
441 G4ErrorMatrixIter m1j = m.begin(); << 407 G4ErrorMatrixIter m1j = m.begin();
442 G4ErrorMatrixIter mj = m.begin(); << 408 G4ErrorMatrixIter mj = m.begin();
443 // j >= k 409 // j >= k
444 for(G4int j = 1; j <= num_row(); j++) << 410 for(G4int j=1;j<=num_row();j++)
445 { 411 {
446 G4ErrorMatrixIter mjk = mj; 412 G4ErrorMatrixIter mjk = mj;
447 G4ErrorMatrixIter mkj = m1j; 413 G4ErrorMatrixIter mkj = m1j;
448 for(G4int k = 1; k <= j; k++) << 414 for(G4int k=1;k<=j;k++)
449 { 415 {
450 *(mjk++) += *sjk; 416 *(mjk++) += *sjk;
451 if(j != k) << 417 if(j!=k) *mkj += *sjk;
452 *mkj += *sjk; <<
453 sjk++; 418 sjk++;
454 mkj += n; 419 mkj += n;
455 } 420 }
456 mj += n; 421 mj += n;
457 m1j++; 422 m1j++;
458 } 423 }
459 return (*this); 424 return (*this);
460 } 425 }
461 426
462 G4ErrorSymMatrix& G4ErrorSymMatrix::operator+= << 427 G4ErrorSymMatrix & G4ErrorSymMatrix::operator+=(const G4ErrorSymMatrix &m2)
463 { 428 {
464 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 429 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
465 SIMPLE_BOP(+=) 430 SIMPLE_BOP(+=)
466 return (*this); 431 return (*this);
467 } 432 }
468 433
469 G4ErrorMatrix& G4ErrorMatrix::operator-=(const << 434 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorSymMatrix &m2)
470 { 435 {
471 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 436 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
472 G4int n = num_col(); << 437 G4int n = num_col();
473 G4ErrorMatrixConstIter sjk = mat2.m.begin(); << 438 G4ErrorMatrixConstIter sjk = m2.m.begin();
474 G4ErrorMatrixIter m1j = m.begin(); << 439 G4ErrorMatrixIter m1j = m.begin();
475 G4ErrorMatrixIter mj = m.begin(); << 440 G4ErrorMatrixIter mj = m.begin();
476 // j >= k 441 // j >= k
477 for(G4int j = 1; j <= num_row(); j++) << 442 for(G4int j=1;j<=num_row();j++)
478 { 443 {
479 G4ErrorMatrixIter mjk = mj; 444 G4ErrorMatrixIter mjk = mj;
480 G4ErrorMatrixIter mkj = m1j; 445 G4ErrorMatrixIter mkj = m1j;
481 for(G4int k = 1; k <= j; k++) << 446 for(G4int k=1;k<=j;k++)
482 { 447 {
483 *(mjk++) -= *sjk; 448 *(mjk++) -= *sjk;
484 if(j != k) << 449 if(j!=k) *mkj -= *sjk;
485 *mkj -= *sjk; <<
486 sjk++; 450 sjk++;
487 mkj += n; 451 mkj += n;
488 } 452 }
489 mj += n; 453 mj += n;
490 m1j++; 454 m1j++;
491 } 455 }
492 return (*this); 456 return (*this);
493 } 457 }
494 458
495 G4ErrorSymMatrix& G4ErrorSymMatrix::operator-= << 459 G4ErrorSymMatrix & G4ErrorSymMatrix::operator-=(const G4ErrorSymMatrix &m2)
496 { 460 {
497 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 461 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
498 SIMPLE_BOP(-=) 462 SIMPLE_BOP(-=)
499 return (*this); 463 return (*this);
500 } 464 }
501 465
502 G4ErrorSymMatrix& G4ErrorSymMatrix::operator/= << 466 G4ErrorSymMatrix & G4ErrorSymMatrix::operator/=(G4double t)
503 { 467 {
504 SIMPLE_UOP(/=) 468 SIMPLE_UOP(/=)
505 return (*this); 469 return (*this);
506 } 470 }
507 471
508 G4ErrorSymMatrix& G4ErrorSymMatrix::operator*= << 472 G4ErrorSymMatrix & G4ErrorSymMatrix::operator*=(G4double t)
509 { 473 {
510 SIMPLE_UOP(*=) 474 SIMPLE_UOP(*=)
511 return (*this); 475 return (*this);
512 } 476 }
513 477
514 G4ErrorMatrix& G4ErrorMatrix::operator=(const << 478 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorSymMatrix &m1)
515 { 479 {
516 if(mat1.nrow * mat1.nrow != size) << 480 if(m1.nrow*m1.nrow != size)
517 { << 481 {
518 size = mat1.nrow * mat1.nrow; << 482 size = m1.nrow * m1.nrow;
519 m.resize(size); << 483 m.resize(size);
520 } << 484 }
521 nrow = mat1.nrow; << 485 nrow = m1.nrow;
522 ncol = mat1.nrow; << 486 ncol = m1.nrow;
523 G4int n = ncol; << 487 G4int n = ncol;
524 G4ErrorMatrixConstIter sjk = mat1.m.begin(); << 488 G4ErrorMatrixConstIter sjk = m1.m.begin();
525 G4ErrorMatrixIter m1j = m.begin(); << 489 G4ErrorMatrixIter m1j = m.begin();
526 G4ErrorMatrixIter mj = m.begin(); << 490 G4ErrorMatrixIter mj = m.begin();
527 // j >= k << 491 // j >= k
528 for(G4int j = 1; j <= num_row(); j++) << 492 for(G4int j=1;j<=num_row();j++)
529 { << 493 {
530 G4ErrorMatrixIter mjk = mj; << 494 G4ErrorMatrixIter mjk = mj;
531 G4ErrorMatrixIter mkj = m1j; << 495 G4ErrorMatrixIter mkj = m1j;
532 for(G4int k = 1; k <= j; k++) << 496 for(G4int k=1;k<=j;k++)
533 { << 497 {
534 *(mjk++) = *sjk; << 498 *(mjk++) = *sjk;
535 if(j != k) << 499 if(j!=k) *mkj = *sjk;
536 *mkj = *sjk; << 500 sjk++;
537 sjk++; << 501 mkj += n;
538 mkj += n; << 502 }
539 } << 503 mj += n;
540 mj += n; << 504 m1j++;
541 m1j++; << 505 }
542 } << 506 return (*this);
543 return (*this); << 507 }
544 } << 508
545 << 509 G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &m1)
546 G4ErrorSymMatrix& G4ErrorSymMatrix::operator=( << 510 {
547 { << 511 if(m1.nrow != nrow)
548 if(&mat1 == this) << 512 {
549 { << 513 nrow = m1.nrow;
550 return *this; << 514 size = m1.size;
551 } << 515 m.resize(size);
552 if(mat1.nrow != nrow) << 516 }
553 { << 517 m = m1.m;
554 nrow = mat1.nrow; << 518 return (*this);
555 size = mat1.size; <<
556 m.resize(size); <<
557 } <<
558 m = mat1.m; <<
559 return (*this); <<
560 } 519 }
561 520
562 // Print the Matrix. 521 // Print the Matrix.
563 522
564 std::ostream& operator<<(std::ostream& os, con << 523 std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q)
565 { 524 {
566 os << G4endl; << 525 s << G4endl;
567 526
568 // Fixed format needs 3 extra characters for 527 // Fixed format needs 3 extra characters for field,
569 // while scientific needs 7 528 // while scientific needs 7
570 529
571 std::size_t width; << 530 G4int width;
572 if(os.flags() & std::ios::fixed) << 531 if(s.flags() & std::ios::fixed)
573 { 532 {
574 width = os.precision() + 3; << 533 width = s.precision()+3;
575 } 534 }
576 else 535 else
577 { 536 {
578 width = os.precision() + 7; << 537 width = s.precision()+7;
579 } 538 }
580 for(G4int irow = 1; irow <= q.num_row(); ++i << 539 for(G4int irow = 1; irow<= q.num_row(); irow++)
581 { 540 {
582 for(G4int icol = 1; icol <= q.num_col(); + << 541 for(G4int icol = 1; icol <= q.num_col(); icol++)
583 { 542 {
584 os.width(width); << 543 s.width(width);
585 os << q(irow, icol) << " "; << 544 s << q(irow,icol) << " ";
586 } 545 }
587 os << G4endl; << 546 s << G4endl;
588 } 547 }
589 return os; << 548 return s;
590 } 549 }
591 550
592 G4ErrorSymMatrix G4ErrorSymMatrix::apply(G4dou << 551 G4ErrorSymMatrix G4ErrorSymMatrix::
593 << 552 apply(G4double (*f)(G4double, G4int, G4int)) const
594 { 553 {
595 G4ErrorSymMatrix mret(num_row()); 554 G4ErrorSymMatrix mret(num_row());
596 G4ErrorMatrixConstIter a = m.cbegin(); << 555 G4ErrorMatrixConstIter a = m.begin();
597 G4ErrorMatrixIter b = mret.m.begin(); << 556 G4ErrorMatrixIter b = mret.m.begin();
598 for(G4int ir = 1; ir <= num_row(); ++ir) << 557 for(G4int ir=1;ir<=num_row();ir++)
599 { 558 {
600 for(G4int ic = 1; ic <= ir; ++ic) << 559 for(G4int ic=1;ic<=ir;ic++)
601 { 560 {
602 *(b++) = (*f)(*(a++), ir, ic); 561 *(b++) = (*f)(*(a++), ir, ic);
603 } 562 }
604 } 563 }
605 return mret; 564 return mret;
606 } 565 }
607 566
608 void G4ErrorSymMatrix::assign(const G4ErrorMat << 567 void G4ErrorSymMatrix::assign (const G4ErrorMatrix &m1)
609 { 568 {
610 if(mat1.nrow != nrow) << 569 if(m1.nrow != nrow)
611 { << 570 {
612 nrow = mat1.nrow; << 571 nrow = m1.nrow;
613 size = nrow * (nrow + 1) / 2; << 572 size = nrow * (nrow+1) / 2;
614 m.resize(size); << 573 m.resize(size);
615 } << 574 }
616 G4ErrorMatrixConstIter a = mat1.m.begin(); << 575 G4ErrorMatrixConstIter a = m1.m.begin();
617 G4ErrorMatrixIter b = m.begin(); << 576 G4ErrorMatrixIter b = m.begin();
618 for(G4int r = 1; r <= nrow; r++) << 577 for(G4int r=1;r<=nrow;r++)
619 { << 578 {
620 G4ErrorMatrixConstIter d = a; << 579 G4ErrorMatrixConstIter d = a;
621 for(G4int c = 1; c <= r; c++) << 580 for(G4int c=1;c<=r;c++)
622 { << 581 {
623 *(b++) = *(d++); << 582 *(b++) = *(d++);
624 } << 583 }
625 a += nrow; << 584 a += nrow;
626 } << 585 }
627 } 586 }
628 587
629 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 588 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorMatrix &m1) const
630 { 589 {
631 G4ErrorSymMatrix mret(mat1.num_row()); << 590 G4ErrorSymMatrix mret(m1.num_row());
632 G4ErrorMatrix temp = mat1 * (*this); << 591 G4ErrorMatrix temp = m1*(*this);
633 592
634 // If mat1*(*this) has correct dimensions, t << 593 // If m1*(*this) has correct dimensions, then so will the m1.T multiplication.
635 // multiplication. So there is no need to ch << 594 // So there is no need to check dimensions again.
636 595
637 G4int n = mat1.num_col(); << 596 G4int n = m1.num_col();
638 G4ErrorMatrixIter mr = mret.m.begin(); << 597 G4ErrorMatrixIter mr = mret.m.begin();
639 G4ErrorMatrixIter tempr1 = temp.m.begin(); 598 G4ErrorMatrixIter tempr1 = temp.m.begin();
640 for(G4int r = 1; r <= mret.num_row(); r++) << 599 for(G4int r=1;r<=mret.num_row();r++)
641 { 600 {
642 G4ErrorMatrixConstIter m1c1 = mat1.m.begin << 601 G4ErrorMatrixConstIter m1c1 = m1.m.begin();
643 for(G4int c = 1; c <= r; c++) << 602 for(G4int c=1;c<=r;c++)
644 { 603 {
645 G4double tmp = 0.0; << 604 G4double tmp = 0.0;
646 G4ErrorMatrixIter tempri = tempr1; << 605 G4ErrorMatrixIter tempri = tempr1;
647 G4ErrorMatrixConstIter m1ci = m1c1; 606 G4ErrorMatrixConstIter m1ci = m1c1;
648 for(G4int i = 1; i <= mat1.num_col(); i+ << 607 for(G4int i=1;i<=m1.num_col();i++)
649 { 608 {
650 tmp += (*(tempri++)) * (*(m1ci++)); << 609 tmp+=(*(tempri++))*(*(m1ci++));
651 } 610 }
652 *(mr++) = tmp; 611 *(mr++) = tmp;
653 m1c1 += n; 612 m1c1 += n;
654 } 613 }
655 tempr1 += n; 614 tempr1 += n;
656 } 615 }
657 return mret; 616 return mret;
658 } 617 }
659 618
660 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 619 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorSymMatrix &m1) const
661 const G4ErrorSymMatrix& mat1) const <<
662 { 620 {
663 G4ErrorSymMatrix mret(mat1.num_row()); << 621 G4ErrorSymMatrix mret(m1.num_row());
664 G4ErrorMatrix temp = mat1 * (*this); << 622 G4ErrorMatrix temp = m1*(*this);
665 G4int n = mat1.num_col(); << 623 G4int n = m1.num_col();
666 G4ErrorMatrixIter mr = mret.m.begin(); << 624 G4ErrorMatrixIter mr = mret.m.begin();
667 G4ErrorMatrixIter tempr1 = temp.m.begin(); 625 G4ErrorMatrixIter tempr1 = temp.m.begin();
668 for(G4int r = 1; r <= mret.num_row(); r++) << 626 for(G4int r=1;r<=mret.num_row();r++)
669 { 627 {
670 G4ErrorMatrixConstIter m1c1 = mat1.m.begin << 628 G4ErrorMatrixConstIter m1c1 = m1.m.begin();
671 G4int c; 629 G4int c;
672 for(c = 1; c <= r; c++) << 630 for(c=1;c<=r;c++)
673 { 631 {
674 G4double tmp = 0.0; << 632 G4double tmp = 0.0;
675 G4ErrorMatrixIter tempri = tempr1; << 633 G4ErrorMatrixIter tempri = tempr1;
676 G4ErrorMatrixConstIter m1ci = m1c1; 634 G4ErrorMatrixConstIter m1ci = m1c1;
677 G4int i; 635 G4int i;
678 for(i = 1; i < c; i++) << 636 for(i=1;i<c;i++)
679 { 637 {
680 tmp += (*(tempri++)) * (*(m1ci++)); << 638 tmp+=(*(tempri++))*(*(m1ci++));
681 } 639 }
682 for(i = c; i <= mat1.num_col(); i++) << 640 for(i=c;i<=m1.num_col();i++)
683 { 641 {
684 tmp += (*(tempri++)) * (*(m1ci)); << 642 tmp+=(*(tempri++))*(*(m1ci));
685 m1ci += i; 643 m1ci += i;
686 } 644 }
687 *(mr++) = tmp; 645 *(mr++) = tmp;
688 m1c1 += c; 646 m1c1 += c;
689 } 647 }
690 tempr1 += n; 648 tempr1 += n;
691 } 649 }
692 return mret; 650 return mret;
693 } 651 }
694 652
695 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT << 653 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT(const G4ErrorMatrix &m1) const
696 { 654 {
697 G4ErrorSymMatrix mret(mat1.num_col()); << 655 G4ErrorSymMatrix mret(m1.num_col());
698 G4ErrorMatrix temp = (*this) * mat1; << 656 G4ErrorMatrix temp = (*this)*m1;
699 G4int n = mat1.num_col(); << 657 G4int n = m1.num_col();
700 G4ErrorMatrixIter mrc = mret.m.begin(); << 658 G4ErrorMatrixIter mrc = mret.m.begin();
701 G4ErrorMatrixIter temp1r = temp.m.begin(); 659 G4ErrorMatrixIter temp1r = temp.m.begin();
702 for(G4int r = 1; r <= mret.num_row(); r++) << 660 for(G4int r=1;r<=mret.num_row();r++)
703 { 661 {
704 G4ErrorMatrixConstIter m11c = mat1.m.begin << 662 G4ErrorMatrixConstIter m11c = m1.m.begin();
705 for(G4int c = 1; c <= r; c++) << 663 for(G4int c=1;c<=r;c++)
706 { 664 {
707 G4double tmp = 0.0; << 665 G4double tmp = 0.0;
708 G4ErrorMatrixIter tempir = temp1r; << 666 G4ErrorMatrixIter tempir = temp1r;
709 G4ErrorMatrixConstIter m1ic = m11c; 667 G4ErrorMatrixConstIter m1ic = m11c;
710 for(G4int i = 1; i <= mat1.num_row(); i+ << 668 for(G4int i=1;i<=m1.num_row();i++)
711 { 669 {
712 tmp += (*(tempir)) * (*(m1ic)); << 670 tmp+=(*(tempir))*(*(m1ic));
713 tempir += n; 671 tempir += n;
714 m1ic += n; 672 m1ic += n;
715 } 673 }
716 *(mrc++) = tmp; 674 *(mrc++) = tmp;
717 m11c++; 675 m11c++;
718 } 676 }
719 temp1r++; 677 temp1r++;
720 } 678 }
721 return mret; 679 return mret;
722 } 680 }
723 681
724 void G4ErrorSymMatrix::invert(G4int& ifail) << 682 void G4ErrorSymMatrix::invert(G4int &ifail)
725 { << 683 {
726 ifail = 0; 684 ifail = 0;
727 685
728 switch(nrow) 686 switch(nrow)
729 { 687 {
730 case 3: << 688 case 3:
731 { 689 {
732 G4double det, temp; 690 G4double det, temp;
733 G4double t1, t2, t3; 691 G4double t1, t2, t3;
734 G4double c11, c12, c13, c22, c23, c33; << 692 G4double c11,c12,c13,c22,c23,c33;
735 c11 = (*(m.begin() + 2)) * (*(m.begin() << 693 c11 = (*(m.begin()+2)) * (*(m.begin()+5))
736 (*(m.begin() + 4)) * (*(m.begin() << 694 - (*(m.begin()+4)) * (*(m.begin()+4));
737 c12 = (*(m.begin() + 4)) * (*(m.begin() << 695 c12 = (*(m.begin()+4)) * (*(m.begin()+3))
738 (*(m.begin() + 1)) * (*(m.begin() << 696 - (*(m.begin()+1)) * (*(m.begin()+5));
739 c13 = (*(m.begin() + 1)) * (*(m.begin() << 697 c13 = (*(m.begin()+1)) * (*(m.begin()+4))
740 (*(m.begin() + 2)) * (*(m.begin() << 698 - (*(m.begin()+2)) * (*(m.begin()+3));
741 c22 = (*(m.begin() + 5)) * (*m.begin()) << 699 c22 = (*(m.begin()+5)) * (*m.begin())
742 (*(m.begin() + 3)) * (*(m.begin() << 700 - (*(m.begin()+3)) * (*(m.begin()+3));
743 c23 = (*(m.begin() + 3)) * (*(m.begin() << 701 c23 = (*(m.begin()+3)) * (*(m.begin()+1))
744 (*(m.begin() + 4)) * (*m.begin()); << 702 - (*(m.begin()+4)) * (*m.begin());
745 c33 = (*m.begin()) * (*(m.begin() + 2)) << 703 c33 = (*m.begin()) * (*(m.begin()+2))
746 (*(m.begin() + 1)) * (*(m.begin() << 704 - (*(m.begin()+1)) * (*(m.begin()+1));
747 t1 = std::fabs(*m.begin()); 705 t1 = std::fabs(*m.begin());
748 t2 = std::fabs(*(m.begin() + 1)); << 706 t2 = std::fabs(*(m.begin()+1));
749 t3 = std::fabs(*(m.begin() + 3)); << 707 t3 = std::fabs(*(m.begin()+3));
750 if(t1 >= t2) << 708 if (t1 >= t2)
751 { 709 {
752 if(t3 >= t1) << 710 if (t3 >= t1)
753 { 711 {
754 temp = *(m.begin() + 3); << 712 temp = *(m.begin()+3);
755 det = c23 * c12 - c22 * c13; << 713 det = c23*c12-c22*c13;
756 } 714 }
757 else 715 else
758 { 716 {
759 temp = *m.begin(); 717 temp = *m.begin();
760 det = c22 * c33 - c23 * c23; << 718 det = c22*c33-c23*c23;
761 } 719 }
762 } 720 }
763 else if(t3 >= t2) << 721 else if (t3 >= t2)
764 { 722 {
765 temp = *(m.begin() + 3); << 723 temp = *(m.begin()+3);
766 det = c23 * c12 - c22 * c13; << 724 det = c23*c12-c22*c13;
767 } 725 }
768 else 726 else
769 { 727 {
770 temp = *(m.begin() + 1); << 728 temp = *(m.begin()+1);
771 det = c13 * c23 - c12 * c33; << 729 det = c13*c23-c12*c33;
772 } 730 }
773 if(det == 0) << 731 if (det==0)
774 { 732 {
775 ifail = 1; 733 ifail = 1;
776 return; 734 return;
777 } 735 }
778 { 736 {
779 G4double ss = temp / det; << 737 G4double s = temp/det;
780 G4ErrorMatrixIter mq = m.begin(); << 738 G4ErrorMatrixIter mm = m.begin();
781 *(mq++) = ss * c11; << 739 *(mm++) = s*c11;
782 *(mq++) = ss * c12; << 740 *(mm++) = s*c12;
783 *(mq++) = ss * c22; << 741 *(mm++) = s*c22;
784 *(mq++) = ss * c13; << 742 *(mm++) = s*c13;
785 *(mq++) = ss * c23; << 743 *(mm++) = s*c23;
786 *(mq) = ss * c33; << 744 *(mm) = s*c33;
787 } 745 }
788 } 746 }
789 break; 747 break;
790 case 2: << 748 case 2:
791 { 749 {
792 G4double det, temp, ss; << 750 G4double det, temp, s;
793 det = (*m.begin()) * (*(m.begin() + 2)) << 751 det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
794 (*(m.begin() + 1)) * (*(m.begin() << 752 if (det==0)
795 if(det == 0) <<
796 { 753 {
797 ifail = 1; 754 ifail = 1;
798 return; 755 return;
799 } 756 }
800 ss = 1.0 / det; << 757 s = 1.0/det;
801 *(m.begin() + 1) *= -ss; << 758 *(m.begin()+1) *= -s;
802 temp = ss * (*(m.begin() + 2 << 759 temp = s*(*(m.begin()+2));
803 *(m.begin() + 2) = ss * (*m.begin()); << 760 *(m.begin()+2) = s*(*m.begin());
804 *m.begin() = temp; << 761 *m.begin() = temp;
805 break; 762 break;
806 } 763 }
807 case 1: << 764 case 1:
808 { 765 {
809 if((*m.begin()) == 0) << 766 if ((*m.begin())==0)
810 { 767 {
811 ifail = 1; 768 ifail = 1;
812 return; 769 return;
813 } 770 }
814 *m.begin() = 1.0 / (*m.begin()); << 771 *m.begin() = 1.0/(*m.begin());
815 break; 772 break;
816 } 773 }
817 case 5: << 774 case 5:
818 { 775 {
819 invert5(ifail); 776 invert5(ifail);
820 return; 777 return;
821 } 778 }
822 case 6: << 779 case 6:
823 { 780 {
824 invert6(ifail); 781 invert6(ifail);
825 return; 782 return;
826 } 783 }
827 case 4: << 784 case 4:
828 { 785 {
829 invert4(ifail); 786 invert4(ifail);
830 return; 787 return;
831 } 788 }
832 default: << 789 default:
833 { 790 {
834 invertBunchKaufman(ifail); 791 invertBunchKaufman(ifail);
835 return; 792 return;
836 } 793 }
837 } 794 }
838 return; // inversion successful << 795 return; // inversion successful
839 } 796 }
840 797
841 G4double G4ErrorSymMatrix::determinant() const 798 G4double G4ErrorSymMatrix::determinant() const
842 { 799 {
843 static const G4int max_array = 20; 800 static const G4int max_array = 20;
844 801
845 // ir must point to an array which is ***1 l 802 // ir must point to an array which is ***1 longer than*** nrow
846 803
847 static std::vector<G4int> ir_vec(max_array + << 804 static std::vector<G4int> ir_vec (max_array+1);
848 if(ir_vec.size() <= static_cast<unsigned int << 805 if (ir_vec.size() <= static_cast<unsigned int>(nrow))
849 { 806 {
850 ir_vec.resize(nrow + 1); << 807 ir_vec.resize(nrow+1);
851 } 808 }
852 G4int* ir = &ir_vec[0]; << 809 G4int * ir = &ir_vec[0];
853 810
854 G4double det; 811 G4double det;
855 G4ErrorMatrix mt(*this); 812 G4ErrorMatrix mt(*this);
856 G4int i = mt.dfact_matrix(det, ir); 813 G4int i = mt.dfact_matrix(det, ir);
857 if(i == 0) << 814 if(i==0) { return det; }
858 { <<
859 return det; <<
860 } <<
861 return 0.0; 815 return 0.0;
862 } 816 }
863 817
864 G4double G4ErrorSymMatrix::trace() const 818 G4double G4ErrorSymMatrix::trace() const
865 { 819 {
866 G4double t = 0.0; << 820 G4double t = 0.0;
867 for(G4int i = 0; i < nrow; i++) << 821 for (G4int i=0; i<nrow; i++)
868 { << 822 { t += *(m.begin() + (i+3)*i/2); }
869 t += *(m.begin() + (i + 3) * i / 2); << 823 return t;
870 } <<
871 return t; <<
872 } 824 }
873 825
874 void G4ErrorSymMatrix::invertBunchKaufman(G4in << 826 void G4ErrorSymMatrix::invertBunchKaufman(G4int &ifail)
875 { 827 {
876 // Bunch-Kaufman diagonal pivoting method 828 // Bunch-Kaufman diagonal pivoting method
877 // It is decribed in J.R. Bunch, L. Kaufman << 829 // It is decribed in J.R. Bunch, L. Kaufman (1977).
878 // "Some Stable Methods for Calculating Iner << 830 // "Some Stable Methods for Calculating Inertia and Solving Symmetric
879 // Linear Systems", Math. Comp. 31, p. 162-1 << 831 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
880 // Charles F. van Loan, "Matrix Computations << 832 // Charles F. van Loan, "Matrix Computations" (the second edition
881 // has a bug.) and implemented in "lapack" 833 // has a bug.) and implemented in "lapack"
882 // Mario Stanke, 09/97 834 // Mario Stanke, 09/97
883 835
884 G4int i, j, k, ss; << 836 G4int i, j, k, s;
885 G4int pivrow; 837 G4int pivrow;
886 838
887 // Establish the two working-space arrays ne 839 // Establish the two working-space arrays needed: x and piv are
888 // used as pointers to arrays of doubles and 840 // used as pointers to arrays of doubles and ints respectively, each
889 // of length nrow. We do not want to reallo 841 // of length nrow. We do not want to reallocate each time through
890 // unless the size needs to grow. We do not 842 // unless the size needs to grow. We do not want to leak memory, even
891 // by having a new without a delete that is 843 // by having a new without a delete that is only done once.
892 << 844
893 static const G4int max_array << 845 static const G4int max_array = 25;
894 static G4ThreadLocal std::vector<G4double>* << 846 static std::vector<G4double> xvec (max_array);
895 if(!xvec) << 847 static std::vector<G4int> pivv (max_array);
896 xvec = new std::vector<G4double>(max_array << 848 typedef std::vector<G4int>::iterator pivIter;
897 static G4ThreadLocal std::vector<G4int>* piv << 849 if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
898 if(!pivv) << 850 if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
899 pivv = new std::vector<G4int>(max_array); << 851 // Note - resize should do nothing if the size is already larger than nrow,
900 typedef std::vector<G4int>::iterator pivIter << 852 // but on VC++ there are indications that it does so we check.
901 if(xvec->size() < static_cast<unsigned int>( << 853 // Note - the data elements in a vector are guaranteed to be contiguous,
902 xvec->resize(nrow); << 854 // so x[i] and piv[i] are optimally fast.
903 if(pivv->size() < static_cast<unsigned int>( << 855 G4ErrorMatrixIter x = xvec.begin();
904 pivv->resize(nrow); << 856 // x[i] is used as helper storage, needs to have at least size nrow.
905 // Note - resize should do nothing if the s << 857 pivIter piv = pivv.begin();
906 // but on VC++ there are indications << 858 // piv[i] is used to store details of exchanges
907 // Note - the data elements in a vector are << 859
908 // so x[i] and piv[i] are optimally f <<
909 G4ErrorMatrixIter x = xvec->begin(); <<
910 // x[i] is used as helper storage, needs to <<
911 pivIter piv = pivv->begin(); <<
912 // piv[i] is used to store details of exchan <<
913 <<
914 G4double temp1, temp2; 860 G4double temp1, temp2;
915 G4ErrorMatrixIter ip, mjj, iq; 861 G4ErrorMatrixIter ip, mjj, iq;
916 G4double lambda, sigma; 862 G4double lambda, sigma;
917 const G4double alpha = .6404; // = (1+sqr << 863 const G4double alpha = .6404; // = (1+sqrt(17))/8
918 const G4double epsilon = 32 * DBL_EPSILON; << 864 const G4double epsilon = 32*DBL_EPSILON;
919 // whenever a sum of two doubles is below or << 865 // whenever a sum of two doubles is below or equal to epsilon
920 // it is set to zero. << 866 // it is set to zero.
921 // this constant could be set to zero but th << 867 // this constant could be set to zero but then the algorithm
922 // doesn't neccessarily detect that a matrix << 868 // doesn't neccessarily detect that a matrix is singular
923 << 869
924 for(i = 0; i < nrow; i++) << 870 for (i = 0; i < nrow; i++)
925 { 871 {
926 piv[i] = i + 1; << 872 piv[i] = i+1;
927 } 873 }
928 << 874
929 ifail = 0; 875 ifail = 0;
930 << 876
931 // compute the factorization P*A*P^T = L * D << 877 // compute the factorization P*A*P^T = L * D * L^T
932 // L is unit lower triangular, D is direct s 878 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
933 // L and D^-1 are stored in A = *this, P is 879 // L and D^-1 are stored in A = *this, P is stored in piv[]
934 << 880
935 for(j = 1; j < nrow; j += ss) // main loop << 881 for (j=1; j < nrow; j+=s) // main loop over columns
936 { 882 {
937 mjj = m.begin() + j * (j - 1) / 2 + j - << 883 mjj = m.begin() + j*(j-1)/2 + j-1;
938 lambda = 0; // compute lambda = max of A( << 884 lambda = 0; // compute lambda = max of A(j+1:n,j)
939 pivrow = j + 1; << 885 pivrow = j+1;
940 ip = m.begin() + (j + 1) * j / 2 + j - << 886 ip = m.begin() + (j+1)*j/2 + j-1;
941 for(i = j + 1; i <= nrow; ip += i++) << 887 for (i=j+1; i <= nrow ; ip += i++)
942 { <<
943 if(std::fabs(*ip) > lambda) <<
944 { <<
945 lambda = std::fabs(*ip); <<
946 pivrow = i; <<
947 } <<
948 } <<
949 if(lambda == 0) <<
950 { <<
951 if(*mjj == 0) <<
952 { <<
953 ifail = 1; <<
954 return; <<
955 } <<
956 ss = 1; <<
957 *mjj = 1. / *mjj; <<
958 } <<
959 else <<
960 { <<
961 if(std::fabs(*mjj) >= lambda * alpha) <<
962 { <<
963 ss = 1; <<
964 pivrow = j; <<
965 } <<
966 else <<
967 { <<
968 sigma = 0; // compute sigma = max A(p <<
969 ip = m.begin() + pivrow * (pivrow - <<
970 for(k = j; k < pivrow; k++) <<
971 { <<
972 if(std::fabs(*ip) > sigma) <<
973 sigma = std::fabs(*ip); <<
974 ip++; <<
975 } <<
976 if(sigma * std::fabs(*mjj) >= alpha * <<
977 { <<
978 ss = 1; <<
979 pivrow = j; <<
980 } <<
981 else if(std::fabs(*(m.begin() + pivrow <<
982 1)) >= alpha * sig <<
983 { <<
984 ss = 1; <<
985 } <<
986 else <<
987 { <<
988 ss = 2; <<
989 } <<
990 } <<
991 if(pivrow == j) // no permutation necce <<
992 { <<
993 piv[j - 1] = pivrow; <<
994 if(*mjj == 0) <<
995 { <<
996 ifail = 1; <<
997 return; <<
998 } <<
999 temp2 = *mjj = 1. / *mjj; // invert D <<
1000 <<
1001 // update A(j+1:n, j+1,n) <<
1002 for(i = j + 1; i <= nrow; i++) <<
1003 { <<
1004 temp1 = *(m.begin() + i * (i - 1) / <<
1005 ip = m.begin() + i * (i - 1) / 2 <<
1006 for(k = j + 1; k <= i; k++) <<
1007 { <<
1008 *ip -= temp1 * *(m.begin() + k * <<
1009 if(std::fabs(*ip) <= epsilon) <<
1010 { <<
1011 *ip = 0; <<
1012 } <<
1013 ip++; <<
1014 } <<
1015 } <<
1016 // update L <<
1017 ip = m.begin() + (j + 1) * j / 2 + j <<
1018 for(i = j + 1; i <= nrow; ip += i++) <<
1019 { <<
1020 *ip *= temp2; <<
1021 } <<
1022 } <<
1023 else if(ss == 1) // 1x1 pivot <<
1024 { <<
1025 piv[j - 1] = pivrow; <<
1026 <<
1027 // interchange rows and columns j and <<
1028 // submatrix (j:n,j:n) <<
1029 ip = m.begin() + pivrow * (pivrow - 1 <<
1030 for(i = j + 1; i < pivrow; i++, ip++) <<
1031 { <<
1032 temp1 = *(m.begin() + i * (i - 1) / <<
1033 *(m.begin() + i * (i - 1) / 2 + j - <<
1034 *ip <<
1035 } <<
1036 temp1 = *mjj; <<
1037 *mjj = *(m.begin() + pivrow * (pivro <<
1038 *(m.begin() + pivrow * (pivrow - 1) / <<
1039 ip = m.begin() + (pivrow + 1) * pivro <<
1040 iq = ip + pivrow - j; <<
1041 for(i = pivrow + 1; i <= nrow; ip += <<
1042 { <<
1043 temp1 = *iq; <<
1044 *iq = *ip; <<
1045 *ip = temp1; <<
1046 } <<
1047 <<
1048 if(*mjj == 0) <<
1049 { <<
1050 ifail = 1; <<
1051 return; <<
1052 } <<
1053 temp2 = *mjj = 1. / *mjj; // invert <<
1054 <<
1055 // update A(j+1:n, j+1:n) <<
1056 for(i = j + 1; i <= nrow; i++) <<
1057 { <<
1058 temp1 = *(m.begin() + i * (i - 1) / <<
1059 ip = m.begin() + i * (i - 1) / 2 <<
1060 for(k = j + 1; k <= i; k++) <<
1061 { 888 {
1062 *ip -= temp1 * *(m.begin() + k * << 889 if (std::fabs(*ip) > lambda)
1063 if(std::fabs(*ip) <= epsilon) <<
1064 { <<
1065 *ip = 0; <<
1066 } <<
1067 ip++; <<
1068 } <<
1069 } <<
1070 // update L <<
1071 ip = m.begin() + (j + 1) * j / 2 + j <<
1072 for(i = j + 1; i <= nrow; ip += i++) <<
1073 { <<
1074 *ip *= temp2; <<
1075 } <<
1076 } <<
1077 else // ss=2, ie use a 2x2 pivot <<
1078 { <<
1079 piv[j - 1] = -pivrow; <<
1080 piv[j] = 0; // that means this i <<
1081 <<
1082 if(j + 1 != pivrow) <<
1083 { <<
1084 // interchange rows and columns j+1 <<
1085 // submatrix (j:n,j:n) <<
1086 ip = m.begin() + pivrow * (pivrow - <<
1087 for(i = j + 2; i < pivrow; i++, ip+ <<
1088 { <<
1089 temp1 = *(m.begin() + i * (i - 1) <<
1090 *(m.begin() + i * (i - 1) / 2 + j <<
1091 *ip <<
1092 } <<
1093 temp1 = *(mjj + j + 1); <<
1094 *(mjj + j + 1) = <<
1095 *(m.begin() + pivrow * (pivrow - <<
1096 *(m.begin() + pivrow * (pivrow - 1) <<
1097 temp1 <<
1098 *(mjj + j) = *(m.begin() + pivrow * <<
1099 *(m.begin() + pivrow * (pivrow - 1) <<
1100 ip = m.begin() + (pivrow + 1) * piv <<
1101 iq = ip + pivrow - (j + 1); <<
1102 for(i = pivrow + 1; i <= nrow; ip + <<
1103 { <<
1104 temp1 = *iq; <<
1105 *iq = *ip; <<
1106 *ip = temp1; <<
1107 } <<
1108 } <<
1109 // invert D(j:j+1,j:j+1) <<
1110 temp2 = *mjj * *(mjj + j + 1) - *(mjj <<
1111 if(temp2 == 0) <<
1112 { <<
1113 G4Exception("G4ErrorSymMatrix::bunc <<
1114 "GEANT4e-Notification", <<
1115 "Error in pivot choice! <<
1116 } <<
1117 temp2 = 1. / temp2; <<
1118 <<
1119 // this quotient is guaranteed to exi <<
1120 // of the pivot <<
1121 <<
1122 temp1 = *mjj; <<
1123 *mjj = *(mjj + j + 1) * tem <<
1124 *(mjj + j + 1) = temp1 * temp2; <<
1125 *(mjj + j) = -*(mjj + j) * temp2; <<
1126 <<
1127 if(j < nrow - 1) // otherwise do not <<
1128 { <<
1129 // update A(j+2:n, j+2:n) <<
1130 for(i = j + 2; i <= nrow; i++) <<
1131 { <<
1132 ip = m.begin() + i * (i - 1) / <<
1133 temp1 = *ip * *mjj + *(ip + 1) * <<
1134 if(std::fabs(temp1) <= epsilon) <<
1135 { <<
1136 temp1 = 0; <<
1137 } <<
1138 temp2 = *ip * *(mjj + j) + *(ip + <<
1139 if(std::fabs(temp2) <= epsilon) <<
1140 { <<
1141 temp2 = 0; <<
1142 } <<
1143 for(k = j + 2; k <= i; k++) <<
1144 { <<
1145 ip = m.begin() + i * (i - 1) / <<
1146 iq = m.begin() + k * (k - 1) / <<
1147 *ip -= temp1 * *iq + temp2 * *( <<
1148 if(std::fabs(*ip) <= epsilon) <<
1149 { 890 {
1150 *ip = 0; << 891 lambda = std::fabs(*ip);
>> 892 pivrow = i;
1151 } 893 }
1152 } <<
1153 } 894 }
1154 // update L << 895 if (lambda == 0 )
1155 for(i = j + 2; i <= nrow; i++) <<
1156 { <<
1157 ip = m.begin() + i * (i - 1) / <<
1158 temp1 = *ip * *mjj + *(ip + 1) * <<
1159 if(std::fabs(temp1) <= epsilon) <<
1160 { 896 {
1161 temp1 = 0; << 897 if (*mjj == 0)
>> 898 {
>> 899 ifail = 1;
>> 900 return;
>> 901 }
>> 902 s=1;
>> 903 *mjj = 1./ *mjj;
1162 } 904 }
1163 *(ip + 1) = *ip * *(mjj + j) + *( << 905 else
1164 if(std::fabs(*(ip + 1)) <= epsilo <<
1165 { 906 {
1166 *(ip + 1) = 0; << 907 if (std::fabs(*mjj) >= lambda*alpha)
>> 908 {
>> 909 s=1;
>> 910 pivrow=j;
>> 911 }
>> 912 else
>> 913 {
>> 914 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
>> 915 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
>> 916 for (k=j; k < pivrow; k++)
>> 917 {
>> 918 if (std::fabs(*ip) > sigma)
>> 919 sigma = std::fabs(*ip);
>> 920 ip++;
>> 921 }
>> 922 if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
>> 923 {
>> 924 s=1;
>> 925 pivrow = j;
>> 926 }
>> 927 else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
>> 928 >= alpha * sigma)
>> 929 { s=1; }
>> 930 else
>> 931 { s=2; }
>> 932 }
>> 933 if (pivrow == j) // no permutation neccessary
>> 934 {
>> 935 piv[j-1] = pivrow;
>> 936 if (*mjj == 0)
>> 937 {
>> 938 ifail=1;
>> 939 return;
>> 940 }
>> 941 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 942
>> 943 // update A(j+1:n, j+1,n)
>> 944 for (i=j+1; i <= nrow; i++)
>> 945 {
>> 946 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 947 ip = m.begin()+i*(i-1)/2+j;
>> 948 for (k=j+1; k<=i; k++)
>> 949 {
>> 950 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 951 if (std::fabs(*ip) <= epsilon)
>> 952 { *ip=0; }
>> 953 ip++;
>> 954 }
>> 955 }
>> 956 // update L
>> 957 ip = m.begin() + (j+1)*j/2 + j-1;
>> 958 for (i=j+1; i <= nrow; ip += i++)
>> 959 {
>> 960 *ip *= temp2;
>> 961 }
>> 962 }
>> 963 else if (s==1) // 1x1 pivot
>> 964 {
>> 965 piv[j-1] = pivrow;
>> 966
>> 967 // interchange rows and columns j and pivrow in
>> 968 // submatrix (j:n,j:n)
>> 969 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
>> 970 for (i=j+1; i < pivrow; i++, ip++)
>> 971 {
>> 972 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
>> 973 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
>> 974 *ip = temp1;
>> 975 }
>> 976 temp1 = *mjj;
>> 977 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
>> 978 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
>> 979 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
>> 980 iq = ip + pivrow-j;
>> 981 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 982 {
>> 983 temp1 = *iq;
>> 984 *iq = *ip;
>> 985 *ip = temp1;
>> 986 }
>> 987
>> 988 if (*mjj == 0)
>> 989 {
>> 990 ifail = 1;
>> 991 return;
>> 992 }
>> 993 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 994
>> 995 // update A(j+1:n, j+1:n)
>> 996 for (i = j+1; i <= nrow; i++)
>> 997 {
>> 998 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 999 ip = m.begin()+i*(i-1)/2+j;
>> 1000 for (k=j+1; k<=i; k++)
>> 1001 {
>> 1002 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 1003 if (std::fabs(*ip) <= epsilon)
>> 1004 { *ip=0; }
>> 1005 ip++;
>> 1006 }
>> 1007 }
>> 1008 // update L
>> 1009 ip = m.begin() + (j+1)*j/2 + j-1;
>> 1010 for (i=j+1; i<=nrow; ip += i++)
>> 1011 {
>> 1012 *ip *= temp2;
>> 1013 }
>> 1014 }
>> 1015 else // s=2, ie use a 2x2 pivot
>> 1016 {
>> 1017 piv[j-1] = -pivrow;
>> 1018 piv[j] = 0; // that means this is the second row of a 2x2 pivot
>> 1019
>> 1020 if (j+1 != pivrow)
>> 1021 {
>> 1022 // interchange rows and columns j+1 and pivrow in
>> 1023 // submatrix (j:n,j:n)
>> 1024 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
>> 1025 for (i=j+2; i < pivrow; i++, ip++)
>> 1026 {
>> 1027 temp1 = *(m.begin() + i*(i-1)/2 + j);
>> 1028 *(m.begin() + i*(i-1)/2 + j) = *ip;
>> 1029 *ip = temp1;
>> 1030 }
>> 1031 temp1 = *(mjj + j + 1);
>> 1032 *(mjj + j + 1) =
>> 1033 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
>> 1034 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
>> 1035 temp1 = *(mjj + j);
>> 1036 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
>> 1037 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
>> 1038 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
>> 1039 iq = ip + pivrow-(j+1);
>> 1040 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 1041 {
>> 1042 temp1 = *iq;
>> 1043 *iq = *ip;
>> 1044 *ip = temp1;
>> 1045 }
>> 1046 }
>> 1047 // invert D(j:j+1,j:j+1)
>> 1048 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
>> 1049 if (temp2 == 0)
>> 1050 {
>> 1051 G4cerr
>> 1052 << "G4ErrorSymMatrix::bunch_invert: error in pivot choice"
>> 1053 << G4endl;
>> 1054 }
>> 1055 temp2 = 1. / temp2;
>> 1056
>> 1057 // this quotient is guaranteed to exist by the choice
>> 1058 // of the pivot
>> 1059
>> 1060 temp1 = *mjj;
>> 1061 *mjj = *(mjj + j + 1) * temp2;
>> 1062 *(mjj + j + 1) = temp1 * temp2;
>> 1063 *(mjj + j) = - *(mjj + j) * temp2;
>> 1064
>> 1065 if (j < nrow-1) // otherwise do nothing
>> 1066 {
>> 1067 // update A(j+2:n, j+2:n)
>> 1068 for (i=j+2; i <= nrow ; i++)
>> 1069 {
>> 1070 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1071 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
>> 1072 if (std::fabs(temp1 ) <= epsilon)
>> 1073 { temp1 = 0; }
>> 1074 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
>> 1075 if (std::fabs(temp2 ) <= epsilon)
>> 1076 { temp2 = 0; }
>> 1077 for (k = j+2; k <= i ; k++)
>> 1078 {
>> 1079 ip = m.begin() + i*(i-1)/2 + k-1;
>> 1080 iq = m.begin() + k*(k-1)/2 + j-1;
>> 1081 *ip -= temp1 * *iq + temp2 * *(iq+1);
>> 1082 if (std::fabs(*ip) <= epsilon)
>> 1083 { *ip = 0; }
>> 1084 }
>> 1085 }
>> 1086 // update L
>> 1087 for (i=j+2; i <= nrow ; i++)
>> 1088 {
>> 1089 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1090 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
>> 1091 if (std::fabs(temp1) <= epsilon)
>> 1092 { temp1 = 0; }
>> 1093 *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
>> 1094 if (std::fabs(*(ip+1)) <= epsilon)
>> 1095 { *(ip+1) = 0; }
>> 1096 *ip = temp1;
>> 1097 }
>> 1098 }
>> 1099 }
1167 } 1100 }
1168 *ip = temp1; << 1101 } // end of main loop over columns
1169 } <<
1170 } <<
1171 } <<
1172 } <<
1173 } // end of main loop over columns <<
1174 1102
1175 if(j == nrow) // the the last pivot is 1x1 << 1103 if (j == nrow) // the the last pivot is 1x1
1176 { 1104 {
1177 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1105 mjj = m.begin() + j*(j-1)/2 + j-1;
1178 if(*mjj == 0) << 1106 if (*mjj == 0)
1179 { 1107 {
1180 ifail = 1; 1108 ifail = 1;
1181 return; 1109 return;
1182 } 1110 }
1183 else 1111 else
1184 { 1112 {
1185 *mjj = 1. / *mjj; 1113 *mjj = 1. / *mjj;
1186 } 1114 }
1187 } // end of last pivot code << 1115 } // end of last pivot code
1188 1116
1189 // computing the inverse from the factoriza 1117 // computing the inverse from the factorization
1190 << 1118
1191 for(j = nrow; j >= 1; j -= ss) // loop ove << 1119 for (j = nrow ; j >= 1 ; j -= s) // loop over columns
1192 { 1120 {
1193 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1121 mjj = m.begin() + j*(j-1)/2 + j-1;
1194 if(piv[j - 1] > 0) // 1x1 pivot, compute << 1122 if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
1195 { << 1123 {
1196 ss = 1; << 1124 s = 1;
1197 if(j < nrow) << 1125 if (j < nrow)
1198 { << 1126 {
1199 ip = m.begin() + (j + 1) * j / 2 + j << 1127 ip = m.begin() + (j+1)*j/2 + j-1;
1200 for(i = 0; i < nrow - j; ip += 1 + j << 1128 for (i=0; i < nrow-j; ip += 1+j+i++)
1201 { << 1129 {
1202 x[i] = *ip; << 1130 x[i] = *ip;
1203 } << 1131 }
1204 for(i = j + 1; i <= nrow; i++) << 1132 for (i=j+1; i<=nrow ; i++)
1205 { << 1133 {
1206 temp2 = 0; << 1134 temp2=0;
1207 ip = m.begin() + i * (i - 1) / 2 << 1135 ip = m.begin() + i*(i-1)/2 + j;
1208 for(k = 0; k <= i - j - 1; k++) << 1136 for (k=0; k <= i-j-1; k++)
1209 { << 1137 { temp2 += *ip++ * x[k]; }
1210 temp2 += *ip++ * x[k]; << 1138 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1211 } << 1139 { temp2 += *ip * x[k]; }
1212 for(ip += i - 1; k < nrow - j; ip + << 1140 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1213 { << 1141 }
1214 temp2 += *ip * x[k]; << 1142 temp2 = 0;
1215 } << 1143 ip = m.begin() + (j+1)*j/2 + j-1;
1216 *(m.begin() + i * (i - 1) / 2 + j - << 1144 for (k=0; k < nrow-j; ip += 1+j+k++)
1217 } << 1145 { temp2 += x[k] * *ip; }
1218 temp2 = 0; << 1146 *mjj -= temp2;
1219 ip = m.begin() + (j + 1) * j / 2 + << 1147 }
1220 for(k = 0; k < nrow - j; ip += 1 + j << 1148 }
1221 { << 1149 else //2x2 pivot, compute columns j and j-1 of the inverse
1222 temp2 += x[k] * *ip; << 1150 {
1223 } << 1151 if (piv[j-1] != 0)
1224 *mjj -= temp2; << 1152 { G4cerr << "error in piv" << piv[j-1] << G4endl; }
1225 } << 1153 s=2;
1226 } << 1154 if (j < nrow)
1227 else // 2x2 pivot, compute columns j and << 1155 {
1228 { << 1156 ip = m.begin() + (j+1)*j/2 + j-1;
1229 if(piv[j - 1] != 0) << 1157 for (i=0; i < nrow-j; ip += 1+j+i++)
1230 { << 1158 {
1231 std::ostringstream message; << 1159 x[i] = *ip;
1232 message << "Error in pivot: " << piv[ << 1160 }
1233 G4Exception("G4ErrorSymMatrix::invert << 1161 for (i=j+1; i<=nrow ; i++)
1234 "GEANT4e-Notification", J << 1162 {
1235 } << 1163 temp2 = 0;
1236 ss = 2; << 1164 ip = m.begin() + i*(i-1)/2 + j;
1237 if(j < nrow) << 1165 for (k=0; k <= i-j-1; k++)
1238 { << 1166 { temp2 += *ip++ * x[k]; }
1239 ip = m.begin() + (j + 1) * j / 2 + j << 1167 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1240 for(i = 0; i < nrow - j; ip += 1 + j << 1168 { temp2 += *ip * x[k]; }
1241 { << 1169 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1242 x[i] = *ip; << 1170 }
1243 } << 1171 temp2 = 0;
1244 for(i = j + 1; i <= nrow; i++) << 1172 ip = m.begin() + (j+1)*j/2 + j-1;
1245 { << 1173 for (k=0; k < nrow-j; ip += 1+j+k++)
1246 temp2 = 0; << 1174 { temp2 += x[k] * *ip; }
1247 ip = m.begin() + i * (i - 1) / 2 << 1175 *mjj -= temp2;
1248 for(k = 0; k <= i - j - 1; k++) << 1176 temp2 = 0;
1249 { << 1177 ip = m.begin() + (j+1)*j/2 + j-2;
1250 temp2 += *ip++ * x[k]; << 1178 for (i=j+1; i <= nrow; ip += i++)
1251 } << 1179 { temp2 += *ip * *(ip+1); }
1252 for(ip += i - 1; k < nrow - j; ip + << 1180 *(mjj-1) -= temp2;
1253 { << 1181 ip = m.begin() + (j+1)*j/2 + j-2;
1254 temp2 += *ip * x[k]; << 1182 for (i=0; i < nrow-j; ip += 1+j+i++)
1255 } << 1183 {
1256 *(m.begin() + i * (i - 1) / 2 + j - << 1184 x[i] = *ip;
1257 } << 1185 }
1258 temp2 = 0; << 1186 for (i=j+1; i <= nrow ; i++)
1259 ip = m.begin() + (j + 1) * j / 2 + << 1187 {
1260 for(k = 0; k < nrow - j; ip += 1 + j << 1188 temp2 = 0;
1261 { << 1189 ip = m.begin() + i*(i-1)/2 + j;
1262 temp2 += x[k] * *ip; << 1190 for (k=0; k <= i-j-1; k++)
1263 } << 1191 { temp2 += *ip++ * x[k]; }
1264 *mjj -= temp2; << 1192 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1265 temp2 = 0; << 1193 { temp2 += *ip * x[k]; }
1266 ip = m.begin() + (j + 1) * j / 2 + << 1194 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1267 for(i = j + 1; i <= nrow; ip += i++) << 1195 }
1268 { << 1196 temp2 = 0;
1269 temp2 += *ip * *(ip + 1); << 1197 ip = m.begin() + (j+1)*j/2 + j-2;
1270 } << 1198 for (k=0; k < nrow-j; ip += 1+j+k++)
1271 *(mjj - 1) -= temp2; << 1199 { temp2 += x[k] * *ip; }
1272 ip = m.begin() + (j + 1) * j / 2 + j << 1200 *(mjj-j) -= temp2;
1273 for(i = 0; i < nrow - j; ip += 1 + j << 1201 }
1274 { << 1202 }
1275 x[i] = *ip; << 1203
1276 } << 1204 // interchange rows and columns j and piv[j-1]
1277 for(i = j + 1; i <= nrow; i++) << 1205 // or rows and columns j and -piv[j-2]
1278 { << 1206
1279 temp2 = 0; << 1207 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1280 ip = m.begin() + i * (i - 1) / 2 << 1208 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1281 for(k = 0; k <= i - j - 1; k++) << 1209 for (i=j+1;i < pivrow; i++, ip++)
1282 { << 1210 {
1283 temp2 += *ip++ * x[k]; << 1211 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1284 } << 1212 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1285 for(ip += i - 1; k < nrow - j; ip + << 1213 *ip = temp1;
1286 { << 1214 }
1287 temp2 += *ip * x[k]; << 1215 temp1 = *mjj;
1288 } << 1216 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1289 *(m.begin() + i * (i - 1) / 2 + j - << 1217 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1290 } << 1218 if (s==2)
1291 temp2 = 0; << 1219 {
1292 ip = m.begin() + (j + 1) * j / 2 + << 1220 temp1 = *(mjj-1);
1293 for(k = 0; k < nrow - j; ip += 1 + j << 1221 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1294 { << 1222 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1295 temp2 += x[k] * *ip; << 1223 }
1296 } << 1224
1297 *(mjj - j) -= temp2; << 1225 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1298 } << 1226 iq = ip + pivrow-j;
1299 } << 1227 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1300 << 1228 {
1301 // interchange rows and columns j and piv << 1229 temp1 = *iq;
1302 // or rows and columns j and -piv[j-2] << 1230 *iq = *ip;
1303 << 1231 *ip = temp1;
1304 pivrow = (piv[j - 1] == 0) ? -piv[j - 2] << 1232 }
1305 ip = m.begin() + pivrow * (pivrow - 1 << 1233 } // end of loop over columns (in computing inverse from factorization)
1306 for(i = j + 1; i < pivrow; i++, ip++) << 1234
1307 { << 1235 return; // inversion successful
1308 temp1 = *(m.begin() + i * (i - 1) / 2 + << 1236 }
1309 *(m.begin() + i * (i - 1) / 2 + j - 1) << 1237
1310 *ip << 1238 G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
1311 } << 1239 G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
1312 temp1 = *mjj; << 1240 G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
1313 *mjj = *(m.begin() + pivrow * (pivrow - << 1241 G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
1314 *(m.begin() + pivrow * (pivrow - 1) / 2 + << 1242 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
1315 if(ss == 2) << 1243 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
1316 { << 1244 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
1317 temp1 = *(mjj - 1); << 1245 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;
1318 *(mjj - 1) = *(m.begin() + pivrow * (pi <<
1319 *(m.begin() + pivrow * (pivrow - 1) / 2 <<
1320 } <<
1321 <<
1322 ip = m.begin() + (pivrow + 1) * pivrow / <<
1323 iq = ip + pivrow - j; <<
1324 for(i = pivrow + 1; i <= nrow; ip += i, i <<
1325 { <<
1326 temp1 = *iq; <<
1327 *iq = *ip; <<
1328 *ip = temp1; <<
1329 } <<
1330 } // end of loop over columns (in computin <<
1331 <<
1332 return; // inversion successful <<
1333 } <<
1334 <<
1335 G4ThreadLocal G4double G4ErrorSymMatrix::posD <<
1336 G4ThreadLocal G4double G4ErrorSymMatrix::posD <<
1337 G4ThreadLocal G4double G4ErrorSymMatrix::adju <<
1338 G4ThreadLocal G4double G4ErrorSymMatrix::adju <<
1339 const G4double G4ErrorSymMatrix::CHOLESKY_THR <<
1340 const G4double G4ErrorSymMatrix::CHOLESKY_THR <<
1341 const G4double G4ErrorSymMatrix::CHOLESKY_CRE <<
1342 const G4double G4ErrorSymMatrix::CHOLESKY_CRE <<
1343 1246
1344 // Aij are indices for a 6x6 symmetric matrix 1247 // Aij are indices for a 6x6 symmetric matrix.
1345 // The indices for 5x5 or 4x4 symmetric m << 1248 // The indices for 5x5 or 4x4 symmetric matrices are the same,
1346 // ignoring all combinations with an inde 1249 // ignoring all combinations with an index which is inapplicable.
1347 1250
1348 #define A00 0 1251 #define A00 0
1349 #define A01 1 1252 #define A01 1
1350 #define A02 3 1253 #define A02 3
1351 #define A03 6 1254 #define A03 6
1352 #define A04 10 1255 #define A04 10
1353 #define A05 15 1256 #define A05 15
1354 1257
1355 #define A10 1 1258 #define A10 1
1356 #define A11 2 1259 #define A11 2
1357 #define A12 4 1260 #define A12 4
1358 #define A13 7 1261 #define A13 7
1359 #define A14 11 1262 #define A14 11
1360 #define A15 16 1263 #define A15 16
1361 1264
1362 #define A20 3 1265 #define A20 3
1363 #define A21 4 1266 #define A21 4
1364 #define A22 5 1267 #define A22 5
1365 #define A23 8 1268 #define A23 8
1366 #define A24 12 1269 #define A24 12
1367 #define A25 17 1270 #define A25 17
1368 1271
1369 #define A30 6 1272 #define A30 6
1370 #define A31 7 1273 #define A31 7
1371 #define A32 8 1274 #define A32 8
1372 #define A33 9 1275 #define A33 9
1373 #define A34 13 1276 #define A34 13
1374 #define A35 18 1277 #define A35 18
1375 1278
1376 #define A40 10 1279 #define A40 10
1377 #define A41 11 1280 #define A41 11
1378 #define A42 12 1281 #define A42 12
1379 #define A43 13 1282 #define A43 13
1380 #define A44 14 1283 #define A44 14
1381 #define A45 19 1284 #define A45 19
1382 1285
1383 #define A50 15 1286 #define A50 15
1384 #define A51 16 1287 #define A51 16
1385 #define A52 17 1288 #define A52 17
1386 #define A53 18 1289 #define A53 18
1387 #define A54 19 1290 #define A54 19
1388 #define A55 20 1291 #define A55 20
1389 1292
1390 void G4ErrorSymMatrix::invert5(G4int& ifail) << 1293 void G4ErrorSymMatrix::invert5(G4int & ifail)
1391 { << 1294 {
1392 if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_ << 1295 if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
1393 { << 1296 {
1394 invertCholesky5(ifail); << 1297 invertCholesky5(ifail);
1395 posDefFraction5x5 = .9 * posDefFraction5x << 1298 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
1396 if(ifail != 0) // Cholesky failed -- inv << 1299 if (ifail!=0) // Cholesky failed -- invert using Haywood
1397 { << 1300 {
1398 invertHaywood5(ifail); << 1301 invertHaywood5(ifail);
1399 } << 1302 }
1400 } << 1303 }
1401 else << 1304 else
1402 { <<
1403 if(posDefFraction5x5 + adjustment5x5 >= C <<
1404 { <<
1405 invertCholesky5(ifail); <<
1406 posDefFraction5x5 = .9 * posDefFraction <<
1407 if(ifail != 0) // Cholesky failed -- i <<
1408 { 1305 {
1409 invertHaywood5(ifail); << 1306 if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
1410 adjustment5x5 = 0; << 1307 {
>> 1308 invertCholesky5(ifail);
>> 1309 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
>> 1310 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1311 {
>> 1312 invertHaywood5(ifail);
>> 1313 adjustment5x5 = 0;
>> 1314 }
>> 1315 }
>> 1316 else
>> 1317 {
>> 1318 invertHaywood5(ifail);
>> 1319 adjustment5x5 += CHOLESKY_CREEP_5x5;
>> 1320 }
1411 } 1321 }
1412 } << 1322 return;
1413 else <<
1414 { <<
1415 invertHaywood5(ifail); <<
1416 adjustment5x5 += CHOLESKY_CREEP_5x5; <<
1417 } <<
1418 } <<
1419 return; <<
1420 } 1323 }
1421 1324
1422 void G4ErrorSymMatrix::invert6(G4int& ifail) << 1325 void G4ErrorSymMatrix::invert6(G4int & ifail)
1423 { 1326 {
1424 if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_ << 1327 if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
1425 { <<
1426 invertCholesky6(ifail); <<
1427 posDefFraction6x6 = .9 * posDefFraction6x <<
1428 if(ifail != 0) // Cholesky failed -- inv <<
1429 { <<
1430 invertHaywood6(ifail); <<
1431 } <<
1432 } <<
1433 else <<
1434 { <<
1435 if(posDefFraction6x6 + adjustment6x6 >= C <<
1436 { <<
1437 invertCholesky6(ifail); <<
1438 posDefFraction6x6 = .9 * posDefFraction <<
1439 if(ifail != 0) // Cholesky failed -- i <<
1440 { 1328 {
1441 invertHaywood6(ifail); << 1329 invertCholesky6(ifail);
1442 adjustment6x6 = 0; << 1330 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
>> 1331 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1332 {
>> 1333 invertHaywood6(ifail);
>> 1334 }
1443 } 1335 }
1444 } << 1336 else
1445 else << 1337 {
1446 { << 1338 if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
1447 invertHaywood6(ifail); << 1339 {
1448 adjustment6x6 += CHOLESKY_CREEP_6x6; << 1340 invertCholesky6(ifail);
1449 } << 1341 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
1450 } << 1342 if (ifail!=0) // Cholesky failed -- invert using Haywood
1451 return; << 1343 {
>> 1344 invertHaywood6(ifail);
>> 1345 adjustment6x6 = 0;
>> 1346 }
>> 1347 }
>> 1348 else
>> 1349 {
>> 1350 invertHaywood6(ifail);
>> 1351 adjustment6x6 += CHOLESKY_CREEP_6x6;
>> 1352 }
>> 1353 }
>> 1354 return;
1452 } 1355 }
1453 1356
1454 void G4ErrorSymMatrix::invertHaywood5(G4int& << 1357 void G4ErrorSymMatrix::invertHaywood5 (G4int & ifail)
1455 { 1358 {
1456 ifail = 0; 1359 ifail = 0;
1457 1360
1458 // Find all NECESSARY 2x2 dets: (25 of the 1361 // Find all NECESSARY 2x2 dets: (25 of them)
1459 1362
1460 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 1363 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
1461 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 1364 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
1462 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 1365 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
1463 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 1366 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
1464 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 1367 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
1465 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 1368 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
1466 G4double Det2_24_01 = m[A20] * m[A41] - m[A << 1369 G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
1467 G4double Det2_24_02 = m[A20] * m[A42] - m[A << 1370 G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
1468 G4double Det2_24_03 = m[A20] * m[A43] - m[A << 1371 G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
1469 G4double Det2_24_04 = m[A20] * m[A44] - m[A << 1372 G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
1470 G4double Det2_24_12 = m[A21] * m[A42] - m[A << 1373 G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
1471 G4double Det2_24_13 = m[A21] * m[A43] - m[A << 1374 G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
1472 G4double Det2_24_14 = m[A21] * m[A44] - m[A << 1375 G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
1473 G4double Det2_24_23 = m[A22] * m[A43] - m[A << 1376 G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
1474 G4double Det2_24_24 = m[A22] * m[A44] - m[A << 1377 G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
1475 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1378 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1476 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1379 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1477 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1380 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1478 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1381 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1479 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1382 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1480 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1383 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1481 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1384 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1482 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1385 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1483 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1386 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1484 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1387 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1485 1388
1486 // Find all NECESSARY 3x3 dets: (30 of th 1389 // Find all NECESSARY 3x3 dets: (30 of them)
1487 1390
1488 G4double Det3_123_012 = << 1391 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
1489 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 1392 + m[A12]*Det2_23_01;
1490 G4double Det3_123_013 = << 1393 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
1491 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 1394 + m[A13]*Det2_23_01;
1492 G4double Det3_123_023 = << 1395 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
1493 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 1396 + m[A13]*Det2_23_02;
1494 G4double Det3_123_123 = << 1397 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
1495 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 1398 + m[A13]*Det2_23_12;
1496 G4double Det3_124_012 = << 1399 G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02
1497 m[A10] * Det2_24_12 - m[A11] * Det2_24_02 << 1400 + m[A12]*Det2_24_01;
1498 G4double Det3_124_013 = << 1401 G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03
1499 m[A10] * Det2_24_13 - m[A11] * Det2_24_03 << 1402 + m[A13]*Det2_24_01;
1500 G4double Det3_124_014 = << 1403 G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04
1501 m[A10] * Det2_24_14 - m[A11] * Det2_24_04 << 1404 + m[A14]*Det2_24_01;
1502 G4double Det3_124_023 = << 1405 G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03
1503 m[A10] * Det2_24_23 - m[A12] * Det2_24_03 << 1406 + m[A13]*Det2_24_02;
1504 G4double Det3_124_024 = << 1407 G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04
1505 m[A10] * Det2_24_24 - m[A12] * Det2_24_04 << 1408 + m[A14]*Det2_24_02;
1506 G4double Det3_124_123 = << 1409 G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13
1507 m[A11] * Det2_24_23 - m[A12] * Det2_24_13 << 1410 + m[A13]*Det2_24_12;
1508 G4double Det3_124_124 = << 1411 G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14
1509 m[A11] * Det2_24_24 - m[A12] * Det2_24_14 << 1412 + m[A14]*Det2_24_12;
1510 G4double Det3_134_012 = << 1413 G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02
1511 m[A10] * Det2_34_12 - m[A11] * Det2_34_02 << 1414 + m[A12]*Det2_34_01;
1512 G4double Det3_134_013 = << 1415 G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03
1513 m[A10] * Det2_34_13 - m[A11] * Det2_34_03 << 1416 + m[A13]*Det2_34_01;
1514 G4double Det3_134_014 = << 1417 G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04
1515 m[A10] * Det2_34_14 - m[A11] * Det2_34_04 << 1418 + m[A14]*Det2_34_01;
1516 G4double Det3_134_023 = << 1419 G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03
1517 m[A10] * Det2_34_23 - m[A12] * Det2_34_03 << 1420 + m[A13]*Det2_34_02;
1518 G4double Det3_134_024 = << 1421 G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04
1519 m[A10] * Det2_34_24 - m[A12] * Det2_34_04 << 1422 + m[A14]*Det2_34_02;
1520 G4double Det3_134_034 = << 1423 G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04
1521 m[A10] * Det2_34_34 - m[A13] * Det2_34_04 << 1424 + m[A14]*Det2_34_03;
1522 G4double Det3_134_123 = << 1425 G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13
1523 m[A11] * Det2_34_23 - m[A12] * Det2_34_13 << 1426 + m[A13]*Det2_34_12;
1524 G4double Det3_134_124 = << 1427 G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14
1525 m[A11] * Det2_34_24 - m[A12] * Det2_34_14 << 1428 + m[A14]*Det2_34_12;
1526 G4double Det3_134_134 = << 1429 G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14
1527 m[A11] * Det2_34_34 - m[A13] * Det2_34_14 << 1430 + m[A14]*Det2_34_13;
1528 G4double Det3_234_012 = << 1431 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1529 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1432 + m[A22]*Det2_34_01;
1530 G4double Det3_234_013 = << 1433 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1531 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1434 + m[A23]*Det2_34_01;
1532 G4double Det3_234_014 = << 1435 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1533 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1436 + m[A24]*Det2_34_01;
1534 G4double Det3_234_023 = << 1437 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1535 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1438 + m[A23]*Det2_34_02;
1536 G4double Det3_234_024 = << 1439 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1537 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1440 + m[A24]*Det2_34_02;
1538 G4double Det3_234_034 = << 1441 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1539 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1442 + m[A24]*Det2_34_03;
1540 G4double Det3_234_123 = << 1443 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1541 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1444 + m[A23]*Det2_34_12;
1542 G4double Det3_234_124 = << 1445 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1543 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1446 + m[A24]*Det2_34_12;
1544 G4double Det3_234_134 = << 1447 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1545 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1448 + m[A24]*Det2_34_13;
1546 G4double Det3_234_234 = << 1449 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1547 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1450 + m[A24]*Det2_34_23;
1548 1451
1549 // Find all NECESSARY 4x4 dets: (15 of th 1452 // Find all NECESSARY 4x4 dets: (15 of them)
1550 1453
1551 G4double Det4_0123_0123 = m[A00] * Det3_123 << 1454 G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023
1552 m[A02] * Det3_123 << 1455 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
1553 G4double Det4_0124_0123 = m[A00] * Det3_124 << 1456 G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023
1554 m[A02] * Det3_124 << 1457 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
1555 G4double Det4_0124_0124 = m[A00] * Det3_124 << 1458 G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024
1556 m[A02] * Det3_124 << 1459 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
1557 G4double Det4_0134_0123 = m[A00] * Det3_134 << 1460 G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023
1558 m[A02] * Det3_134 << 1461 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
1559 G4double Det4_0134_0124 = m[A00] * Det3_134 << 1462 G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024
1560 m[A02] * Det3_134 << 1463 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
1561 G4double Det4_0134_0134 = m[A00] * Det3_134 << 1464 G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034
1562 m[A03] * Det3_134 << 1465 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
1563 G4double Det4_0234_0123 = m[A00] * Det3_234 << 1466 G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023
1564 m[A02] * Det3_234 << 1467 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
1565 G4double Det4_0234_0124 = m[A00] * Det3_234 << 1468 G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024
1566 m[A02] * Det3_234 << 1469 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
1567 G4double Det4_0234_0134 = m[A00] * Det3_234 << 1470 G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034
1568 m[A03] * Det3_234 << 1471 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
1569 G4double Det4_0234_0234 = m[A00] * Det3_234 << 1472 G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034
1570 m[A03] * Det3_234 << 1473 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
1571 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1474 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1572 m[A12] * Det3_234 << 1475 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1573 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1476 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1574 m[A12] * Det3_234 << 1477 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1575 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1478 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1576 m[A13] * Det3_234 << 1479 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1577 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1480 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1578 m[A13] * Det3_234 << 1481 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1579 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1482 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1580 m[A13] * Det3_234 << 1483 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1581 1484
1582 // Find the 5x5 det: 1485 // Find the 5x5 det:
1583 1486
1584 G4double det = m[A00] * Det4_1234_1234 - m[ << 1487 G4double det = m[A00]*Det4_1234_1234
1585 m[A02] * Det4_1234_0134 - m[ << 1488 - m[A01]*Det4_1234_0234
1586 m[A04] * Det4_1234_0123; << 1489 + m[A02]*Det4_1234_0134
>> 1490 - m[A03]*Det4_1234_0124
>> 1491 + m[A04]*Det4_1234_0123;
1587 1492
1588 if(det == 0) << 1493 if ( det == 0 )
1589 { << 1494 {
1590 ifail = 1; 1495 ifail = 1;
1591 return; 1496 return;
1592 } << 1497 }
1593 1498
1594 G4double oneOverDet = 1.0 / det; << 1499 G4double oneOverDet = 1.0/det;
1595 G4double mn1OverDet = -oneOverDet; << 1500 G4double mn1OverDet = - oneOverDet;
1596 1501
1597 m[A00] = Det4_1234_1234 * oneOverDet; << 1502 m[A00] = Det4_1234_1234 * oneOverDet;
1598 m[A01] = Det4_1234_0234 * mn1OverDet; << 1503 m[A01] = Det4_1234_0234 * mn1OverDet;
1599 m[A02] = Det4_1234_0134 * oneOverDet; << 1504 m[A02] = Det4_1234_0134 * oneOverDet;
1600 m[A03] = Det4_1234_0124 * mn1OverDet; << 1505 m[A03] = Det4_1234_0124 * mn1OverDet;
1601 m[A04] = Det4_1234_0123 * oneOverDet; << 1506 m[A04] = Det4_1234_0123 * oneOverDet;
1602 << 1507
1603 m[A11] = Det4_0234_0234 * oneOverDet; << 1508 m[A11] = Det4_0234_0234 * oneOverDet;
1604 m[A12] = Det4_0234_0134 * mn1OverDet; << 1509 m[A12] = Det4_0234_0134 * mn1OverDet;
1605 m[A13] = Det4_0234_0124 * oneOverDet; << 1510 m[A13] = Det4_0234_0124 * oneOverDet;
1606 m[A14] = Det4_0234_0123 * mn1OverDet; << 1511 m[A14] = Det4_0234_0123 * mn1OverDet;
1607 << 1512
1608 m[A22] = Det4_0134_0134 * oneOverDet; << 1513 m[A22] = Det4_0134_0134 * oneOverDet;
1609 m[A23] = Det4_0134_0124 * mn1OverDet; << 1514 m[A23] = Det4_0134_0124 * mn1OverDet;
1610 m[A24] = Det4_0134_0123 * oneOverDet; << 1515 m[A24] = Det4_0134_0123 * oneOverDet;
1611 1516
1612 m[A33] = Det4_0124_0124 * oneOverDet; << 1517 m[A33] = Det4_0124_0124 * oneOverDet;
1613 m[A34] = Det4_0124_0123 * mn1OverDet; << 1518 m[A34] = Det4_0124_0123 * mn1OverDet;
1614 1519
1615 m[A44] = Det4_0123_0123 * oneOverDet; << 1520 m[A44] = Det4_0123_0123 * oneOverDet;
1616 1521
1617 return; 1522 return;
1618 } 1523 }
1619 1524
1620 void G4ErrorSymMatrix::invertHaywood6(G4int& << 1525 void G4ErrorSymMatrix::invertHaywood6 (G4int & ifail)
1621 { 1526 {
1622 ifail = 0; 1527 ifail = 0;
1623 1528
1624 // Find all NECESSARY 2x2 dets: (39 of the 1529 // Find all NECESSARY 2x2 dets: (39 of them)
1625 1530
1626 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1531 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1627 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1532 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1628 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1533 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1629 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1534 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1630 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1535 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1631 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1536 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1632 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1537 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1633 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1538 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1634 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1539 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1635 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1540 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1636 G4double Det2_35_01 = m[A30] * m[A51] - m[A << 1541 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1637 G4double Det2_35_02 = m[A30] * m[A52] - m[A << 1542 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1638 G4double Det2_35_03 = m[A30] * m[A53] - m[A << 1543 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1639 G4double Det2_35_04 = m[A30] * m[A54] - m[A << 1544 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1640 G4double Det2_35_05 = m[A30] * m[A55] - m[A << 1545 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1641 G4double Det2_35_12 = m[A31] * m[A52] - m[A << 1546 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1642 G4double Det2_35_13 = m[A31] * m[A53] - m[A << 1547 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1643 G4double Det2_35_14 = m[A31] * m[A54] - m[A << 1548 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1644 G4double Det2_35_15 = m[A31] * m[A55] - m[A << 1549 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1645 G4double Det2_35_23 = m[A32] * m[A53] - m[A << 1550 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1646 G4double Det2_35_24 = m[A32] * m[A54] - m[A << 1551 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1647 G4double Det2_35_25 = m[A32] * m[A55] - m[A << 1552 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1648 G4double Det2_35_34 = m[A33] * m[A54] - m[A << 1553 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1649 G4double Det2_35_35 = m[A33] * m[A55] - m[A << 1554 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1650 G4double Det2_45_01 = m[A40] * m[A51] - m[A << 1555 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1651 G4double Det2_45_02 = m[A40] * m[A52] - m[A << 1556 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1652 G4double Det2_45_03 = m[A40] * m[A53] - m[A << 1557 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1653 G4double Det2_45_04 = m[A40] * m[A54] - m[A << 1558 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1654 G4double Det2_45_05 = m[A40] * m[A55] - m[A << 1559 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1655 G4double Det2_45_12 = m[A41] * m[A52] - m[A << 1560 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1656 G4double Det2_45_13 = m[A41] * m[A53] - m[A << 1561 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1657 G4double Det2_45_14 = m[A41] * m[A54] - m[A << 1562 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1658 G4double Det2_45_15 = m[A41] * m[A55] - m[A << 1563 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1659 G4double Det2_45_23 = m[A42] * m[A53] - m[A << 1564 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1660 G4double Det2_45_24 = m[A42] * m[A54] - m[A << 1565 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1661 G4double Det2_45_25 = m[A42] * m[A55] - m[A << 1566 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1662 G4double Det2_45_34 = m[A43] * m[A54] - m[A << 1567 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1663 G4double Det2_45_35 = m[A43] * m[A55] - m[A << 1568 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1664 G4double Det2_45_45 = m[A44] * m[A55] - m[A << 1569 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1665 1570
1666 // Find all NECESSARY 3x3 dets: (65 of the 1571 // Find all NECESSARY 3x3 dets: (65 of them)
1667 1572
1668 G4double Det3_234_012 = << 1573 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1669 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1574 + m[A22]*Det2_34_01;
1670 G4double Det3_234_013 = << 1575 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1671 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1576 + m[A23]*Det2_34_01;
1672 G4double Det3_234_014 = << 1577 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1673 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1578 + m[A24]*Det2_34_01;
1674 G4double Det3_234_023 = << 1579 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1675 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1580 + m[A23]*Det2_34_02;
1676 G4double Det3_234_024 = << 1581 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1677 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1582 + m[A24]*Det2_34_02;
1678 G4double Det3_234_034 = << 1583 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1679 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1584 + m[A24]*Det2_34_03;
1680 G4double Det3_234_123 = << 1585 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1681 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1586 + m[A23]*Det2_34_12;
1682 G4double Det3_234_124 = << 1587 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1683 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1588 + m[A24]*Det2_34_12;
1684 G4double Det3_234_134 = << 1589 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1685 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1590 + m[A24]*Det2_34_13;
1686 G4double Det3_234_234 = << 1591 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1687 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1592 + m[A24]*Det2_34_23;
1688 G4double Det3_235_012 = << 1593 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
1689 m[A20] * Det2_35_12 - m[A21] * Det2_35_02 << 1594 + m[A22]*Det2_35_01;
1690 G4double Det3_235_013 = << 1595 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
1691 m[A20] * Det2_35_13 - m[A21] * Det2_35_03 << 1596 + m[A23]*Det2_35_01;
1692 G4double Det3_235_014 = << 1597 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
1693 m[A20] * Det2_35_14 - m[A21] * Det2_35_04 << 1598 + m[A24]*Det2_35_01;
1694 G4double Det3_235_015 = << 1599 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
1695 m[A20] * Det2_35_15 - m[A21] * Det2_35_05 << 1600 + m[A25]*Det2_35_01;
1696 G4double Det3_235_023 = << 1601 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
1697 m[A20] * Det2_35_23 - m[A22] * Det2_35_03 << 1602 + m[A23]*Det2_35_02;
1698 G4double Det3_235_024 = << 1603 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
1699 m[A20] * Det2_35_24 - m[A22] * Det2_35_04 << 1604 + m[A24]*Det2_35_02;
1700 G4double Det3_235_025 = << 1605 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
1701 m[A20] * Det2_35_25 - m[A22] * Det2_35_05 << 1606 + m[A25]*Det2_35_02;
1702 G4double Det3_235_034 = << 1607 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
1703 m[A20] * Det2_35_34 - m[A23] * Det2_35_04 << 1608 + m[A24]*Det2_35_03;
1704 G4double Det3_235_035 = << 1609 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
1705 m[A20] * Det2_35_35 - m[A23] * Det2_35_05 << 1610 + m[A25]*Det2_35_03;
1706 G4double Det3_235_123 = << 1611 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
1707 m[A21] * Det2_35_23 - m[A22] * Det2_35_13 << 1612 + m[A23]*Det2_35_12;
1708 G4double Det3_235_124 = << 1613 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
1709 m[A21] * Det2_35_24 - m[A22] * Det2_35_14 << 1614 + m[A24]*Det2_35_12;
1710 G4double Det3_235_125 = << 1615 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
1711 m[A21] * Det2_35_25 - m[A22] * Det2_35_15 << 1616 + m[A25]*Det2_35_12;
1712 G4double Det3_235_134 = << 1617 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
1713 m[A21] * Det2_35_34 - m[A23] * Det2_35_14 << 1618 + m[A24]*Det2_35_13;
1714 G4double Det3_235_135 = << 1619 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
1715 m[A21] * Det2_35_35 - m[A23] * Det2_35_15 << 1620 + m[A25]*Det2_35_13;
1716 G4double Det3_235_234 = << 1621 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
1717 m[A22] * Det2_35_34 - m[A23] * Det2_35_24 << 1622 + m[A24]*Det2_35_23;
1718 G4double Det3_235_235 = << 1623 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
1719 m[A22] * Det2_35_35 - m[A23] * Det2_35_25 << 1624 + m[A25]*Det2_35_23;
1720 G4double Det3_245_012 = << 1625 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
1721 m[A20] * Det2_45_12 - m[A21] * Det2_45_02 << 1626 + m[A22]*Det2_45_01;
1722 G4double Det3_245_013 = << 1627 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
1723 m[A20] * Det2_45_13 - m[A21] * Det2_45_03 << 1628 + m[A23]*Det2_45_01;
1724 G4double Det3_245_014 = << 1629 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
1725 m[A20] * Det2_45_14 - m[A21] * Det2_45_04 << 1630 + m[A24]*Det2_45_01;
1726 G4double Det3_245_015 = << 1631 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
1727 m[A20] * Det2_45_15 - m[A21] * Det2_45_05 << 1632 + m[A25]*Det2_45_01;
1728 G4double Det3_245_023 = << 1633 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
1729 m[A20] * Det2_45_23 - m[A22] * Det2_45_03 << 1634 + m[A23]*Det2_45_02;
1730 G4double Det3_245_024 = << 1635 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
1731 m[A20] * Det2_45_24 - m[A22] * Det2_45_04 << 1636 + m[A24]*Det2_45_02;
1732 G4double Det3_245_025 = << 1637 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
1733 m[A20] * Det2_45_25 - m[A22] * Det2_45_05 << 1638 + m[A25]*Det2_45_02;
1734 G4double Det3_245_034 = << 1639 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
1735 m[A20] * Det2_45_34 - m[A23] * Det2_45_04 << 1640 + m[A24]*Det2_45_03;
1736 G4double Det3_245_035 = << 1641 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
1737 m[A20] * Det2_45_35 - m[A23] * Det2_45_05 << 1642 + m[A25]*Det2_45_03;
1738 G4double Det3_245_045 = << 1643 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
1739 m[A20] * Det2_45_45 - m[A24] * Det2_45_05 << 1644 + m[A25]*Det2_45_04;
1740 G4double Det3_245_123 = << 1645 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
1741 m[A21] * Det2_45_23 - m[A22] * Det2_45_13 << 1646 + m[A23]*Det2_45_12;
1742 G4double Det3_245_124 = << 1647 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
1743 m[A21] * Det2_45_24 - m[A22] * Det2_45_14 << 1648 + m[A24]*Det2_45_12;
1744 G4double Det3_245_125 = << 1649 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
1745 m[A21] * Det2_45_25 - m[A22] * Det2_45_15 << 1650 + m[A25]*Det2_45_12;
1746 G4double Det3_245_134 = << 1651 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
1747 m[A21] * Det2_45_34 - m[A23] * Det2_45_14 << 1652 + m[A24]*Det2_45_13;
1748 G4double Det3_245_135 = << 1653 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
1749 m[A21] * Det2_45_35 - m[A23] * Det2_45_15 << 1654 + m[A25]*Det2_45_13;
1750 G4double Det3_245_145 = << 1655 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
1751 m[A21] * Det2_45_45 - m[A24] * Det2_45_15 << 1656 + m[A25]*Det2_45_14;
1752 G4double Det3_245_234 = << 1657 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
1753 m[A22] * Det2_45_34 - m[A23] * Det2_45_24 << 1658 + m[A24]*Det2_45_23;
1754 G4double Det3_245_235 = << 1659 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
1755 m[A22] * Det2_45_35 - m[A23] * Det2_45_25 << 1660 + m[A25]*Det2_45_23;
1756 G4double Det3_245_245 = << 1661 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
1757 m[A22] * Det2_45_45 - m[A24] * Det2_45_25 << 1662 + m[A25]*Det2_45_24;
1758 G4double Det3_345_012 = << 1663 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
1759 m[A30] * Det2_45_12 - m[A31] * Det2_45_02 << 1664 + m[A32]*Det2_45_01;
1760 G4double Det3_345_013 = << 1665 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
1761 m[A30] * Det2_45_13 - m[A31] * Det2_45_03 << 1666 + m[A33]*Det2_45_01;
1762 G4double Det3_345_014 = << 1667 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
1763 m[A30] * Det2_45_14 - m[A31] * Det2_45_04 << 1668 + m[A34]*Det2_45_01;
1764 G4double Det3_345_015 = << 1669 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
1765 m[A30] * Det2_45_15 - m[A31] * Det2_45_05 << 1670 + m[A35]*Det2_45_01;
1766 G4double Det3_345_023 = << 1671 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
1767 m[A30] * Det2_45_23 - m[A32] * Det2_45_03 << 1672 + m[A33]*Det2_45_02;
1768 G4double Det3_345_024 = << 1673 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
1769 m[A30] * Det2_45_24 - m[A32] * Det2_45_04 << 1674 + m[A34]*Det2_45_02;
1770 G4double Det3_345_025 = << 1675 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
1771 m[A30] * Det2_45_25 - m[A32] * Det2_45_05 << 1676 + m[A35]*Det2_45_02;
1772 G4double Det3_345_034 = << 1677 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
1773 m[A30] * Det2_45_34 - m[A33] * Det2_45_04 << 1678 + m[A34]*Det2_45_03;
1774 G4double Det3_345_035 = << 1679 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
1775 m[A30] * Det2_45_35 - m[A33] * Det2_45_05 << 1680 + m[A35]*Det2_45_03;
1776 G4double Det3_345_045 = << 1681 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
1777 m[A30] * Det2_45_45 - m[A34] * Det2_45_05 << 1682 + m[A35]*Det2_45_04;
1778 G4double Det3_345_123 = << 1683 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
1779 m[A31] * Det2_45_23 - m[A32] * Det2_45_13 << 1684 + m[A33]*Det2_45_12;
1780 G4double Det3_345_124 = << 1685 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
1781 m[A31] * Det2_45_24 - m[A32] * Det2_45_14 << 1686 + m[A34]*Det2_45_12;
1782 G4double Det3_345_125 = << 1687 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
1783 m[A31] * Det2_45_25 - m[A32] * Det2_45_15 << 1688 + m[A35]*Det2_45_12;
1784 G4double Det3_345_134 = << 1689 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
1785 m[A31] * Det2_45_34 - m[A33] * Det2_45_14 << 1690 + m[A34]*Det2_45_13;
1786 G4double Det3_345_135 = << 1691 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
1787 m[A31] * Det2_45_35 - m[A33] * Det2_45_15 << 1692 + m[A35]*Det2_45_13;
1788 G4double Det3_345_145 = << 1693 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
1789 m[A31] * Det2_45_45 - m[A34] * Det2_45_15 << 1694 + m[A35]*Det2_45_14;
1790 G4double Det3_345_234 = << 1695 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
1791 m[A32] * Det2_45_34 - m[A33] * Det2_45_24 << 1696 + m[A34]*Det2_45_23;
1792 G4double Det3_345_235 = << 1697 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
1793 m[A32] * Det2_45_35 - m[A33] * Det2_45_25 << 1698 + m[A35]*Det2_45_23;
1794 G4double Det3_345_245 = << 1699 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
1795 m[A32] * Det2_45_45 - m[A34] * Det2_45_25 << 1700 + m[A35]*Det2_45_24;
1796 G4double Det3_345_345 = << 1701 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
1797 m[A33] * Det2_45_45 - m[A34] * Det2_45_35 << 1702 + m[A35]*Det2_45_34;
1798 1703
1799 // Find all NECESSARY 4x4 dets: (55 of the 1704 // Find all NECESSARY 4x4 dets: (55 of them)
1800 1705
1801 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1706 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1802 m[A12] * Det3_234 << 1707 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1803 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1708 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1804 m[A12] * Det3_234 << 1709 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1805 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1710 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1806 m[A13] * Det3_234 << 1711 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1807 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1712 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1808 m[A13] * Det3_234 << 1713 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1809 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1714 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1810 m[A13] * Det3_234 << 1715 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1811 G4double Det4_1235_0123 = m[A10] * Det3_235 << 1716 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1812 m[A12] * Det3_235 << 1717 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1813 G4double Det4_1235_0124 = m[A10] * Det3_235 << 1718 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1814 m[A12] * Det3_235 << 1719 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1815 G4double Det4_1235_0125 = m[A10] * Det3_235 << 1720 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1816 m[A12] * Det3_235 << 1721 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1817 G4double Det4_1235_0134 = m[A10] * Det3_235 << 1722 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1818 m[A13] * Det3_235 << 1723 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1819 G4double Det4_1235_0135 = m[A10] * Det3_235 << 1724 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1820 m[A13] * Det3_235 << 1725 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1821 G4double Det4_1235_0234 = m[A10] * Det3_235 << 1726 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1822 m[A13] * Det3_235 << 1727 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1823 G4double Det4_1235_0235 = m[A10] * Det3_235 << 1728 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1824 m[A13] * Det3_235 << 1729 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1825 G4double Det4_1235_1234 = m[A11] * Det3_235 << 1730 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1826 m[A13] * Det3_235 << 1731 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1827 G4double Det4_1235_1235 = m[A11] * Det3_235 << 1732 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1828 m[A13] * Det3_235 << 1733 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1829 G4double Det4_1245_0123 = m[A10] * Det3_245 << 1734 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1830 m[A12] * Det3_245 << 1735 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1831 G4double Det4_1245_0124 = m[A10] * Det3_245 << 1736 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1832 m[A12] * Det3_245 << 1737 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1833 G4double Det4_1245_0125 = m[A10] * Det3_245 << 1738 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1834 m[A12] * Det3_245 << 1739 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1835 G4double Det4_1245_0134 = m[A10] * Det3_245 << 1740 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1836 m[A13] * Det3_245 << 1741 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1837 G4double Det4_1245_0135 = m[A10] * Det3_245 << 1742 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1838 m[A13] * Det3_245 << 1743 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1839 G4double Det4_1245_0145 = m[A10] * Det3_245 << 1744 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1840 m[A14] * Det3_245 << 1745 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1841 G4double Det4_1245_0234 = m[A10] * Det3_245 << 1746 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1842 m[A13] * Det3_245 << 1747 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1843 G4double Det4_1245_0235 = m[A10] * Det3_245 << 1748 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1844 m[A13] * Det3_245 << 1749 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1845 G4double Det4_1245_0245 = m[A10] * Det3_245 << 1750 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1846 m[A14] * Det3_245 << 1751 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1847 G4double Det4_1245_1234 = m[A11] * Det3_245 << 1752 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1848 m[A13] * Det3_245 << 1753 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1849 G4double Det4_1245_1235 = m[A11] * Det3_245 << 1754 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1850 m[A13] * Det3_245 << 1755 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1851 G4double Det4_1245_1245 = m[A11] * Det3_245 << 1756 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1852 m[A14] * Det3_245 << 1757 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1853 G4double Det4_1345_0123 = m[A10] * Det3_345 << 1758 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1854 m[A12] * Det3_345 << 1759 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1855 G4double Det4_1345_0124 = m[A10] * Det3_345 << 1760 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1856 m[A12] * Det3_345 << 1761 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1857 G4double Det4_1345_0125 = m[A10] * Det3_345 << 1762 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1858 m[A12] * Det3_345 << 1763 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1859 G4double Det4_1345_0134 = m[A10] * Det3_345 << 1764 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1860 m[A13] * Det3_345 << 1765 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1861 G4double Det4_1345_0135 = m[A10] * Det3_345 << 1766 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1862 m[A13] * Det3_345 << 1767 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1863 G4double Det4_1345_0145 = m[A10] * Det3_345 << 1768 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1864 m[A14] * Det3_345 << 1769 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1865 G4double Det4_1345_0234 = m[A10] * Det3_345 << 1770 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1866 m[A13] * Det3_345 << 1771 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1867 G4double Det4_1345_0235 = m[A10] * Det3_345 << 1772 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1868 m[A13] * Det3_345 << 1773 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1869 G4double Det4_1345_0245 = m[A10] * Det3_345 << 1774 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1870 m[A14] * Det3_345 << 1775 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1871 G4double Det4_1345_0345 = m[A10] * Det3_345 << 1776 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1872 m[A14] * Det3_345 << 1777 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1873 G4double Det4_1345_1234 = m[A11] * Det3_345 << 1778 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1874 m[A13] * Det3_345 << 1779 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1875 G4double Det4_1345_1235 = m[A11] * Det3_345 << 1780 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1876 m[A13] * Det3_345 << 1781 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1877 G4double Det4_1345_1245 = m[A11] * Det3_345 << 1782 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1878 m[A14] * Det3_345 << 1783 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1879 G4double Det4_1345_1345 = m[A11] * Det3_345 << 1784 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1880 m[A14] * Det3_345 << 1785 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1881 G4double Det4_2345_0123 = m[A20] * Det3_345 << 1786 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1882 m[A22] * Det3_345 << 1787 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1883 G4double Det4_2345_0124 = m[A20] * Det3_345 << 1788 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1884 m[A22] * Det3_345 << 1789 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1885 G4double Det4_2345_0125 = m[A20] * Det3_345 << 1790 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1886 m[A22] * Det3_345 << 1791 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1887 G4double Det4_2345_0134 = m[A20] * Det3_345 << 1792 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1888 m[A23] * Det3_345 << 1793 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1889 G4double Det4_2345_0135 = m[A20] * Det3_345 << 1794 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1890 m[A23] * Det3_345 << 1795 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1891 G4double Det4_2345_0145 = m[A20] * Det3_345 << 1796 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1892 m[A24] * Det3_345 << 1797 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1893 G4double Det4_2345_0234 = m[A20] * Det3_345 << 1798 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1894 m[A23] * Det3_345 << 1799 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1895 G4double Det4_2345_0235 = m[A20] * Det3_345 << 1800 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1896 m[A23] * Det3_345 << 1801 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1897 G4double Det4_2345_0245 = m[A20] * Det3_345 << 1802 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1898 m[A24] * Det3_345 << 1803 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1899 G4double Det4_2345_0345 = m[A20] * Det3_345 << 1804 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1900 m[A24] * Det3_345 << 1805 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1901 G4double Det4_2345_1234 = m[A21] * Det3_345 << 1806 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1902 m[A23] * Det3_345 << 1807 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1903 G4double Det4_2345_1235 = m[A21] * Det3_345 << 1808 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1904 m[A23] * Det3_345 << 1809 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1905 G4double Det4_2345_1245 = m[A21] * Det3_345 << 1810 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1906 m[A24] * Det3_345 << 1811 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1907 G4double Det4_2345_1345 = m[A21] * Det3_345 << 1812 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1908 m[A24] * Det3_345 << 1813 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1909 G4double Det4_2345_2345 = m[A22] * Det3_345 << 1814 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1910 m[A24] * Det3_345 << 1815 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1911 1816
1912 // Find all NECESSARY 5x5 dets: (19 of the 1817 // Find all NECESSARY 5x5 dets: (19 of them)
1913 1818
1914 G4double Det5_01234_01234 = << 1819 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1915 m[A00] * Det4_1234_1234 - m[A01] * Det4_1 << 1820 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1916 m[A02] * Det4_1234_0134 - m[A03] * Det4_1 << 1821 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1917 G4double Det5_01235_01234 = << 1822 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1918 m[A00] * Det4_1235_1234 - m[A01] * Det4_1 << 1823 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1919 m[A02] * Det4_1235_0134 - m[A03] * Det4_1 << 1824 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1920 G4double Det5_01235_01235 = << 1825 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1921 m[A00] * Det4_1235_1235 - m[A01] * Det4_1 << 1826 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1922 m[A02] * Det4_1235_0135 - m[A03] * Det4_1 << 1827 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1923 G4double Det5_01245_01234 = << 1828 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1924 m[A00] * Det4_1245_1234 - m[A01] * Det4_1 << 1829 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1925 m[A02] * Det4_1245_0134 - m[A03] * Det4_1 << 1830 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1926 G4double Det5_01245_01235 = << 1831 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1927 m[A00] * Det4_1245_1235 - m[A01] * Det4_1 << 1832 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1928 m[A02] * Det4_1245_0135 - m[A03] * Det4_1 << 1833 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1929 G4double Det5_01245_01245 = << 1834 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1930 m[A00] * Det4_1245_1245 - m[A01] * Det4_1 << 1835 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1931 m[A02] * Det4_1245_0145 - m[A04] * Det4_1 << 1836 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1932 G4double Det5_01345_01234 = << 1837 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1933 m[A00] * Det4_1345_1234 - m[A01] * Det4_1 << 1838 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1934 m[A02] * Det4_1345_0134 - m[A03] * Det4_1 << 1839 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1935 G4double Det5_01345_01235 = << 1840 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1936 m[A00] * Det4_1345_1235 - m[A01] * Det4_1 << 1841 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1937 m[A02] * Det4_1345_0135 - m[A03] * Det4_1 << 1842 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1938 G4double Det5_01345_01245 = << 1843 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1939 m[A00] * Det4_1345_1245 - m[A01] * Det4_1 << 1844 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1940 m[A02] * Det4_1345_0145 - m[A04] * Det4_1 << 1845 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1941 G4double Det5_01345_01345 = << 1846 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1942 m[A00] * Det4_1345_1345 - m[A01] * Det4_1 << 1847 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1943 m[A03] * Det4_1345_0145 - m[A04] * Det4_1 << 1848 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1944 G4double Det5_02345_01234 = << 1849 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1945 m[A00] * Det4_2345_1234 - m[A01] * Det4_2 << 1850 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1946 m[A02] * Det4_2345_0134 - m[A03] * Det4_2 << 1851 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1947 G4double Det5_02345_01235 = << 1852 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1948 m[A00] * Det4_2345_1235 - m[A01] * Det4_2 << 1853 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1949 m[A02] * Det4_2345_0135 - m[A03] * Det4_2 << 1854 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1950 G4double Det5_02345_01245 = << 1855 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1951 m[A00] * Det4_2345_1245 - m[A01] * Det4_2 << 1856 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1952 m[A02] * Det4_2345_0145 - m[A04] * Det4_2 << 1857 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1953 G4double Det5_02345_01345 = << 1858 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1954 m[A00] * Det4_2345_1345 - m[A01] * Det4_2 << 1859 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1955 m[A03] * Det4_2345_0145 - m[A04] * Det4_2 << 1860 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1956 G4double Det5_02345_02345 = << 1861
1957 m[A00] * Det4_2345_2345 - m[A02] * Det4_2 << 1862 // Find the determinant
1958 m[A03] * Det4_2345_0245 - m[A04] * Det4_2 << 1863
1959 G4double Det5_12345_01234 = << 1864 G4double det = m[A00]*Det5_12345_12345
1960 m[A10] * Det4_2345_1234 - m[A11] * Det4_2 << 1865 - m[A01]*Det5_12345_02345
1961 m[A12] * Det4_2345_0134 - m[A13] * Det4_2 << 1866 + m[A02]*Det5_12345_01345
1962 G4double Det5_12345_01235 = << 1867 - m[A03]*Det5_12345_01245
1963 m[A10] * Det4_2345_1235 - m[A11] * Det4_2 << 1868 + m[A04]*Det5_12345_01235
1964 m[A12] * Det4_2345_0135 - m[A13] * Det4_2 << 1869 - m[A05]*Det5_12345_01234;
1965 G4double Det5_12345_01245 = <<
1966 m[A10] * Det4_2345_1245 - m[A11] * Det4_2 <<
1967 m[A12] * Det4_2345_0145 - m[A14] * Det4_2 <<
1968 G4double Det5_12345_01345 = <<
1969 m[A10] * Det4_2345_1345 - m[A11] * Det4_2 <<
1970 m[A13] * Det4_2345_0145 - m[A14] * Det4_2 <<
1971 G4double Det5_12345_02345 = <<
1972 m[A10] * Det4_2345_2345 - m[A12] * Det4_2 <<
1973 m[A13] * Det4_2345_0245 - m[A14] * Det4_2 <<
1974 G4double Det5_12345_12345 = <<
1975 m[A11] * Det4_2345_2345 - m[A12] * Det4_2 <<
1976 m[A13] * Det4_2345_1245 - m[A14] * Det4_2 <<
1977 <<
1978 // Find the determinant <<
1979 <<
1980 G4double det = m[A00] * Det5_12345_12345 - <<
1981 m[A02] * Det5_12345_01345 - <<
1982 m[A04] * Det5_12345_01235 - <<
1983 1870
1984 if(det == 0) << 1871 if ( det == 0 )
1985 { << 1872 {
1986 ifail = 1; 1873 ifail = 1;
1987 return; 1874 return;
1988 } << 1875 }
1989 1876
1990 G4double oneOverDet = 1.0 / det; << 1877 G4double oneOverDet = 1.0/det;
1991 G4double mn1OverDet = -oneOverDet; << 1878 G4double mn1OverDet = - oneOverDet;
1992 1879
1993 m[A00] = Det5_12345_12345 * oneOverDet; << 1880 m[A00] = Det5_12345_12345*oneOverDet;
1994 m[A01] = Det5_12345_02345 * mn1OverDet; << 1881 m[A01] = Det5_12345_02345*mn1OverDet;
1995 m[A02] = Det5_12345_01345 * oneOverDet; << 1882 m[A02] = Det5_12345_01345*oneOverDet;
1996 m[A03] = Det5_12345_01245 * mn1OverDet; << 1883 m[A03] = Det5_12345_01245*mn1OverDet;
1997 m[A04] = Det5_12345_01235 * oneOverDet; << 1884 m[A04] = Det5_12345_01235*oneOverDet;
1998 m[A05] = Det5_12345_01234 * mn1OverDet; << 1885 m[A05] = Det5_12345_01234*mn1OverDet;
1999 << 1886
2000 m[A11] = Det5_02345_02345 * oneOverDet; << 1887 m[A11] = Det5_02345_02345*oneOverDet;
2001 m[A12] = Det5_02345_01345 * mn1OverDet; << 1888 m[A12] = Det5_02345_01345*mn1OverDet;
2002 m[A13] = Det5_02345_01245 * oneOverDet; << 1889 m[A13] = Det5_02345_01245*oneOverDet;
2003 m[A14] = Det5_02345_01235 * mn1OverDet; << 1890 m[A14] = Det5_02345_01235*mn1OverDet;
2004 m[A15] = Det5_02345_01234 * oneOverDet; << 1891 m[A15] = Det5_02345_01234*oneOverDet;
2005 << 1892
2006 m[A22] = Det5_01345_01345 * oneOverDet; << 1893 m[A22] = Det5_01345_01345*oneOverDet;
2007 m[A23] = Det5_01345_01245 * mn1OverDet; << 1894 m[A23] = Det5_01345_01245*mn1OverDet;
2008 m[A24] = Det5_01345_01235 * oneOverDet; << 1895 m[A24] = Det5_01345_01235*oneOverDet;
2009 m[A25] = Det5_01345_01234 * mn1OverDet; << 1896 m[A25] = Det5_01345_01234*mn1OverDet;
2010 << 1897
2011 m[A33] = Det5_01245_01245 * oneOverDet; << 1898 m[A33] = Det5_01245_01245*oneOverDet;
2012 m[A34] = Det5_01245_01235 * mn1OverDet; << 1899 m[A34] = Det5_01245_01235*mn1OverDet;
2013 m[A35] = Det5_01245_01234 * oneOverDet; << 1900 m[A35] = Det5_01245_01234*oneOverDet;
2014 1901
2015 m[A44] = Det5_01235_01235 * oneOverDet; << 1902 m[A44] = Det5_01235_01235*oneOverDet;
2016 m[A45] = Det5_01235_01234 * mn1OverDet; << 1903 m[A45] = Det5_01235_01234*mn1OverDet;
2017 1904
2018 m[A55] = Det5_01234_01234 * oneOverDet; << 1905 m[A55] = Det5_01234_01234*oneOverDet;
2019 1906
2020 return; 1907 return;
2021 } 1908 }
2022 1909
2023 void G4ErrorSymMatrix::invertCholesky5(G4int& << 1910 void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
2024 { 1911 {
2025 // Invert by << 1912
>> 1913 // Invert by
2026 // 1914 //
2027 // a) decomposing M = G*G^T with G lower tr 1915 // a) decomposing M = G*G^T with G lower triangular
2028 // (if M is not positive definite this w 1916 // (if M is not positive definite this will fail, leaving this unchanged)
2029 // b) inverting G to form H 1917 // b) inverting G to form H
2030 // c) multiplying H^T * H to get M^-1. 1918 // c) multiplying H^T * H to get M^-1.
2031 // 1919 //
2032 // If the matrix is pos. def. it is inverte 1920 // If the matrix is pos. def. it is inverted and 1 is returned.
2033 // If the matrix is not pos. def. it remain 1921 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2034 1922
2035 G4double h10; // below-diagonal elements o << 1923 G4double h10; // below-diagonal elements of H
2036 G4double h20, h21; 1924 G4double h20, h21;
2037 G4double h30, h31, h32; 1925 G4double h30, h31, h32;
2038 G4double h40, h41, h42, h43; 1926 G4double h40, h41, h42, h43;
2039 1927
2040 G4double h00, h11, h22, h33, h44; // 1/dia << 1928 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
2041 // diago << 1929 // diagonal elements of H
2042 1930
2043 G4double g10; // below-diagonal elements o << 1931 G4double g10; // below-diagonal elements of G
2044 G4double g20, g21; 1932 G4double g20, g21;
2045 G4double g30, g31, g32; 1933 G4double g30, g31, g32;
2046 G4double g40, g41, g42, g43; 1934 G4double g40, g41, g42, g43;
2047 1935
2048 ifail = 1; // We start by assuing we won't 1936 ifail = 1; // We start by assuing we won't succeed...
2049 1937
2050 // Form G -- compute diagonal members of H 1938 // Form G -- compute diagonal members of H directly rather than of G
2051 //------- 1939 //-------
2052 1940
2053 // Scale first column by 1/sqrt(A00) 1941 // Scale first column by 1/sqrt(A00)
2054 1942
2055 h00 = m[A00]; << 1943 h00 = m[A00];
2056 if(h00 <= 0) << 1944 if (h00 <= 0) { return; }
2057 { <<
2058 return; <<
2059 } <<
2060 h00 = 1.0 / std::sqrt(h00); 1945 h00 = 1.0 / std::sqrt(h00);
2061 1946
2062 g10 = m[A10] * h00; 1947 g10 = m[A10] * h00;
2063 g20 = m[A20] * h00; 1948 g20 = m[A20] * h00;
2064 g30 = m[A30] * h00; 1949 g30 = m[A30] * h00;
2065 g40 = m[A40] * h00; 1950 g40 = m[A40] * h00;
2066 1951
2067 // Form G11 (actually, just h11) 1952 // Form G11 (actually, just h11)
2068 1953
2069 h11 = m[A11] - (g10 * g10); 1954 h11 = m[A11] - (g10 * g10);
2070 if(h11 <= 0) << 1955 if (h11 <= 0) { return; }
2071 { <<
2072 return; <<
2073 } <<
2074 h11 = 1.0 / std::sqrt(h11); 1956 h11 = 1.0 / std::sqrt(h11);
2075 1957
2076 // Subtract inter-column column dot product 1958 // Subtract inter-column column dot products from rest of column 1 and
2077 // scale to get column 1 of G << 1959 // scale to get column 1 of G
2078 1960
2079 g21 = (m[A21] - (g10 * g20)) * h11; 1961 g21 = (m[A21] - (g10 * g20)) * h11;
2080 g31 = (m[A31] - (g10 * g30)) * h11; 1962 g31 = (m[A31] - (g10 * g30)) * h11;
2081 g41 = (m[A41] - (g10 * g40)) * h11; 1963 g41 = (m[A41] - (g10 * g40)) * h11;
2082 1964
2083 // Form G22 (actually, just h22) 1965 // Form G22 (actually, just h22)
2084 1966
2085 h22 = m[A22] - (g20 * g20) - (g21 * g21); 1967 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2086 if(h22 <= 0) << 1968 if (h22 <= 0) { return; }
2087 { <<
2088 return; <<
2089 } <<
2090 h22 = 1.0 / std::sqrt(h22); 1969 h22 = 1.0 / std::sqrt(h22);
2091 1970
2092 // Subtract inter-column column dot product 1971 // Subtract inter-column column dot products from rest of column 2 and
2093 // scale to get column 2 of G << 1972 // scale to get column 2 of G
2094 1973
2095 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 1974 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2096 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 1975 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2097 1976
2098 // Form G33 (actually, just h33) 1977 // Form G33 (actually, just h33)
2099 1978
2100 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 1979 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2101 if(h33 <= 0) << 1980 if (h33 <= 0) { return; }
2102 { <<
2103 return; <<
2104 } <<
2105 h33 = 1.0 / std::sqrt(h33); 1981 h33 = 1.0 / std::sqrt(h33);
2106 1982
2107 // Subtract inter-column column dot product 1983 // Subtract inter-column column dot product from A43 and scale to get G43
2108 1984
2109 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 1985 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2110 1986
2111 // Finally form h44 - if this is possible i 1987 // Finally form h44 - if this is possible inversion succeeds
2112 1988
2113 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 1989 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2114 if(h44 <= 0) << 1990 if (h44 <= 0) { return; }
2115 { <<
2116 return; <<
2117 } <<
2118 h44 = 1.0 / std::sqrt(h44); 1991 h44 = 1.0 / std::sqrt(h44);
2119 1992
2120 // Form H = 1/G -- diagonal members of H ar 1993 // Form H = 1/G -- diagonal members of H are already correct
2121 //------------- 1994 //-------------
2122 1995
2123 // The order here is dictated by speed cons 1996 // The order here is dictated by speed considerations
2124 1997
2125 h43 = -h33 * g43 * h44; << 1998 h43 = -h33 * g43 * h44;
2126 h32 = -h22 * g32 * h33; << 1999 h32 = -h22 * g32 * h33;
2127 h42 = -h22 * (g32 * h43 + g42 * h44); 2000 h42 = -h22 * (g32 * h43 + g42 * h44);
2128 h21 = -h11 * g21 * h22; << 2001 h21 = -h11 * g21 * h22;
2129 h31 = -h11 * (g21 * h32 + g31 * h33); 2002 h31 = -h11 * (g21 * h32 + g31 * h33);
2130 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2003 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2131 h10 = -h00 * g10 * h11; << 2004 h10 = -h00 * g10 * h11;
2132 h20 = -h00 * (g10 * h21 + g20 * h22); 2005 h20 = -h00 * (g10 * h21 + g20 * h22);
2133 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2006 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2134 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2007 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2135 2008
2136 // Change this to its inverse = H^T*H 2009 // Change this to its inverse = H^T*H
2137 //------------------------------------ 2010 //------------------------------------
2138 2011
2139 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 2012 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2140 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 2013 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2141 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2014 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
2142 m[A02] = h20 * h22 + h30 * h32 + h40 * h42; 2015 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2143 m[A12] = h21 * h22 + h31 * h32 + h41 * h42; 2016 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2144 m[A22] = h22 * h22 + h32 * h32 + h42 * h42; 2017 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2145 m[A03] = h30 * h33 + h40 * h43; 2018 m[A03] = h30 * h33 + h40 * h43;
2146 m[A13] = h31 * h33 + h41 * h43; 2019 m[A13] = h31 * h33 + h41 * h43;
2147 m[A23] = h32 * h33 + h42 * h43; 2020 m[A23] = h32 * h33 + h42 * h43;
2148 m[A33] = h33 * h33 + h43 * h43; 2021 m[A33] = h33 * h33 + h43 * h43;
2149 m[A04] = h40 * h44; 2022 m[A04] = h40 * h44;
2150 m[A14] = h41 * h44; 2023 m[A14] = h41 * h44;
2151 m[A24] = h42 * h44; 2024 m[A24] = h42 * h44;
2152 m[A34] = h43 * h44; 2025 m[A34] = h43 * h44;
2153 m[A44] = h44 * h44; 2026 m[A44] = h44 * h44;
2154 2027
2155 ifail = 0; 2028 ifail = 0;
2156 return; 2029 return;
2157 } 2030 }
2158 2031
2159 void G4ErrorSymMatrix::invertCholesky6(G4int& << 2032 void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
2160 { 2033 {
2161 // Invert by << 2034 // Invert by
2162 // 2035 //
2163 // a) decomposing M = G*G^T with G lower tr 2036 // a) decomposing M = G*G^T with G lower triangular
2164 // (if M is not positive definite this w 2037 // (if M is not positive definite this will fail, leaving this unchanged)
2165 // b) inverting G to form H 2038 // b) inverting G to form H
2166 // c) multiplying H^T * H to get M^-1. 2039 // c) multiplying H^T * H to get M^-1.
2167 // 2040 //
2168 // If the matrix is pos. def. it is inverte 2041 // If the matrix is pos. def. it is inverted and 1 is returned.
2169 // If the matrix is not pos. def. it remain 2042 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2170 2043
2171 G4double h10; // below-diagonal elements o << 2044 G4double h10; // below-diagonal elements of H
2172 G4double h20, h21; 2045 G4double h20, h21;
2173 G4double h30, h31, h32; 2046 G4double h30, h31, h32;
2174 G4double h40, h41, h42, h43; 2047 G4double h40, h41, h42, h43;
2175 G4double h50, h51, h52, h53, h54; 2048 G4double h50, h51, h52, h53, h54;
2176 2049
2177 G4double h00, h11, h22, h33, h44, h55; // << 2050 G4double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G =
2178 // 2051 // diagonal elements of H
2179 2052
2180 G4double g10; // below-diagonal elements o << 2053 G4double g10; // below-diagonal elements of G
2181 G4double g20, g21; 2054 G4double g20, g21;
2182 G4double g30, g31, g32; 2055 G4double g30, g31, g32;
2183 G4double g40, g41, g42, g43; 2056 G4double g40, g41, g42, g43;
2184 G4double g50, g51, g52, g53, g54; 2057 G4double g50, g51, g52, g53, g54;
2185 2058
2186 ifail = 1; // We start by assuing we won't 2059 ifail = 1; // We start by assuing we won't succeed...
2187 2060
2188 // Form G -- compute diagonal members of H 2061 // Form G -- compute diagonal members of H directly rather than of G
2189 //------- 2062 //-------
2190 2063
2191 // Scale first column by 1/sqrt(A00) 2064 // Scale first column by 1/sqrt(A00)
2192 2065
2193 h00 = m[A00]; << 2066 h00 = m[A00];
2194 if(h00 <= 0) << 2067 if (h00 <= 0) { return; }
2195 { <<
2196 return; <<
2197 } <<
2198 h00 = 1.0 / std::sqrt(h00); 2068 h00 = 1.0 / std::sqrt(h00);
2199 2069
2200 g10 = m[A10] * h00; 2070 g10 = m[A10] * h00;
2201 g20 = m[A20] * h00; 2071 g20 = m[A20] * h00;
2202 g30 = m[A30] * h00; 2072 g30 = m[A30] * h00;
2203 g40 = m[A40] * h00; 2073 g40 = m[A40] * h00;
2204 g50 = m[A50] * h00; 2074 g50 = m[A50] * h00;
2205 2075
2206 // Form G11 (actually, just h11) 2076 // Form G11 (actually, just h11)
2207 2077
2208 h11 = m[A11] - (g10 * g10); 2078 h11 = m[A11] - (g10 * g10);
2209 if(h11 <= 0) << 2079 if (h11 <= 0) { return; }
2210 { <<
2211 return; <<
2212 } <<
2213 h11 = 1.0 / std::sqrt(h11); 2080 h11 = 1.0 / std::sqrt(h11);
2214 2081
2215 // Subtract inter-column column dot product 2082 // Subtract inter-column column dot products from rest of column 1 and
2216 // scale to get column 1 of G << 2083 // scale to get column 1 of G
2217 2084
2218 g21 = (m[A21] - (g10 * g20)) * h11; 2085 g21 = (m[A21] - (g10 * g20)) * h11;
2219 g31 = (m[A31] - (g10 * g30)) * h11; 2086 g31 = (m[A31] - (g10 * g30)) * h11;
2220 g41 = (m[A41] - (g10 * g40)) * h11; 2087 g41 = (m[A41] - (g10 * g40)) * h11;
2221 g51 = (m[A51] - (g10 * g50)) * h11; 2088 g51 = (m[A51] - (g10 * g50)) * h11;
2222 2089
2223 // Form G22 (actually, just h22) 2090 // Form G22 (actually, just h22)
2224 2091
2225 h22 = m[A22] - (g20 * g20) - (g21 * g21); 2092 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2226 if(h22 <= 0) << 2093 if (h22 <= 0) { return; }
2227 { <<
2228 return; <<
2229 } <<
2230 h22 = 1.0 / std::sqrt(h22); 2094 h22 = 1.0 / std::sqrt(h22);
2231 2095
2232 // Subtract inter-column column dot product 2096 // Subtract inter-column column dot products from rest of column 2 and
2233 // scale to get column 2 of G << 2097 // scale to get column 2 of G
2234 2098
2235 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 2099 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2236 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 2100 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2237 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) 2101 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2238 2102
2239 // Form G33 (actually, just h33) 2103 // Form G33 (actually, just h33)
2240 2104
2241 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 2105 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2242 if(h33 <= 0) << 2106 if (h33 <= 0) { return; }
2243 { <<
2244 return; <<
2245 } <<
2246 h33 = 1.0 / std::sqrt(h33); 2107 h33 = 1.0 / std::sqrt(h33);
2247 2108
2248 // Subtract inter-column column dot product 2109 // Subtract inter-column column dot products from rest of column 3 and
2249 // scale to get column 3 of G << 2110 // scale to get column 3 of G
2250 2111
2251 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 2112 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2252 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - 2113 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2253 2114
2254 // Form G44 (actually, just h44) 2115 // Form G44 (actually, just h44)
2255 2116
2256 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 2117 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2257 if(h44 <= 0) << 2118 if (h44 <= 0) { return; }
2258 { <<
2259 return; <<
2260 } <<
2261 h44 = 1.0 / std::sqrt(h44); 2119 h44 = 1.0 / std::sqrt(h44);
2262 2120
2263 // Subtract inter-column column dot product 2121 // Subtract inter-column column dot product from M54 and scale to get G54
2264 2122
2265 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - 2123 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2266 2124
2267 // Finally form h55 - if this is possible i 2125 // Finally form h55 - if this is possible inversion succeeds
2268 2126
2269 h55 = m[A55] - (g50 * g50) - (g51 * g51) - << 2127 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
2270 (g54 * g54); << 2128 if (h55 <= 0) { return; }
2271 if(h55 <= 0) <<
2272 { <<
2273 return; <<
2274 } <<
2275 h55 = 1.0 / std::sqrt(h55); 2129 h55 = 1.0 / std::sqrt(h55);
2276 2130
2277 // Form H = 1/G -- diagonal members of H ar 2131 // Form H = 1/G -- diagonal members of H are already correct
2278 //------------- 2132 //-------------
2279 2133
2280 // The order here is dictated by speed cons 2134 // The order here is dictated by speed considerations
2281 2135
2282 h54 = -h44 * g54 * h55; << 2136 h54 = -h44 * g54 * h55;
2283 h43 = -h33 * g43 * h44; << 2137 h43 = -h33 * g43 * h44;
2284 h53 = -h33 * (g43 * h54 + g53 * h55); 2138 h53 = -h33 * (g43 * h54 + g53 * h55);
2285 h32 = -h22 * g32 * h33; << 2139 h32 = -h22 * g32 * h33;
2286 h42 = -h22 * (g32 * h43 + g42 * h44); 2140 h42 = -h22 * (g32 * h43 + g42 * h44);
2287 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * 2141 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
2288 h21 = -h11 * g21 * h22; << 2142 h21 = -h11 * g21 * h22;
2289 h31 = -h11 * (g21 * h32 + g31 * h33); 2143 h31 = -h11 * (g21 * h32 + g31 * h33);
2290 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2144 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2291 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * 2145 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2292 h10 = -h00 * g10 * h11; << 2146 h10 = -h00 * g10 * h11;
2293 h20 = -h00 * (g10 * h21 + g20 * h22); 2147 h20 = -h00 * (g10 * h21 + g20 * h22);
2294 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2148 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2295 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2149 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2296 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * 2150 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2297 2151
2298 // Change this to its inverse = H^T*H 2152 // Change this to its inverse = H^T*H
2299 //------------------------------------ 2153 //------------------------------------
2300 2154
2301 m[A00] = << 2155 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
2302 h00 * h00 + h10 * h10 + h20 * h20 + h30 * <<
2303 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 2156 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2304 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2157 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2305 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 2158 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2306 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 2159 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2307 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 2160 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
2308 m[A03] = h30 * h33 + h40 * h43 + h50 * h53; 2161 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2309 m[A13] = h31 * h33 + h41 * h43 + h51 * h53; 2162 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2310 m[A23] = h32 * h33 + h42 * h43 + h52 * h53; 2163 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2311 m[A33] = h33 * h33 + h43 * h43 + h53 * h53; 2164 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2312 m[A04] = h40 * h44 + h50 * h54; 2165 m[A04] = h40 * h44 + h50 * h54;
2313 m[A14] = h41 * h44 + h51 * h54; 2166 m[A14] = h41 * h44 + h51 * h54;
2314 m[A24] = h42 * h44 + h52 * h54; 2167 m[A24] = h42 * h44 + h52 * h54;
2315 m[A34] = h43 * h44 + h53 * h54; 2168 m[A34] = h43 * h44 + h53 * h54;
2316 m[A44] = h44 * h44 + h54 * h54; 2169 m[A44] = h44 * h44 + h54 * h54;
2317 m[A05] = h50 * h55; 2170 m[A05] = h50 * h55;
2318 m[A15] = h51 * h55; 2171 m[A15] = h51 * h55;
2319 m[A25] = h52 * h55; 2172 m[A25] = h52 * h55;
2320 m[A35] = h53 * h55; 2173 m[A35] = h53 * h55;
2321 m[A45] = h54 * h55; 2174 m[A45] = h54 * h55;
2322 m[A55] = h55 * h55; 2175 m[A55] = h55 * h55;
2323 2176
2324 ifail = 0; 2177 ifail = 0;
2325 return; 2178 return;
2326 } 2179 }
2327 2180
2328 void G4ErrorSymMatrix::invert4(G4int& ifail) << 2181 void G4ErrorSymMatrix::invert4 (G4int & ifail)
2329 { 2182 {
2330 ifail = 0; 2183 ifail = 0;
2331 2184
2332 // Find all NECESSARY 2x2 dets: (14 of the 2185 // Find all NECESSARY 2x2 dets: (14 of them)
2333 2186
2334 G4double Det2_12_01 = m[A10] * m[A21] - m[A << 2187 G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
2335 G4double Det2_12_02 = m[A10] * m[A22] - m[A << 2188 G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
2336 G4double Det2_12_12 = m[A11] * m[A22] - m[A << 2189 G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
2337 G4double Det2_13_01 = m[A10] * m[A31] - m[A << 2190 G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
2338 G4double Det2_13_02 = m[A10] * m[A32] - m[A << 2191 G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
2339 G4double Det2_13_03 = m[A10] * m[A33] - m[A << 2192 G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
2340 G4double Det2_13_12 = m[A11] * m[A32] - m[A << 2193 G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
2341 G4double Det2_13_13 = m[A11] * m[A33] - m[A << 2194 G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
2342 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 2195 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
2343 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 2196 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
2344 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 2197 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
2345 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 2198 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
2346 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 2199 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
2347 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 2200 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
2348 2201
2349 // Find all NECESSARY 3x3 dets: (10 of th 2202 // Find all NECESSARY 3x3 dets: (10 of them)
2350 2203
2351 G4double Det3_012_012 = << 2204 G4double Det3_012_012 = m[A00]*Det2_12_12 - m[A01]*Det2_12_02
2352 m[A00] * Det2_12_12 - m[A01] * Det2_12_02 << 2205 + m[A02]*Det2_12_01;
2353 G4double Det3_013_012 = << 2206 G4double Det3_013_012 = m[A00]*Det2_13_12 - m[A01]*Det2_13_02
2354 m[A00] * Det2_13_12 - m[A01] * Det2_13_02 << 2207 + m[A02]*Det2_13_01;
2355 G4double Det3_013_013 = << 2208 G4double Det3_013_013 = m[A00]*Det2_13_13 - m[A01]*Det2_13_03
2356 m[A00] * Det2_13_13 - m[A01] * Det2_13_03 << 2209 + m[A03]*Det2_13_01;
2357 G4double Det3_023_012 = << 2210 G4double Det3_023_012 = m[A00]*Det2_23_12 - m[A01]*Det2_23_02
2358 m[A00] * Det2_23_12 - m[A01] * Det2_23_02 << 2211 + m[A02]*Det2_23_01;
2359 G4double Det3_023_013 = << 2212 G4double Det3_023_013 = m[A00]*Det2_23_13 - m[A01]*Det2_23_03
2360 m[A00] * Det2_23_13 - m[A01] * Det2_23_03 << 2213 + m[A03]*Det2_23_01;
2361 G4double Det3_023_023 = << 2214 G4double Det3_023_023 = m[A00]*Det2_23_23 - m[A02]*Det2_23_03
2362 m[A00] * Det2_23_23 - m[A02] * Det2_23_03 << 2215 + m[A03]*Det2_23_02;
2363 G4double Det3_123_012 = << 2216 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
2364 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 2217 + m[A12]*Det2_23_01;
2365 G4double Det3_123_013 = << 2218 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
2366 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 2219 + m[A13]*Det2_23_01;
2367 G4double Det3_123_023 = << 2220 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
2368 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 2221 + m[A13]*Det2_23_02;
2369 G4double Det3_123_123 = << 2222 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
2370 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 2223 + m[A13]*Det2_23_12;
2371 2224
2372 // Find the 4x4 det: 2225 // Find the 4x4 det:
2373 2226
2374 G4double det = m[A00] * Det3_123_123 - m[A0 << 2227 G4double det = m[A00]*Det3_123_123
2375 m[A02] * Det3_123_013 - m[A0 << 2228 - m[A01]*Det3_123_023
>> 2229 + m[A02]*Det3_123_013
>> 2230 - m[A03]*Det3_123_012;
2376 2231
2377 if(det == 0) << 2232 if ( det == 0 )
2378 { << 2233 {
2379 ifail = 1; 2234 ifail = 1;
2380 return; 2235 return;
2381 } << 2236 }
>> 2237
>> 2238 G4double oneOverDet = 1.0/det;
>> 2239 G4double mn1OverDet = - oneOverDet;
2382 2240
2383 G4double oneOverDet = 1.0 / det; << 2241 m[A00] = Det3_123_123 * oneOverDet;
2384 G4double mn1OverDet = -oneOverDet; << 2242 m[A01] = Det3_123_023 * mn1OverDet;
>> 2243 m[A02] = Det3_123_013 * oneOverDet;
>> 2244 m[A03] = Det3_123_012 * mn1OverDet;
2385 2245
2386 m[A00] = Det3_123_123 * oneOverDet; <<
2387 m[A01] = Det3_123_023 * mn1OverDet; <<
2388 m[A02] = Det3_123_013 * oneOverDet; <<
2389 m[A03] = Det3_123_012 * mn1OverDet; <<
2390 <<
2391 m[A11] = Det3_023_023 * oneOverDet; <<
2392 m[A12] = Det3_023_013 * mn1OverDet; <<
2393 m[A13] = Det3_023_012 * oneOverDet; <<
2394 2246
2395 m[A22] = Det3_013_013 * oneOverDet; << 2247 m[A11] = Det3_023_023 * oneOverDet;
2396 m[A23] = Det3_013_012 * mn1OverDet; << 2248 m[A12] = Det3_023_013 * mn1OverDet;
>> 2249 m[A13] = Det3_023_012 * oneOverDet;
2397 2250
2398 m[A33] = Det3_012_012 * oneOverDet; << 2251 m[A22] = Det3_013_013 * oneOverDet;
>> 2252 m[A23] = Det3_013_012 * mn1OverDet;
>> 2253
>> 2254 m[A33] = Det3_012_012 * oneOverDet;
2399 2255
2400 return; 2256 return;
2401 } 2257 }
2402 2258
2403 void G4ErrorSymMatrix::invertHaywood4(G4int& << 2259 void G4ErrorSymMatrix::invertHaywood4 (G4int & ifail)
2404 { 2260 {
2405 invert4(ifail); // For the 4x4 case, the m << 2261 invert4(ifail); // For the 4x4 case, the method we use for invert is already
2406 // the Haywood method. << 2262 // the Haywood method.
2407 } 2263 }
>> 2264
2408 2265