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