| 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 69766 2013-05-14 14:33:55Z 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)
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)
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)
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)
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 { <<
516 if(mat1.nrow * mat1.nrow != size) <<
517 { <<
518 size = mat1.nrow * mat1.nrow; <<
519 m.resize(size); <<
520 } <<
521 nrow = mat1.nrow; <<
522 ncol = mat1.nrow; <<
523 G4int n = ncol; <<
524 G4ErrorMatrixConstIter sjk = mat1.m.begin(); <<
525 G4ErrorMatrixIter m1j = m.begin(); <<
526 G4ErrorMatrixIter mj = m.begin(); <<
527 // j >= k <<
528 for(G4int j = 1; j <= num_row(); j++) <<
529 { <<
530 G4ErrorMatrixIter mjk = mj; <<
531 G4ErrorMatrixIter mkj = m1j; <<
532 for(G4int k = 1; k <= j; k++) <<
533 { <<
534 *(mjk++) = *sjk; <<
535 if(j != k) <<
536 *mkj = *sjk; <<
537 sjk++; <<
538 mkj += n; <<
539 } <<
540 mj += n; <<
541 m1j++; <<
542 } <<
543 return (*this); <<
544 } <<
545 <<
546 G4ErrorSymMatrix& G4ErrorSymMatrix::operator=( <<
547 { 479 {
548 if(&mat1 == this) << 480 if(mat1.nrow*mat1.nrow != size)
549 { << 481 {
550 return *this; << 482 size = mat1.nrow * mat1.nrow;
551 } << 483 m.resize(size);
552 if(mat1.nrow != nrow) << 484 }
553 { << 485 nrow = mat1.nrow;
554 nrow = mat1.nrow; << 486 ncol = mat1.nrow;
555 size = mat1.size; << 487 G4int n = ncol;
556 m.resize(size); << 488 G4ErrorMatrixConstIter sjk = mat1.m.begin();
557 } << 489 G4ErrorMatrixIter m1j = m.begin();
558 m = mat1.m; << 490 G4ErrorMatrixIter mj = m.begin();
559 return (*this); << 491 // j >= k
>> 492 for(G4int j=1;j<=num_row();j++)
>> 493 {
>> 494 G4ErrorMatrixIter mjk = mj;
>> 495 G4ErrorMatrixIter mkj = m1j;
>> 496 for(G4int k=1;k<=j;k++)
>> 497 {
>> 498 *(mjk++) = *sjk;
>> 499 if(j!=k) *mkj = *sjk;
>> 500 sjk++;
>> 501 mkj += n;
>> 502 }
>> 503 mj += n;
>> 504 m1j++;
>> 505 }
>> 506 return (*this);
>> 507 }
>> 508
>> 509 G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &mat1)
>> 510 {
>> 511 if (&mat1 == this) { return *this; }
>> 512 if(mat1.nrow != nrow)
>> 513 {
>> 514 nrow = mat1.nrow;
>> 515 size = mat1.size;
>> 516 m.resize(size);
>> 517 }
>> 518 m = mat1.m;
>> 519 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 std::vector<G4double> xvec (max_array);
895 if(!xvec) << 848 static std::vector<G4int> pivv (max_array);
896 xvec = new std::vector<G4double>(max_array << 849 typedef std::vector<G4int>::iterator pivIter;
897 static G4ThreadLocal std::vector<G4int>* piv << 850 if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
898 if(!pivv) << 851 if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
899 pivv = new std::vector<G4int>(max_array); << 852 // Note - resize should do nothing if the size is already larger than nrow,
900 typedef std::vector<G4int>::iterator pivIter << 853 // but on VC++ there are indications that it does so we check.
901 if(xvec->size() < static_cast<unsigned int>( << 854 // Note - the data elements in a vector are guaranteed to be contiguous,
902 xvec->resize(nrow); << 855 // so x[i] and piv[i] are optimally fast.
903 if(pivv->size() < static_cast<unsigned int>( << 856 G4ErrorMatrixIter x = xvec.begin();
904 pivv->resize(nrow); << 857 // x[i] is used as helper storage, needs to have at least size nrow.
905 // Note - resize should do nothing if the s << 858 pivIter piv = pivv.begin();
906 // but on VC++ there are indications << 859 // piv[i] is used to store details of exchanges
907 // Note - the data elements in a vector are << 860
908 // so x[i] and piv[i] are optimally f <<
909 G4ErrorMatrixIter x = xvec->begin(); <<
910 // x[i] is used as helper storage, needs to <<
911 pivIter piv = pivv->begin(); <<
912 // piv[i] is used to store details of exchan <<
913 <<
914 G4double temp1, temp2; 861 G4double temp1, temp2;
915 G4ErrorMatrixIter ip, mjj, iq; 862 G4ErrorMatrixIter ip, mjj, iq;
916 G4double lambda, sigma; 863 G4double lambda, sigma;
917 const G4double alpha = .6404; // = (1+sqr << 864 const G4double alpha = .6404; // = (1+sqrt(17))/8
918 const G4double epsilon = 32 * DBL_EPSILON; << 865 const G4double epsilon = 32*DBL_EPSILON;
919 // whenever a sum of two doubles is below or << 866 // whenever a sum of two doubles is below or equal to epsilon
920 // it is set to zero. << 867 // it is set to zero.
921 // this constant could be set to zero but th << 868 // this constant could be set to zero but then the algorithm
922 // doesn't neccessarily detect that a matrix << 869 // doesn't neccessarily detect that a matrix is singular
923 << 870
924 for(i = 0; i < nrow; i++) << 871 for (i = 0; i < nrow; i++)
925 { 872 {
926 piv[i] = i + 1; << 873 piv[i] = i+1;
927 } 874 }
928 << 875
929 ifail = 0; 876 ifail = 0;
930 << 877
931 // compute the factorization P*A*P^T = L * D << 878 // compute the factorization P*A*P^T = L * D * L^T
932 // L is unit lower triangular, D is direct s 879 // 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 880 // L and D^-1 are stored in A = *this, P is stored in piv[]
934 << 881
935 for(j = 1; j < nrow; j += ss) // main loop << 882 for (j=1; j < nrow; j+=ss) // main loop over columns
936 { 883 {
937 mjj = m.begin() + j * (j - 1) / 2 + j - << 884 mjj = m.begin() + j*(j-1)/2 + j-1;
938 lambda = 0; // compute lambda = max of A( << 885 lambda = 0; // compute lambda = max of A(j+1:n,j)
939 pivrow = j + 1; << 886 pivrow = j+1;
940 ip = m.begin() + (j + 1) * j / 2 + j - << 887 ip = m.begin() + (j+1)*j/2 + j-1;
941 for(i = j + 1; i <= nrow; ip += i++) << 888 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 { 889 {
1008 *ip -= temp1 * *(m.begin() + k * << 890 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 { 891 {
1150 *ip = 0; << 892 lambda = std::fabs(*ip);
>> 893 pivrow = i;
1151 } 894 }
1152 } <<
1153 } 895 }
1154 // update L << 896 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 { 897 {
1161 temp1 = 0; << 898 if (*mjj == 0)
>> 899 {
>> 900 ifail = 1;
>> 901 return;
>> 902 }
>> 903 ss=1;
>> 904 *mjj = 1./ *mjj;
1162 } 905 }
1163 *(ip + 1) = *ip * *(mjj + j) + *( << 906 else
1164 if(std::fabs(*(ip + 1)) <= epsilo <<
1165 { 907 {
1166 *(ip + 1) = 0; << 908 if (std::fabs(*mjj) >= lambda*alpha)
>> 909 {
>> 910 ss=1;
>> 911 pivrow=j;
>> 912 }
>> 913 else
>> 914 {
>> 915 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
>> 916 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
>> 917 for (k=j; k < pivrow; k++)
>> 918 {
>> 919 if (std::fabs(*ip) > sigma)
>> 920 sigma = std::fabs(*ip);
>> 921 ip++;
>> 922 }
>> 923 if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
>> 924 {
>> 925 ss=1;
>> 926 pivrow = j;
>> 927 }
>> 928 else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
>> 929 >= alpha * sigma)
>> 930 { ss=1; }
>> 931 else
>> 932 { ss=2; }
>> 933 }
>> 934 if (pivrow == j) // no permutation neccessary
>> 935 {
>> 936 piv[j-1] = pivrow;
>> 937 if (*mjj == 0)
>> 938 {
>> 939 ifail=1;
>> 940 return;
>> 941 }
>> 942 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 943
>> 944 // update A(j+1:n, j+1,n)
>> 945 for (i=j+1; i <= nrow; i++)
>> 946 {
>> 947 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 948 ip = m.begin()+i*(i-1)/2+j;
>> 949 for (k=j+1; k<=i; k++)
>> 950 {
>> 951 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 952 if (std::fabs(*ip) <= epsilon)
>> 953 { *ip=0; }
>> 954 ip++;
>> 955 }
>> 956 }
>> 957 // update L
>> 958 ip = m.begin() + (j+1)*j/2 + j-1;
>> 959 for (i=j+1; i <= nrow; ip += i++)
>> 960 {
>> 961 *ip *= temp2;
>> 962 }
>> 963 }
>> 964 else if (ss==1) // 1x1 pivot
>> 965 {
>> 966 piv[j-1] = pivrow;
>> 967
>> 968 // interchange rows and columns j and pivrow in
>> 969 // submatrix (j:n,j:n)
>> 970 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
>> 971 for (i=j+1; i < pivrow; i++, ip++)
>> 972 {
>> 973 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
>> 974 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
>> 975 *ip = temp1;
>> 976 }
>> 977 temp1 = *mjj;
>> 978 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
>> 979 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
>> 980 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
>> 981 iq = ip + pivrow-j;
>> 982 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 983 {
>> 984 temp1 = *iq;
>> 985 *iq = *ip;
>> 986 *ip = temp1;
>> 987 }
>> 988
>> 989 if (*mjj == 0)
>> 990 {
>> 991 ifail = 1;
>> 992 return;
>> 993 }
>> 994 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 995
>> 996 // update A(j+1:n, j+1:n)
>> 997 for (i = j+1; i <= nrow; i++)
>> 998 {
>> 999 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 1000 ip = m.begin()+i*(i-1)/2+j;
>> 1001 for (k=j+1; k<=i; k++)
>> 1002 {
>> 1003 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 1004 if (std::fabs(*ip) <= epsilon)
>> 1005 { *ip=0; }
>> 1006 ip++;
>> 1007 }
>> 1008 }
>> 1009 // update L
>> 1010 ip = m.begin() + (j+1)*j/2 + j-1;
>> 1011 for (i=j+1; i<=nrow; ip += i++)
>> 1012 {
>> 1013 *ip *= temp2;
>> 1014 }
>> 1015 }
>> 1016 else // ss=2, ie use a 2x2 pivot
>> 1017 {
>> 1018 piv[j-1] = -pivrow;
>> 1019 piv[j] = 0; // that means this is the second row of a 2x2 pivot
>> 1020
>> 1021 if (j+1 != pivrow)
>> 1022 {
>> 1023 // interchange rows and columns j+1 and pivrow in
>> 1024 // submatrix (j:n,j:n)
>> 1025 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
>> 1026 for (i=j+2; i < pivrow; i++, ip++)
>> 1027 {
>> 1028 temp1 = *(m.begin() + i*(i-1)/2 + j);
>> 1029 *(m.begin() + i*(i-1)/2 + j) = *ip;
>> 1030 *ip = temp1;
>> 1031 }
>> 1032 temp1 = *(mjj + j + 1);
>> 1033 *(mjj + j + 1) =
>> 1034 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
>> 1035 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
>> 1036 temp1 = *(mjj + j);
>> 1037 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
>> 1038 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
>> 1039 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
>> 1040 iq = ip + pivrow-(j+1);
>> 1041 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 1042 {
>> 1043 temp1 = *iq;
>> 1044 *iq = *ip;
>> 1045 *ip = temp1;
>> 1046 }
>> 1047 }
>> 1048 // invert D(j:j+1,j:j+1)
>> 1049 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
>> 1050 if (temp2 == 0)
>> 1051 {
>> 1052 G4cerr
>> 1053 << "G4ErrorSymMatrix::bunch_invert: error in pivot choice"
>> 1054 << G4endl;
>> 1055 }
>> 1056 temp2 = 1. / temp2;
>> 1057
>> 1058 // this quotient is guaranteed to exist by the choice
>> 1059 // of the pivot
>> 1060
>> 1061 temp1 = *mjj;
>> 1062 *mjj = *(mjj + j + 1) * temp2;
>> 1063 *(mjj + j + 1) = temp1 * temp2;
>> 1064 *(mjj + j) = - *(mjj + j) * temp2;
>> 1065
>> 1066 if (j < nrow-1) // otherwise do nothing
>> 1067 {
>> 1068 // update A(j+2:n, j+2:n)
>> 1069 for (i=j+2; i <= nrow ; i++)
>> 1070 {
>> 1071 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1072 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
>> 1073 if (std::fabs(temp1 ) <= epsilon)
>> 1074 { temp1 = 0; }
>> 1075 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
>> 1076 if (std::fabs(temp2 ) <= epsilon)
>> 1077 { temp2 = 0; }
>> 1078 for (k = j+2; k <= i ; k++)
>> 1079 {
>> 1080 ip = m.begin() + i*(i-1)/2 + k-1;
>> 1081 iq = m.begin() + k*(k-1)/2 + j-1;
>> 1082 *ip -= temp1 * *iq + temp2 * *(iq+1);
>> 1083 if (std::fabs(*ip) <= epsilon)
>> 1084 { *ip = 0; }
>> 1085 }
>> 1086 }
>> 1087 // update L
>> 1088 for (i=j+2; i <= nrow ; i++)
>> 1089 {
>> 1090 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1091 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
>> 1092 if (std::fabs(temp1) <= epsilon)
>> 1093 { temp1 = 0; }
>> 1094 *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
>> 1095 if (std::fabs(*(ip+1)) <= epsilon)
>> 1096 { *(ip+1) = 0; }
>> 1097 *ip = temp1;
>> 1098 }
>> 1099 }
>> 1100 }
1167 } 1101 }
1168 *ip = temp1; << 1102 } // end of main loop over columns
1169 } <<
1170 } <<
1171 } <<
1172 } <<
1173 } // end of main loop over columns <<
1174 1103
1175 if(j == nrow) // the the last pivot is 1x1 << 1104 if (j == nrow) // the the last pivot is 1x1
1176 { 1105 {
1177 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1106 mjj = m.begin() + j*(j-1)/2 + j-1;
1178 if(*mjj == 0) << 1107 if (*mjj == 0)
1179 { 1108 {
1180 ifail = 1; 1109 ifail = 1;
1181 return; 1110 return;
1182 } 1111 }
1183 else 1112 else
1184 { 1113 {
1185 *mjj = 1. / *mjj; 1114 *mjj = 1. / *mjj;
1186 } 1115 }
1187 } // end of last pivot code << 1116 } // end of last pivot code
1188 1117
1189 // computing the inverse from the factoriza 1118 // computing the inverse from the factorization
1190 << 1119
1191 for(j = nrow; j >= 1; j -= ss) // loop ove << 1120 for (j = nrow ; j >= 1 ; j -= ss) // loop over columns
1192 { 1121 {
1193 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1122 mjj = m.begin() + j*(j-1)/2 + j-1;
1194 if(piv[j - 1] > 0) // 1x1 pivot, compute << 1123 if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
1195 { << 1124 {
1196 ss = 1; << 1125 ss = 1;
1197 if(j < nrow) << 1126 if (j < nrow)
1198 { << 1127 {
1199 ip = m.begin() + (j + 1) * j / 2 + j << 1128 ip = m.begin() + (j+1)*j/2 + j-1;
1200 for(i = 0; i < nrow - j; ip += 1 + j << 1129 for (i=0; i < nrow-j; ip += 1+j+i++)
1201 { << 1130 {
1202 x[i] = *ip; << 1131 x[i] = *ip;
1203 } << 1132 }
1204 for(i = j + 1; i <= nrow; i++) << 1133 for (i=j+1; i<=nrow ; i++)
1205 { << 1134 {
1206 temp2 = 0; << 1135 temp2=0;
1207 ip = m.begin() + i * (i - 1) / 2 << 1136 ip = m.begin() + i*(i-1)/2 + j;
1208 for(k = 0; k <= i - j - 1; k++) << 1137 for (k=0; k <= i-j-1; k++)
1209 { << 1138 { temp2 += *ip++ * x[k]; }
1210 temp2 += *ip++ * x[k]; << 1139 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1211 } << 1140 { temp2 += *ip * x[k]; }
1212 for(ip += i - 1; k < nrow - j; ip + << 1141 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1213 { << 1142 }
1214 temp2 += *ip * x[k]; << 1143 temp2 = 0;
1215 } << 1144 ip = m.begin() + (j+1)*j/2 + j-1;
1216 *(m.begin() + i * (i - 1) / 2 + j - << 1145 for (k=0; k < nrow-j; ip += 1+j+k++)
1217 } << 1146 { temp2 += x[k] * *ip; }
1218 temp2 = 0; << 1147 *mjj -= temp2;
1219 ip = m.begin() + (j + 1) * j / 2 + << 1148 }
1220 for(k = 0; k < nrow - j; ip += 1 + j << 1149 }
1221 { << 1150 else //2x2 pivot, compute columns j and j-1 of the inverse
1222 temp2 += x[k] * *ip; << 1151 {
1223 } << 1152 if (piv[j-1] != 0)
1224 *mjj -= temp2; << 1153 { G4cerr << "error in piv" << piv[j-1] << G4endl; }
1225 } << 1154 ss=2;
1226 } << 1155 if (j < nrow)
1227 else // 2x2 pivot, compute columns j and << 1156 {
1228 { << 1157 ip = m.begin() + (j+1)*j/2 + j-1;
1229 if(piv[j - 1] != 0) << 1158 for (i=0; i < nrow-j; ip += 1+j+i++)
1230 { << 1159 {
1231 std::ostringstream message; << 1160 x[i] = *ip;
1232 message << "Error in pivot: " << piv[ << 1161 }
1233 G4Exception("G4ErrorSymMatrix::invert << 1162 for (i=j+1; i<=nrow ; i++)
1234 "GEANT4e-Notification", J << 1163 {
1235 } << 1164 temp2 = 0;
1236 ss = 2; << 1165 ip = m.begin() + i*(i-1)/2 + j;
1237 if(j < nrow) << 1166 for (k=0; k <= i-j-1; k++)
1238 { << 1167 { temp2 += *ip++ * x[k]; }
1239 ip = m.begin() + (j + 1) * j / 2 + j << 1168 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1240 for(i = 0; i < nrow - j; ip += 1 + j << 1169 { temp2 += *ip * x[k]; }
1241 { << 1170 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1242 x[i] = *ip; << 1171 }
1243 } << 1172 temp2 = 0;
1244 for(i = j + 1; i <= nrow; i++) << 1173 ip = m.begin() + (j+1)*j/2 + j-1;
1245 { << 1174 for (k=0; k < nrow-j; ip += 1+j+k++)
1246 temp2 = 0; << 1175 { temp2 += x[k] * *ip; }
1247 ip = m.begin() + i * (i - 1) / 2 << 1176 *mjj -= temp2;
1248 for(k = 0; k <= i - j - 1; k++) << 1177 temp2 = 0;
1249 { << 1178 ip = m.begin() + (j+1)*j/2 + j-2;
1250 temp2 += *ip++ * x[k]; << 1179 for (i=j+1; i <= nrow; ip += i++)
1251 } << 1180 { temp2 += *ip * *(ip+1); }
1252 for(ip += i - 1; k < nrow - j; ip + << 1181 *(mjj-1) -= temp2;
1253 { << 1182 ip = m.begin() + (j+1)*j/2 + j-2;
1254 temp2 += *ip * x[k]; << 1183 for (i=0; i < nrow-j; ip += 1+j+i++)
1255 } << 1184 {
1256 *(m.begin() + i * (i - 1) / 2 + j - << 1185 x[i] = *ip;
1257 } << 1186 }
1258 temp2 = 0; << 1187 for (i=j+1; i <= nrow ; i++)
1259 ip = m.begin() + (j + 1) * j / 2 + << 1188 {
1260 for(k = 0; k < nrow - j; ip += 1 + j << 1189 temp2 = 0;
1261 { << 1190 ip = m.begin() + i*(i-1)/2 + j;
1262 temp2 += x[k] * *ip; << 1191 for (k=0; k <= i-j-1; k++)
1263 } << 1192 { temp2 += *ip++ * x[k]; }
1264 *mjj -= temp2; << 1193 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1265 temp2 = 0; << 1194 { temp2 += *ip * x[k]; }
1266 ip = m.begin() + (j + 1) * j / 2 + << 1195 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1267 for(i = j + 1; i <= nrow; ip += i++) << 1196 }
1268 { << 1197 temp2 = 0;
1269 temp2 += *ip * *(ip + 1); << 1198 ip = m.begin() + (j+1)*j/2 + j-2;
1270 } << 1199 for (k=0; k < nrow-j; ip += 1+j+k++)
1271 *(mjj - 1) -= temp2; << 1200 { temp2 += x[k] * *ip; }
1272 ip = m.begin() + (j + 1) * j / 2 + j << 1201 *(mjj-j) -= temp2;
1273 for(i = 0; i < nrow - j; ip += 1 + j << 1202 }
1274 { << 1203 }
1275 x[i] = *ip; << 1204
1276 } << 1205 // interchange rows and columns j and piv[j-1]
1277 for(i = j + 1; i <= nrow; i++) << 1206 // or rows and columns j and -piv[j-2]
1278 { << 1207
1279 temp2 = 0; << 1208 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1280 ip = m.begin() + i * (i - 1) / 2 << 1209 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1281 for(k = 0; k <= i - j - 1; k++) << 1210 for (i=j+1;i < pivrow; i++, ip++)
1282 { << 1211 {
1283 temp2 += *ip++ * x[k]; << 1212 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1284 } << 1213 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1285 for(ip += i - 1; k < nrow - j; ip + << 1214 *ip = temp1;
1286 { << 1215 }
1287 temp2 += *ip * x[k]; << 1216 temp1 = *mjj;
1288 } << 1217 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1289 *(m.begin() + i * (i - 1) / 2 + j - << 1218 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1290 } << 1219 if (ss==2)
1291 temp2 = 0; << 1220 {
1292 ip = m.begin() + (j + 1) * j / 2 + << 1221 temp1 = *(mjj-1);
1293 for(k = 0; k < nrow - j; ip += 1 + j << 1222 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1294 { << 1223 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1295 temp2 += x[k] * *ip; << 1224 }
1296 } << 1225
1297 *(mjj - j) -= temp2; << 1226 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1298 } << 1227 iq = ip + pivrow-j;
1299 } << 1228 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1300 << 1229 {
1301 // interchange rows and columns j and piv << 1230 temp1 = *iq;
1302 // or rows and columns j and -piv[j-2] << 1231 *iq = *ip;
1303 << 1232 *ip = temp1;
1304 pivrow = (piv[j - 1] == 0) ? -piv[j - 2] << 1233 }
1305 ip = m.begin() + pivrow * (pivrow - 1 << 1234 } // end of loop over columns (in computing inverse from factorization)
1306 for(i = j + 1; i < pivrow; i++, ip++) << 1235
1307 { << 1236 return; // inversion successful
1308 temp1 = *(m.begin() + i * (i - 1) / 2 + << 1237 }
1309 *(m.begin() + i * (i - 1) / 2 + j - 1) << 1238
1310 *ip << 1239 G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
1311 } << 1240 G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
1312 temp1 = *mjj; << 1241 G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
1313 *mjj = *(m.begin() + pivrow * (pivrow - << 1242 G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
1314 *(m.begin() + pivrow * (pivrow - 1) / 2 + << 1243 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
1315 if(ss == 2) << 1244 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
1316 { << 1245 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
1317 temp1 = *(mjj - 1); << 1246 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;
1318 *(mjj - 1) = *(m.begin() + pivrow * (pi <<
1319 *(m.begin() + pivrow * (pivrow - 1) / 2 <<
1320 } <<
1321 <<
1322 ip = m.begin() + (pivrow + 1) * pivrow / <<
1323 iq = ip + pivrow - j; <<
1324 for(i = pivrow + 1; i <= nrow; ip += i, i <<
1325 { <<
1326 temp1 = *iq; <<
1327 *iq = *ip; <<
1328 *ip = temp1; <<
1329 } <<
1330 } // end of loop over columns (in computin <<
1331 <<
1332 return; // inversion successful <<
1333 } <<
1334 <<
1335 G4ThreadLocal G4double G4ErrorSymMatrix::posD <<
1336 G4ThreadLocal G4double G4ErrorSymMatrix::posD <<
1337 G4ThreadLocal G4double G4ErrorSymMatrix::adju <<
1338 G4ThreadLocal G4double G4ErrorSymMatrix::adju <<
1339 const G4double G4ErrorSymMatrix::CHOLESKY_THR <<
1340 const G4double G4ErrorSymMatrix::CHOLESKY_THR <<
1341 const G4double G4ErrorSymMatrix::CHOLESKY_CRE <<
1342 const G4double G4ErrorSymMatrix::CHOLESKY_CRE <<
1343 1247
1344 // Aij are indices for a 6x6 symmetric matrix 1248 // Aij are indices for a 6x6 symmetric matrix.
1345 // The indices for 5x5 or 4x4 symmetric m << 1249 // The indices for 5x5 or 4x4 symmetric matrices are the same,
1346 // ignoring all combinations with an inde 1250 // ignoring all combinations with an index which is inapplicable.
1347 1251
1348 #define A00 0 1252 #define A00 0
1349 #define A01 1 1253 #define A01 1
1350 #define A02 3 1254 #define A02 3
1351 #define A03 6 1255 #define A03 6
1352 #define A04 10 1256 #define A04 10
1353 #define A05 15 1257 #define A05 15
1354 1258
1355 #define A10 1 1259 #define A10 1
1356 #define A11 2 1260 #define A11 2
1357 #define A12 4 1261 #define A12 4
1358 #define A13 7 1262 #define A13 7
1359 #define A14 11 1263 #define A14 11
1360 #define A15 16 1264 #define A15 16
1361 1265
1362 #define A20 3 1266 #define A20 3
1363 #define A21 4 1267 #define A21 4
1364 #define A22 5 1268 #define A22 5
1365 #define A23 8 1269 #define A23 8
1366 #define A24 12 1270 #define A24 12
1367 #define A25 17 1271 #define A25 17
1368 1272
1369 #define A30 6 1273 #define A30 6
1370 #define A31 7 1274 #define A31 7
1371 #define A32 8 1275 #define A32 8
1372 #define A33 9 1276 #define A33 9
1373 #define A34 13 1277 #define A34 13
1374 #define A35 18 1278 #define A35 18
1375 1279
1376 #define A40 10 1280 #define A40 10
1377 #define A41 11 1281 #define A41 11
1378 #define A42 12 1282 #define A42 12
1379 #define A43 13 1283 #define A43 13
1380 #define A44 14 1284 #define A44 14
1381 #define A45 19 1285 #define A45 19
1382 1286
1383 #define A50 15 1287 #define A50 15
1384 #define A51 16 1288 #define A51 16
1385 #define A52 17 1289 #define A52 17
1386 #define A53 18 1290 #define A53 18
1387 #define A54 19 1291 #define A54 19
1388 #define A55 20 1292 #define A55 20
1389 1293
1390 void G4ErrorSymMatrix::invert5(G4int& ifail) << 1294 void G4ErrorSymMatrix::invert5(G4int & ifail)
1391 { << 1295 {
1392 if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_ << 1296 if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
1393 { << 1297 {
1394 invertCholesky5(ifail); << 1298 invertCholesky5(ifail);
1395 posDefFraction5x5 = .9 * posDefFraction5x << 1299 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
1396 if(ifail != 0) // Cholesky failed -- inv << 1300 if (ifail!=0) // Cholesky failed -- invert using Haywood
1397 { << 1301 {
1398 invertHaywood5(ifail); << 1302 invertHaywood5(ifail);
1399 } << 1303 }
1400 } << 1304 }
1401 else << 1305 else
1402 { <<
1403 if(posDefFraction5x5 + adjustment5x5 >= C <<
1404 { <<
1405 invertCholesky5(ifail); <<
1406 posDefFraction5x5 = .9 * posDefFraction <<
1407 if(ifail != 0) // Cholesky failed -- i <<
1408 { 1306 {
1409 invertHaywood5(ifail); << 1307 if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
1410 adjustment5x5 = 0; << 1308 {
>> 1309 invertCholesky5(ifail);
>> 1310 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
>> 1311 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1312 {
>> 1313 invertHaywood5(ifail);
>> 1314 adjustment5x5 = 0;
>> 1315 }
>> 1316 }
>> 1317 else
>> 1318 {
>> 1319 invertHaywood5(ifail);
>> 1320 adjustment5x5 += CHOLESKY_CREEP_5x5;
>> 1321 }
1411 } 1322 }
1412 } << 1323 return;
1413 else <<
1414 { <<
1415 invertHaywood5(ifail); <<
1416 adjustment5x5 += CHOLESKY_CREEP_5x5; <<
1417 } <<
1418 } <<
1419 return; <<
1420 } 1324 }
1421 1325
1422 void G4ErrorSymMatrix::invert6(G4int& ifail) << 1326 void G4ErrorSymMatrix::invert6(G4int & ifail)
1423 { 1327 {
1424 if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_ << 1328 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 { 1329 {
1441 invertHaywood6(ifail); << 1330 invertCholesky6(ifail);
1442 adjustment6x6 = 0; << 1331 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
>> 1332 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1333 {
>> 1334 invertHaywood6(ifail);
>> 1335 }
1443 } 1336 }
1444 } << 1337 else
1445 else << 1338 {
1446 { << 1339 if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
1447 invertHaywood6(ifail); << 1340 {
1448 adjustment6x6 += CHOLESKY_CREEP_6x6; << 1341 invertCholesky6(ifail);
1449 } << 1342 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
1450 } << 1343 if (ifail!=0) // Cholesky failed -- invert using Haywood
1451 return; << 1344 {
>> 1345 invertHaywood6(ifail);
>> 1346 adjustment6x6 = 0;
>> 1347 }
>> 1348 }
>> 1349 else
>> 1350 {
>> 1351 invertHaywood6(ifail);
>> 1352 adjustment6x6 += CHOLESKY_CREEP_6x6;
>> 1353 }
>> 1354 }
>> 1355 return;
1452 } 1356 }
1453 1357
1454 void G4ErrorSymMatrix::invertHaywood5(G4int& << 1358 void G4ErrorSymMatrix::invertHaywood5 (G4int & ifail)
1455 { 1359 {
1456 ifail = 0; 1360 ifail = 0;
1457 1361
1458 // Find all NECESSARY 2x2 dets: (25 of the 1362 // Find all NECESSARY 2x2 dets: (25 of them)
1459 1363
1460 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 1364 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
1461 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 1365 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
1462 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 1366 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
1463 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 1367 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
1464 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 1368 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
1465 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 1369 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
1466 G4double Det2_24_01 = m[A20] * m[A41] - m[A << 1370 G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
1467 G4double Det2_24_02 = m[A20] * m[A42] - m[A << 1371 G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
1468 G4double Det2_24_03 = m[A20] * m[A43] - m[A << 1372 G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
1469 G4double Det2_24_04 = m[A20] * m[A44] - m[A << 1373 G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
1470 G4double Det2_24_12 = m[A21] * m[A42] - m[A << 1374 G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
1471 G4double Det2_24_13 = m[A21] * m[A43] - m[A << 1375 G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
1472 G4double Det2_24_14 = m[A21] * m[A44] - m[A << 1376 G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
1473 G4double Det2_24_23 = m[A22] * m[A43] - m[A << 1377 G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
1474 G4double Det2_24_24 = m[A22] * m[A44] - m[A << 1378 G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
1475 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1379 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1476 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1380 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1477 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1381 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1478 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1382 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1479 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1383 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1480 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1384 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1481 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1385 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1482 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1386 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1483 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1387 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1484 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1388 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1485 1389
1486 // Find all NECESSARY 3x3 dets: (30 of th 1390 // Find all NECESSARY 3x3 dets: (30 of them)
1487 1391
1488 G4double Det3_123_012 = << 1392 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 << 1393 + m[A12]*Det2_23_01;
1490 G4double Det3_123_013 = << 1394 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 << 1395 + m[A13]*Det2_23_01;
1492 G4double Det3_123_023 = << 1396 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 << 1397 + m[A13]*Det2_23_02;
1494 G4double Det3_123_123 = << 1398 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 << 1399 + m[A13]*Det2_23_12;
1496 G4double Det3_124_012 = << 1400 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 << 1401 + m[A12]*Det2_24_01;
1498 G4double Det3_124_013 = << 1402 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 << 1403 + m[A13]*Det2_24_01;
1500 G4double Det3_124_014 = << 1404 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 << 1405 + m[A14]*Det2_24_01;
1502 G4double Det3_124_023 = << 1406 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 << 1407 + m[A13]*Det2_24_02;
1504 G4double Det3_124_024 = << 1408 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 << 1409 + m[A14]*Det2_24_02;
1506 G4double Det3_124_123 = << 1410 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 << 1411 + m[A13]*Det2_24_12;
1508 G4double Det3_124_124 = << 1412 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 << 1413 + m[A14]*Det2_24_12;
1510 G4double Det3_134_012 = << 1414 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 << 1415 + m[A12]*Det2_34_01;
1512 G4double Det3_134_013 = << 1416 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 << 1417 + m[A13]*Det2_34_01;
1514 G4double Det3_134_014 = << 1418 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 << 1419 + m[A14]*Det2_34_01;
1516 G4double Det3_134_023 = << 1420 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 << 1421 + m[A13]*Det2_34_02;
1518 G4double Det3_134_024 = << 1422 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 << 1423 + m[A14]*Det2_34_02;
1520 G4double Det3_134_034 = << 1424 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 << 1425 + m[A14]*Det2_34_03;
1522 G4double Det3_134_123 = << 1426 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 << 1427 + m[A13]*Det2_34_12;
1524 G4double Det3_134_124 = << 1428 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 << 1429 + m[A14]*Det2_34_12;
1526 G4double Det3_134_134 = << 1430 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 << 1431 + m[A14]*Det2_34_13;
1528 G4double Det3_234_012 = << 1432 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 << 1433 + m[A22]*Det2_34_01;
1530 G4double Det3_234_013 = << 1434 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 << 1435 + m[A23]*Det2_34_01;
1532 G4double Det3_234_014 = << 1436 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 << 1437 + m[A24]*Det2_34_01;
1534 G4double Det3_234_023 = << 1438 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 << 1439 + m[A23]*Det2_34_02;
1536 G4double Det3_234_024 = << 1440 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 << 1441 + m[A24]*Det2_34_02;
1538 G4double Det3_234_034 = << 1442 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 << 1443 + m[A24]*Det2_34_03;
1540 G4double Det3_234_123 = << 1444 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 << 1445 + m[A23]*Det2_34_12;
1542 G4double Det3_234_124 = << 1446 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 << 1447 + m[A24]*Det2_34_12;
1544 G4double Det3_234_134 = << 1448 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 << 1449 + m[A24]*Det2_34_13;
1546 G4double Det3_234_234 = << 1450 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 << 1451 + m[A24]*Det2_34_23;
1548 1452
1549 // Find all NECESSARY 4x4 dets: (15 of th 1453 // Find all NECESSARY 4x4 dets: (15 of them)
1550 1454
1551 G4double Det4_0123_0123 = m[A00] * Det3_123 << 1455 G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023
1552 m[A02] * Det3_123 << 1456 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
1553 G4double Det4_0124_0123 = m[A00] * Det3_124 << 1457 G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023
1554 m[A02] * Det3_124 << 1458 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
1555 G4double Det4_0124_0124 = m[A00] * Det3_124 << 1459 G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024
1556 m[A02] * Det3_124 << 1460 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
1557 G4double Det4_0134_0123 = m[A00] * Det3_134 << 1461 G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023
1558 m[A02] * Det3_134 << 1462 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
1559 G4double Det4_0134_0124 = m[A00] * Det3_134 << 1463 G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024
1560 m[A02] * Det3_134 << 1464 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
1561 G4double Det4_0134_0134 = m[A00] * Det3_134 << 1465 G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034
1562 m[A03] * Det3_134 << 1466 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
1563 G4double Det4_0234_0123 = m[A00] * Det3_234 << 1467 G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023
1564 m[A02] * Det3_234 << 1468 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
1565 G4double Det4_0234_0124 = m[A00] * Det3_234 << 1469 G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024
1566 m[A02] * Det3_234 << 1470 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
1567 G4double Det4_0234_0134 = m[A00] * Det3_234 << 1471 G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034
1568 m[A03] * Det3_234 << 1472 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
1569 G4double Det4_0234_0234 = m[A00] * Det3_234 << 1473 G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034
1570 m[A03] * Det3_234 << 1474 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
1571 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1475 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1572 m[A12] * Det3_234 << 1476 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1573 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1477 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1574 m[A12] * Det3_234 << 1478 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1575 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1479 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1576 m[A13] * Det3_234 << 1480 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1577 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1481 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1578 m[A13] * Det3_234 << 1482 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1579 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1483 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1580 m[A13] * Det3_234 << 1484 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1581 1485
1582 // Find the 5x5 det: 1486 // Find the 5x5 det:
1583 1487
1584 G4double det = m[A00] * Det4_1234_1234 - m[ << 1488 G4double det = m[A00]*Det4_1234_1234
1585 m[A02] * Det4_1234_0134 - m[ << 1489 - m[A01]*Det4_1234_0234
1586 m[A04] * Det4_1234_0123; << 1490 + m[A02]*Det4_1234_0134
>> 1491 - m[A03]*Det4_1234_0124
>> 1492 + m[A04]*Det4_1234_0123;
1587 1493
1588 if(det == 0) << 1494 if ( det == 0 )
1589 { << 1495 {
1590 ifail = 1; 1496 ifail = 1;
1591 return; 1497 return;
1592 } << 1498 }
1593 1499
1594 G4double oneOverDet = 1.0 / det; << 1500 G4double oneOverDet = 1.0/det;
1595 G4double mn1OverDet = -oneOverDet; << 1501 G4double mn1OverDet = - oneOverDet;
1596 1502
1597 m[A00] = Det4_1234_1234 * oneOverDet; << 1503 m[A00] = Det4_1234_1234 * oneOverDet;
1598 m[A01] = Det4_1234_0234 * mn1OverDet; << 1504 m[A01] = Det4_1234_0234 * mn1OverDet;
1599 m[A02] = Det4_1234_0134 * oneOverDet; << 1505 m[A02] = Det4_1234_0134 * oneOverDet;
1600 m[A03] = Det4_1234_0124 * mn1OverDet; << 1506 m[A03] = Det4_1234_0124 * mn1OverDet;
1601 m[A04] = Det4_1234_0123 * oneOverDet; << 1507 m[A04] = Det4_1234_0123 * oneOverDet;
1602 << 1508
1603 m[A11] = Det4_0234_0234 * oneOverDet; << 1509 m[A11] = Det4_0234_0234 * oneOverDet;
1604 m[A12] = Det4_0234_0134 * mn1OverDet; << 1510 m[A12] = Det4_0234_0134 * mn1OverDet;
1605 m[A13] = Det4_0234_0124 * oneOverDet; << 1511 m[A13] = Det4_0234_0124 * oneOverDet;
1606 m[A14] = Det4_0234_0123 * mn1OverDet; << 1512 m[A14] = Det4_0234_0123 * mn1OverDet;
1607 << 1513
1608 m[A22] = Det4_0134_0134 * oneOverDet; << 1514 m[A22] = Det4_0134_0134 * oneOverDet;
1609 m[A23] = Det4_0134_0124 * mn1OverDet; << 1515 m[A23] = Det4_0134_0124 * mn1OverDet;
1610 m[A24] = Det4_0134_0123 * oneOverDet; << 1516 m[A24] = Det4_0134_0123 * oneOverDet;
1611 1517
1612 m[A33] = Det4_0124_0124 * oneOverDet; << 1518 m[A33] = Det4_0124_0124 * oneOverDet;
1613 m[A34] = Det4_0124_0123 * mn1OverDet; << 1519 m[A34] = Det4_0124_0123 * mn1OverDet;
1614 1520
1615 m[A44] = Det4_0123_0123 * oneOverDet; << 1521 m[A44] = Det4_0123_0123 * oneOverDet;
1616 1522
1617 return; 1523 return;
1618 } 1524 }
1619 1525
1620 void G4ErrorSymMatrix::invertHaywood6(G4int& << 1526 void G4ErrorSymMatrix::invertHaywood6 (G4int & ifail)
1621 { 1527 {
1622 ifail = 0; 1528 ifail = 0;
1623 1529
1624 // Find all NECESSARY 2x2 dets: (39 of the 1530 // Find all NECESSARY 2x2 dets: (39 of them)
1625 1531
1626 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1532 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1627 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1533 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1628 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1534 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1629 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1535 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1630 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1536 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1631 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1537 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1632 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1538 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1633 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1539 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1634 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1540 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1635 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1541 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1636 G4double Det2_35_01 = m[A30] * m[A51] - m[A << 1542 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1637 G4double Det2_35_02 = m[A30] * m[A52] - m[A << 1543 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1638 G4double Det2_35_03 = m[A30] * m[A53] - m[A << 1544 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1639 G4double Det2_35_04 = m[A30] * m[A54] - m[A << 1545 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1640 G4double Det2_35_05 = m[A30] * m[A55] - m[A << 1546 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1641 G4double Det2_35_12 = m[A31] * m[A52] - m[A << 1547 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1642 G4double Det2_35_13 = m[A31] * m[A53] - m[A << 1548 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1643 G4double Det2_35_14 = m[A31] * m[A54] - m[A << 1549 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1644 G4double Det2_35_15 = m[A31] * m[A55] - m[A << 1550 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1645 G4double Det2_35_23 = m[A32] * m[A53] - m[A << 1551 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1646 G4double Det2_35_24 = m[A32] * m[A54] - m[A << 1552 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1647 G4double Det2_35_25 = m[A32] * m[A55] - m[A << 1553 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1648 G4double Det2_35_34 = m[A33] * m[A54] - m[A << 1554 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1649 G4double Det2_35_35 = m[A33] * m[A55] - m[A << 1555 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1650 G4double Det2_45_01 = m[A40] * m[A51] - m[A << 1556 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1651 G4double Det2_45_02 = m[A40] * m[A52] - m[A << 1557 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1652 G4double Det2_45_03 = m[A40] * m[A53] - m[A << 1558 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1653 G4double Det2_45_04 = m[A40] * m[A54] - m[A << 1559 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1654 G4double Det2_45_05 = m[A40] * m[A55] - m[A << 1560 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1655 G4double Det2_45_12 = m[A41] * m[A52] - m[A << 1561 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1656 G4double Det2_45_13 = m[A41] * m[A53] - m[A << 1562 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1657 G4double Det2_45_14 = m[A41] * m[A54] - m[A << 1563 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1658 G4double Det2_45_15 = m[A41] * m[A55] - m[A << 1564 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1659 G4double Det2_45_23 = m[A42] * m[A53] - m[A << 1565 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1660 G4double Det2_45_24 = m[A42] * m[A54] - m[A << 1566 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1661 G4double Det2_45_25 = m[A42] * m[A55] - m[A << 1567 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1662 G4double Det2_45_34 = m[A43] * m[A54] - m[A << 1568 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1663 G4double Det2_45_35 = m[A43] * m[A55] - m[A << 1569 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1664 G4double Det2_45_45 = m[A44] * m[A55] - m[A << 1570 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1665 1571
1666 // Find all NECESSARY 3x3 dets: (65 of the 1572 // Find all NECESSARY 3x3 dets: (65 of them)
1667 1573
1668 G4double Det3_234_012 = << 1574 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 << 1575 + m[A22]*Det2_34_01;
1670 G4double Det3_234_013 = << 1576 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 << 1577 + m[A23]*Det2_34_01;
1672 G4double Det3_234_014 = << 1578 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 << 1579 + m[A24]*Det2_34_01;
1674 G4double Det3_234_023 = << 1580 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 << 1581 + m[A23]*Det2_34_02;
1676 G4double Det3_234_024 = << 1582 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 << 1583 + m[A24]*Det2_34_02;
1678 G4double Det3_234_034 = << 1584 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 << 1585 + m[A24]*Det2_34_03;
1680 G4double Det3_234_123 = << 1586 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 << 1587 + m[A23]*Det2_34_12;
1682 G4double Det3_234_124 = << 1588 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 << 1589 + m[A24]*Det2_34_12;
1684 G4double Det3_234_134 = << 1590 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 << 1591 + m[A24]*Det2_34_13;
1686 G4double Det3_234_234 = << 1592 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 << 1593 + m[A24]*Det2_34_23;
1688 G4double Det3_235_012 = << 1594 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 << 1595 + m[A22]*Det2_35_01;
1690 G4double Det3_235_013 = << 1596 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 << 1597 + m[A23]*Det2_35_01;
1692 G4double Det3_235_014 = << 1598 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 << 1599 + m[A24]*Det2_35_01;
1694 G4double Det3_235_015 = << 1600 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 << 1601 + m[A25]*Det2_35_01;
1696 G4double Det3_235_023 = << 1602 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 << 1603 + m[A23]*Det2_35_02;
1698 G4double Det3_235_024 = << 1604 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 << 1605 + m[A24]*Det2_35_02;
1700 G4double Det3_235_025 = << 1606 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 << 1607 + m[A25]*Det2_35_02;
1702 G4double Det3_235_034 = << 1608 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 << 1609 + m[A24]*Det2_35_03;
1704 G4double Det3_235_035 = << 1610 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 << 1611 + m[A25]*Det2_35_03;
1706 G4double Det3_235_123 = << 1612 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 << 1613 + m[A23]*Det2_35_12;
1708 G4double Det3_235_124 = << 1614 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 << 1615 + m[A24]*Det2_35_12;
1710 G4double Det3_235_125 = << 1616 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 << 1617 + m[A25]*Det2_35_12;
1712 G4double Det3_235_134 = << 1618 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 << 1619 + m[A24]*Det2_35_13;
1714 G4double Det3_235_135 = << 1620 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 << 1621 + m[A25]*Det2_35_13;
1716 G4double Det3_235_234 = << 1622 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 << 1623 + m[A24]*Det2_35_23;
1718 G4double Det3_235_235 = << 1624 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 << 1625 + m[A25]*Det2_35_23;
1720 G4double Det3_245_012 = << 1626 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 << 1627 + m[A22]*Det2_45_01;
1722 G4double Det3_245_013 = << 1628 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 << 1629 + m[A23]*Det2_45_01;
1724 G4double Det3_245_014 = << 1630 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 << 1631 + m[A24]*Det2_45_01;
1726 G4double Det3_245_015 = << 1632 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 << 1633 + m[A25]*Det2_45_01;
1728 G4double Det3_245_023 = << 1634 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 << 1635 + m[A23]*Det2_45_02;
1730 G4double Det3_245_024 = << 1636 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 << 1637 + m[A24]*Det2_45_02;
1732 G4double Det3_245_025 = << 1638 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 << 1639 + m[A25]*Det2_45_02;
1734 G4double Det3_245_034 = << 1640 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 << 1641 + m[A24]*Det2_45_03;
1736 G4double Det3_245_035 = << 1642 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 << 1643 + m[A25]*Det2_45_03;
1738 G4double Det3_245_045 = << 1644 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 << 1645 + m[A25]*Det2_45_04;
1740 G4double Det3_245_123 = << 1646 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 << 1647 + m[A23]*Det2_45_12;
1742 G4double Det3_245_124 = << 1648 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 << 1649 + m[A24]*Det2_45_12;
1744 G4double Det3_245_125 = << 1650 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 << 1651 + m[A25]*Det2_45_12;
1746 G4double Det3_245_134 = << 1652 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 << 1653 + m[A24]*Det2_45_13;
1748 G4double Det3_245_135 = << 1654 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 << 1655 + m[A25]*Det2_45_13;
1750 G4double Det3_245_145 = << 1656 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 << 1657 + m[A25]*Det2_45_14;
1752 G4double Det3_245_234 = << 1658 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 << 1659 + m[A24]*Det2_45_23;
1754 G4double Det3_245_235 = << 1660 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 << 1661 + m[A25]*Det2_45_23;
1756 G4double Det3_245_245 = << 1662 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 << 1663 + m[A25]*Det2_45_24;
1758 G4double Det3_345_012 = << 1664 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 << 1665 + m[A32]*Det2_45_01;
1760 G4double Det3_345_013 = << 1666 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 << 1667 + m[A33]*Det2_45_01;
1762 G4double Det3_345_014 = << 1668 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 << 1669 + m[A34]*Det2_45_01;
1764 G4double Det3_345_015 = << 1670 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 << 1671 + m[A35]*Det2_45_01;
1766 G4double Det3_345_023 = << 1672 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 << 1673 + m[A33]*Det2_45_02;
1768 G4double Det3_345_024 = << 1674 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 << 1675 + m[A34]*Det2_45_02;
1770 G4double Det3_345_025 = << 1676 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 << 1677 + m[A35]*Det2_45_02;
1772 G4double Det3_345_034 = << 1678 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 << 1679 + m[A34]*Det2_45_03;
1774 G4double Det3_345_035 = << 1680 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 << 1681 + m[A35]*Det2_45_03;
1776 G4double Det3_345_045 = << 1682 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 << 1683 + m[A35]*Det2_45_04;
1778 G4double Det3_345_123 = << 1684 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 << 1685 + m[A33]*Det2_45_12;
1780 G4double Det3_345_124 = << 1686 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 << 1687 + m[A34]*Det2_45_12;
1782 G4double Det3_345_125 = << 1688 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 << 1689 + m[A35]*Det2_45_12;
1784 G4double Det3_345_134 = << 1690 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 << 1691 + m[A34]*Det2_45_13;
1786 G4double Det3_345_135 = << 1692 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 << 1693 + m[A35]*Det2_45_13;
1788 G4double Det3_345_145 = << 1694 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 << 1695 + m[A35]*Det2_45_14;
1790 G4double Det3_345_234 = << 1696 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 << 1697 + m[A34]*Det2_45_23;
1792 G4double Det3_345_235 = << 1698 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 << 1699 + m[A35]*Det2_45_23;
1794 G4double Det3_345_245 = << 1700 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 << 1701 + m[A35]*Det2_45_24;
1796 G4double Det3_345_345 = << 1702 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 << 1703 + m[A35]*Det2_45_34;
1798 1704
1799 // Find all NECESSARY 4x4 dets: (55 of the 1705 // Find all NECESSARY 4x4 dets: (55 of them)
1800 1706
1801 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1707 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1802 m[A12] * Det3_234 << 1708 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1803 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1709 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1804 m[A12] * Det3_234 << 1710 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1805 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1711 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1806 m[A13] * Det3_234 << 1712 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1807 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1713 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1808 m[A13] * Det3_234 << 1714 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1809 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1715 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1810 m[A13] * Det3_234 << 1716 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1811 G4double Det4_1235_0123 = m[A10] * Det3_235 << 1717 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1812 m[A12] * Det3_235 << 1718 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1813 G4double Det4_1235_0124 = m[A10] * Det3_235 << 1719 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1814 m[A12] * Det3_235 << 1720 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1815 G4double Det4_1235_0125 = m[A10] * Det3_235 << 1721 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1816 m[A12] * Det3_235 << 1722 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1817 G4double Det4_1235_0134 = m[A10] * Det3_235 << 1723 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1818 m[A13] * Det3_235 << 1724 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1819 G4double Det4_1235_0135 = m[A10] * Det3_235 << 1725 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1820 m[A13] * Det3_235 << 1726 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1821 G4double Det4_1235_0234 = m[A10] * Det3_235 << 1727 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1822 m[A13] * Det3_235 << 1728 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1823 G4double Det4_1235_0235 = m[A10] * Det3_235 << 1729 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1824 m[A13] * Det3_235 << 1730 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1825 G4double Det4_1235_1234 = m[A11] * Det3_235 << 1731 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1826 m[A13] * Det3_235 << 1732 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1827 G4double Det4_1235_1235 = m[A11] * Det3_235 << 1733 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1828 m[A13] * Det3_235 << 1734 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1829 G4double Det4_1245_0123 = m[A10] * Det3_245 << 1735 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1830 m[A12] * Det3_245 << 1736 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1831 G4double Det4_1245_0124 = m[A10] * Det3_245 << 1737 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1832 m[A12] * Det3_245 << 1738 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1833 G4double Det4_1245_0125 = m[A10] * Det3_245 << 1739 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1834 m[A12] * Det3_245 << 1740 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1835 G4double Det4_1245_0134 = m[A10] * Det3_245 << 1741 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1836 m[A13] * Det3_245 << 1742 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1837 G4double Det4_1245_0135 = m[A10] * Det3_245 << 1743 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1838 m[A13] * Det3_245 << 1744 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1839 G4double Det4_1245_0145 = m[A10] * Det3_245 << 1745 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1840 m[A14] * Det3_245 << 1746 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1841 G4double Det4_1245_0234 = m[A10] * Det3_245 << 1747 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1842 m[A13] * Det3_245 << 1748 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1843 G4double Det4_1245_0235 = m[A10] * Det3_245 << 1749 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1844 m[A13] * Det3_245 << 1750 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1845 G4double Det4_1245_0245 = m[A10] * Det3_245 << 1751 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1846 m[A14] * Det3_245 << 1752 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1847 G4double Det4_1245_1234 = m[A11] * Det3_245 << 1753 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1848 m[A13] * Det3_245 << 1754 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1849 G4double Det4_1245_1235 = m[A11] * Det3_245 << 1755 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1850 m[A13] * Det3_245 << 1756 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1851 G4double Det4_1245_1245 = m[A11] * Det3_245 << 1757 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1852 m[A14] * Det3_245 << 1758 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1853 G4double Det4_1345_0123 = m[A10] * Det3_345 << 1759 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1854 m[A12] * Det3_345 << 1760 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1855 G4double Det4_1345_0124 = m[A10] * Det3_345 << 1761 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1856 m[A12] * Det3_345 << 1762 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1857 G4double Det4_1345_0125 = m[A10] * Det3_345 << 1763 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1858 m[A12] * Det3_345 << 1764 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1859 G4double Det4_1345_0134 = m[A10] * Det3_345 << 1765 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1860 m[A13] * Det3_345 << 1766 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1861 G4double Det4_1345_0135 = m[A10] * Det3_345 << 1767 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1862 m[A13] * Det3_345 << 1768 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1863 G4double Det4_1345_0145 = m[A10] * Det3_345 << 1769 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1864 m[A14] * Det3_345 << 1770 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1865 G4double Det4_1345_0234 = m[A10] * Det3_345 << 1771 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1866 m[A13] * Det3_345 << 1772 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1867 G4double Det4_1345_0235 = m[A10] * Det3_345 << 1773 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1868 m[A13] * Det3_345 << 1774 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1869 G4double Det4_1345_0245 = m[A10] * Det3_345 << 1775 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1870 m[A14] * Det3_345 << 1776 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1871 G4double Det4_1345_0345 = m[A10] * Det3_345 << 1777 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1872 m[A14] * Det3_345 << 1778 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1873 G4double Det4_1345_1234 = m[A11] * Det3_345 << 1779 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1874 m[A13] * Det3_345 << 1780 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1875 G4double Det4_1345_1235 = m[A11] * Det3_345 << 1781 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1876 m[A13] * Det3_345 << 1782 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1877 G4double Det4_1345_1245 = m[A11] * Det3_345 << 1783 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1878 m[A14] * Det3_345 << 1784 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1879 G4double Det4_1345_1345 = m[A11] * Det3_345 << 1785 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1880 m[A14] * Det3_345 << 1786 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1881 G4double Det4_2345_0123 = m[A20] * Det3_345 << 1787 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1882 m[A22] * Det3_345 << 1788 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1883 G4double Det4_2345_0124 = m[A20] * Det3_345 << 1789 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1884 m[A22] * Det3_345 << 1790 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1885 G4double Det4_2345_0125 = m[A20] * Det3_345 << 1791 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1886 m[A22] * Det3_345 << 1792 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1887 G4double Det4_2345_0134 = m[A20] * Det3_345 << 1793 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1888 m[A23] * Det3_345 << 1794 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1889 G4double Det4_2345_0135 = m[A20] * Det3_345 << 1795 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1890 m[A23] * Det3_345 << 1796 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1891 G4double Det4_2345_0145 = m[A20] * Det3_345 << 1797 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1892 m[A24] * Det3_345 << 1798 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1893 G4double Det4_2345_0234 = m[A20] * Det3_345 << 1799 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1894 m[A23] * Det3_345 << 1800 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1895 G4double Det4_2345_0235 = m[A20] * Det3_345 << 1801 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1896 m[A23] * Det3_345 << 1802 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1897 G4double Det4_2345_0245 = m[A20] * Det3_345 << 1803 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1898 m[A24] * Det3_345 << 1804 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1899 G4double Det4_2345_0345 = m[A20] * Det3_345 << 1805 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1900 m[A24] * Det3_345 << 1806 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1901 G4double Det4_2345_1234 = m[A21] * Det3_345 << 1807 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1902 m[A23] * Det3_345 << 1808 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1903 G4double Det4_2345_1235 = m[A21] * Det3_345 << 1809 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1904 m[A23] * Det3_345 << 1810 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1905 G4double Det4_2345_1245 = m[A21] * Det3_345 << 1811 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1906 m[A24] * Det3_345 << 1812 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1907 G4double Det4_2345_1345 = m[A21] * Det3_345 << 1813 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1908 m[A24] * Det3_345 << 1814 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1909 G4double Det4_2345_2345 = m[A22] * Det3_345 << 1815 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1910 m[A24] * Det3_345 << 1816 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1911 1817
1912 // Find all NECESSARY 5x5 dets: (19 of the 1818 // Find all NECESSARY 5x5 dets: (19 of them)
1913 1819
1914 G4double Det5_01234_01234 = << 1820 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1915 m[A00] * Det4_1234_1234 - m[A01] * Det4_1 << 1821 + 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 << 1822 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1917 G4double Det5_01235_01234 = << 1823 + 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 << 1824 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1919 m[A02] * Det4_1235_0134 - m[A03] * Det4_1 << 1825 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1920 G4double Det5_01235_01235 = << 1826 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1921 m[A00] * Det4_1235_1235 - m[A01] * Det4_1 << 1827 + 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 << 1828 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1923 G4double Det5_01245_01234 = << 1829 + 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 << 1830 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1925 m[A02] * Det4_1245_0134 - m[A03] * Det4_1 << 1831 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1926 G4double Det5_01245_01235 = << 1832 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1927 m[A00] * Det4_1245_1235 - m[A01] * Det4_1 << 1833 + 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 << 1834 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1929 G4double Det5_01245_01245 = << 1835 + 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 << 1836 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1931 m[A02] * Det4_1245_0145 - m[A04] * Det4_1 << 1837 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1932 G4double Det5_01345_01234 = << 1838 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1933 m[A00] * Det4_1345_1234 - m[A01] * Det4_1 << 1839 + 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 << 1840 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1935 G4double Det5_01345_01235 = << 1841 + 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 << 1842 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1937 m[A02] * Det4_1345_0135 - m[A03] * Det4_1 << 1843 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1938 G4double Det5_01345_01245 = << 1844 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1939 m[A00] * Det4_1345_1245 - m[A01] * Det4_1 << 1845 + 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 << 1846 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1941 G4double Det5_01345_01345 = << 1847 + 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 << 1848 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1943 m[A03] * Det4_1345_0145 - m[A04] * Det4_1 << 1849 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1944 G4double Det5_02345_01234 = << 1850 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1945 m[A00] * Det4_2345_1234 - m[A01] * Det4_2 << 1851 + 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 << 1852 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1947 G4double Det5_02345_01235 = << 1853 + 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 << 1854 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1949 m[A02] * Det4_2345_0135 - m[A03] * Det4_2 << 1855 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1950 G4double Det5_02345_01245 = << 1856 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1951 m[A00] * Det4_2345_1245 - m[A01] * Det4_2 << 1857 + 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 << 1858 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1953 G4double Det5_02345_01345 = << 1859 + 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 << 1860 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1955 m[A03] * Det4_2345_0145 - m[A04] * Det4_2 << 1861 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1956 G4double Det5_02345_02345 = << 1862
1957 m[A00] * Det4_2345_2345 - m[A02] * Det4_2 << 1863 // Find the determinant
1958 m[A03] * Det4_2345_0245 - m[A04] * Det4_2 << 1864
1959 G4double Det5_12345_01234 = << 1865 G4double det = m[A00]*Det5_12345_12345
1960 m[A10] * Det4_2345_1234 - m[A11] * Det4_2 << 1866 - m[A01]*Det5_12345_02345
1961 m[A12] * Det4_2345_0134 - m[A13] * Det4_2 << 1867 + m[A02]*Det5_12345_01345
1962 G4double Det5_12345_01235 = << 1868 - m[A03]*Det5_12345_01245
1963 m[A10] * Det4_2345_1235 - m[A11] * Det4_2 << 1869 + m[A04]*Det5_12345_01235
1964 m[A12] * Det4_2345_0135 - m[A13] * Det4_2 << 1870 - 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 1871
1984 if(det == 0) << 1872 if ( det == 0 )
1985 { << 1873 {
1986 ifail = 1; 1874 ifail = 1;
1987 return; 1875 return;
1988 } << 1876 }
1989 1877
1990 G4double oneOverDet = 1.0 / det; << 1878 G4double oneOverDet = 1.0/det;
1991 G4double mn1OverDet = -oneOverDet; << 1879 G4double mn1OverDet = - oneOverDet;
1992 1880
1993 m[A00] = Det5_12345_12345 * oneOverDet; << 1881 m[A00] = Det5_12345_12345*oneOverDet;
1994 m[A01] = Det5_12345_02345 * mn1OverDet; << 1882 m[A01] = Det5_12345_02345*mn1OverDet;
1995 m[A02] = Det5_12345_01345 * oneOverDet; << 1883 m[A02] = Det5_12345_01345*oneOverDet;
1996 m[A03] = Det5_12345_01245 * mn1OverDet; << 1884 m[A03] = Det5_12345_01245*mn1OverDet;
1997 m[A04] = Det5_12345_01235 * oneOverDet; << 1885 m[A04] = Det5_12345_01235*oneOverDet;
1998 m[A05] = Det5_12345_01234 * mn1OverDet; << 1886 m[A05] = Det5_12345_01234*mn1OverDet;
1999 << 1887
2000 m[A11] = Det5_02345_02345 * oneOverDet; << 1888 m[A11] = Det5_02345_02345*oneOverDet;
2001 m[A12] = Det5_02345_01345 * mn1OverDet; << 1889 m[A12] = Det5_02345_01345*mn1OverDet;
2002 m[A13] = Det5_02345_01245 * oneOverDet; << 1890 m[A13] = Det5_02345_01245*oneOverDet;
2003 m[A14] = Det5_02345_01235 * mn1OverDet; << 1891 m[A14] = Det5_02345_01235*mn1OverDet;
2004 m[A15] = Det5_02345_01234 * oneOverDet; << 1892 m[A15] = Det5_02345_01234*oneOverDet;
2005 << 1893
2006 m[A22] = Det5_01345_01345 * oneOverDet; << 1894 m[A22] = Det5_01345_01345*oneOverDet;
2007 m[A23] = Det5_01345_01245 * mn1OverDet; << 1895 m[A23] = Det5_01345_01245*mn1OverDet;
2008 m[A24] = Det5_01345_01235 * oneOverDet; << 1896 m[A24] = Det5_01345_01235*oneOverDet;
2009 m[A25] = Det5_01345_01234 * mn1OverDet; << 1897 m[A25] = Det5_01345_01234*mn1OverDet;
2010 << 1898
2011 m[A33] = Det5_01245_01245 * oneOverDet; << 1899 m[A33] = Det5_01245_01245*oneOverDet;
2012 m[A34] = Det5_01245_01235 * mn1OverDet; << 1900 m[A34] = Det5_01245_01235*mn1OverDet;
2013 m[A35] = Det5_01245_01234 * oneOverDet; << 1901 m[A35] = Det5_01245_01234*oneOverDet;
2014 1902
2015 m[A44] = Det5_01235_01235 * oneOverDet; << 1903 m[A44] = Det5_01235_01235*oneOverDet;
2016 m[A45] = Det5_01235_01234 * mn1OverDet; << 1904 m[A45] = Det5_01235_01234*mn1OverDet;
2017 1905
2018 m[A55] = Det5_01234_01234 * oneOverDet; << 1906 m[A55] = Det5_01234_01234*oneOverDet;
2019 1907
2020 return; 1908 return;
2021 } 1909 }
2022 1910
2023 void G4ErrorSymMatrix::invertCholesky5(G4int& << 1911 void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
2024 { 1912 {
2025 // Invert by << 1913
>> 1914 // Invert by
2026 // 1915 //
2027 // a) decomposing M = G*G^T with G lower tr 1916 // a) decomposing M = G*G^T with G lower triangular
2028 // (if M is not positive definite this w 1917 // (if M is not positive definite this will fail, leaving this unchanged)
2029 // b) inverting G to form H 1918 // b) inverting G to form H
2030 // c) multiplying H^T * H to get M^-1. 1919 // c) multiplying H^T * H to get M^-1.
2031 // 1920 //
2032 // If the matrix is pos. def. it is inverte 1921 // If the matrix is pos. def. it is inverted and 1 is returned.
2033 // If the matrix is not pos. def. it remain 1922 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2034 1923
2035 G4double h10; // below-diagonal elements o << 1924 G4double h10; // below-diagonal elements of H
2036 G4double h20, h21; 1925 G4double h20, h21;
2037 G4double h30, h31, h32; 1926 G4double h30, h31, h32;
2038 G4double h40, h41, h42, h43; 1927 G4double h40, h41, h42, h43;
2039 1928
2040 G4double h00, h11, h22, h33, h44; // 1/dia << 1929 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
2041 // diago << 1930 // diagonal elements of H
2042 1931
2043 G4double g10; // below-diagonal elements o << 1932 G4double g10; // below-diagonal elements of G
2044 G4double g20, g21; 1933 G4double g20, g21;
2045 G4double g30, g31, g32; 1934 G4double g30, g31, g32;
2046 G4double g40, g41, g42, g43; 1935 G4double g40, g41, g42, g43;
2047 1936
2048 ifail = 1; // We start by assuing we won't 1937 ifail = 1; // We start by assuing we won't succeed...
2049 1938
2050 // Form G -- compute diagonal members of H 1939 // Form G -- compute diagonal members of H directly rather than of G
2051 //------- 1940 //-------
2052 1941
2053 // Scale first column by 1/sqrt(A00) 1942 // Scale first column by 1/sqrt(A00)
2054 1943
2055 h00 = m[A00]; << 1944 h00 = m[A00];
2056 if(h00 <= 0) << 1945 if (h00 <= 0) { return; }
2057 { <<
2058 return; <<
2059 } <<
2060 h00 = 1.0 / std::sqrt(h00); 1946 h00 = 1.0 / std::sqrt(h00);
2061 1947
2062 g10 = m[A10] * h00; 1948 g10 = m[A10] * h00;
2063 g20 = m[A20] * h00; 1949 g20 = m[A20] * h00;
2064 g30 = m[A30] * h00; 1950 g30 = m[A30] * h00;
2065 g40 = m[A40] * h00; 1951 g40 = m[A40] * h00;
2066 1952
2067 // Form G11 (actually, just h11) 1953 // Form G11 (actually, just h11)
2068 1954
2069 h11 = m[A11] - (g10 * g10); 1955 h11 = m[A11] - (g10 * g10);
2070 if(h11 <= 0) << 1956 if (h11 <= 0) { return; }
2071 { <<
2072 return; <<
2073 } <<
2074 h11 = 1.0 / std::sqrt(h11); 1957 h11 = 1.0 / std::sqrt(h11);
2075 1958
2076 // Subtract inter-column column dot product 1959 // Subtract inter-column column dot products from rest of column 1 and
2077 // scale to get column 1 of G << 1960 // scale to get column 1 of G
2078 1961
2079 g21 = (m[A21] - (g10 * g20)) * h11; 1962 g21 = (m[A21] - (g10 * g20)) * h11;
2080 g31 = (m[A31] - (g10 * g30)) * h11; 1963 g31 = (m[A31] - (g10 * g30)) * h11;
2081 g41 = (m[A41] - (g10 * g40)) * h11; 1964 g41 = (m[A41] - (g10 * g40)) * h11;
2082 1965
2083 // Form G22 (actually, just h22) 1966 // Form G22 (actually, just h22)
2084 1967
2085 h22 = m[A22] - (g20 * g20) - (g21 * g21); 1968 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2086 if(h22 <= 0) << 1969 if (h22 <= 0) { return; }
2087 { <<
2088 return; <<
2089 } <<
2090 h22 = 1.0 / std::sqrt(h22); 1970 h22 = 1.0 / std::sqrt(h22);
2091 1971
2092 // Subtract inter-column column dot product 1972 // Subtract inter-column column dot products from rest of column 2 and
2093 // scale to get column 2 of G << 1973 // scale to get column 2 of G
2094 1974
2095 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 1975 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2096 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 1976 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2097 1977
2098 // Form G33 (actually, just h33) 1978 // Form G33 (actually, just h33)
2099 1979
2100 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 1980 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2101 if(h33 <= 0) << 1981 if (h33 <= 0) { return; }
2102 { <<
2103 return; <<
2104 } <<
2105 h33 = 1.0 / std::sqrt(h33); 1982 h33 = 1.0 / std::sqrt(h33);
2106 1983
2107 // Subtract inter-column column dot product 1984 // Subtract inter-column column dot product from A43 and scale to get G43
2108 1985
2109 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 1986 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2110 1987
2111 // Finally form h44 - if this is possible i 1988 // Finally form h44 - if this is possible inversion succeeds
2112 1989
2113 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 1990 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2114 if(h44 <= 0) << 1991 if (h44 <= 0) { return; }
2115 { <<
2116 return; <<
2117 } <<
2118 h44 = 1.0 / std::sqrt(h44); 1992 h44 = 1.0 / std::sqrt(h44);
2119 1993
2120 // Form H = 1/G -- diagonal members of H ar 1994 // Form H = 1/G -- diagonal members of H are already correct
2121 //------------- 1995 //-------------
2122 1996
2123 // The order here is dictated by speed cons 1997 // The order here is dictated by speed considerations
2124 1998
2125 h43 = -h33 * g43 * h44; << 1999 h43 = -h33 * g43 * h44;
2126 h32 = -h22 * g32 * h33; << 2000 h32 = -h22 * g32 * h33;
2127 h42 = -h22 * (g32 * h43 + g42 * h44); 2001 h42 = -h22 * (g32 * h43 + g42 * h44);
2128 h21 = -h11 * g21 * h22; << 2002 h21 = -h11 * g21 * h22;
2129 h31 = -h11 * (g21 * h32 + g31 * h33); 2003 h31 = -h11 * (g21 * h32 + g31 * h33);
2130 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2004 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2131 h10 = -h00 * g10 * h11; << 2005 h10 = -h00 * g10 * h11;
2132 h20 = -h00 * (g10 * h21 + g20 * h22); 2006 h20 = -h00 * (g10 * h21 + g20 * h22);
2133 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2007 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2134 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2008 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2135 2009
2136 // Change this to its inverse = H^T*H 2010 // Change this to its inverse = H^T*H
2137 //------------------------------------ 2011 //------------------------------------
2138 2012
2139 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 2013 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2140 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 2014 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2141 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2015 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
2142 m[A02] = h20 * h22 + h30 * h32 + h40 * h42; 2016 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2143 m[A12] = h21 * h22 + h31 * h32 + h41 * h42; 2017 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2144 m[A22] = h22 * h22 + h32 * h32 + h42 * h42; 2018 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2145 m[A03] = h30 * h33 + h40 * h43; 2019 m[A03] = h30 * h33 + h40 * h43;
2146 m[A13] = h31 * h33 + h41 * h43; 2020 m[A13] = h31 * h33 + h41 * h43;
2147 m[A23] = h32 * h33 + h42 * h43; 2021 m[A23] = h32 * h33 + h42 * h43;
2148 m[A33] = h33 * h33 + h43 * h43; 2022 m[A33] = h33 * h33 + h43 * h43;
2149 m[A04] = h40 * h44; 2023 m[A04] = h40 * h44;
2150 m[A14] = h41 * h44; 2024 m[A14] = h41 * h44;
2151 m[A24] = h42 * h44; 2025 m[A24] = h42 * h44;
2152 m[A34] = h43 * h44; 2026 m[A34] = h43 * h44;
2153 m[A44] = h44 * h44; 2027 m[A44] = h44 * h44;
2154 2028
2155 ifail = 0; 2029 ifail = 0;
2156 return; 2030 return;
2157 } 2031 }
2158 2032
2159 void G4ErrorSymMatrix::invertCholesky6(G4int& << 2033 void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
2160 { 2034 {
2161 // Invert by << 2035 // Invert by
2162 // 2036 //
2163 // a) decomposing M = G*G^T with G lower tr 2037 // a) decomposing M = G*G^T with G lower triangular
2164 // (if M is not positive definite this w 2038 // (if M is not positive definite this will fail, leaving this unchanged)
2165 // b) inverting G to form H 2039 // b) inverting G to form H
2166 // c) multiplying H^T * H to get M^-1. 2040 // c) multiplying H^T * H to get M^-1.
2167 // 2041 //
2168 // If the matrix is pos. def. it is inverte 2042 // If the matrix is pos. def. it is inverted and 1 is returned.
2169 // If the matrix is not pos. def. it remain 2043 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2170 2044
2171 G4double h10; // below-diagonal elements o << 2045 G4double h10; // below-diagonal elements of H
2172 G4double h20, h21; 2046 G4double h20, h21;
2173 G4double h30, h31, h32; 2047 G4double h30, h31, h32;
2174 G4double h40, h41, h42, h43; 2048 G4double h40, h41, h42, h43;
2175 G4double h50, h51, h52, h53, h54; 2049 G4double h50, h51, h52, h53, h54;
2176 2050
2177 G4double h00, h11, h22, h33, h44, h55; // << 2051 G4double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G =
2178 // 2052 // diagonal elements of H
2179 2053
2180 G4double g10; // below-diagonal elements o << 2054 G4double g10; // below-diagonal elements of G
2181 G4double g20, g21; 2055 G4double g20, g21;
2182 G4double g30, g31, g32; 2056 G4double g30, g31, g32;
2183 G4double g40, g41, g42, g43; 2057 G4double g40, g41, g42, g43;
2184 G4double g50, g51, g52, g53, g54; 2058 G4double g50, g51, g52, g53, g54;
2185 2059
2186 ifail = 1; // We start by assuing we won't 2060 ifail = 1; // We start by assuing we won't succeed...
2187 2061
2188 // Form G -- compute diagonal members of H 2062 // Form G -- compute diagonal members of H directly rather than of G
2189 //------- 2063 //-------
2190 2064
2191 // Scale first column by 1/sqrt(A00) 2065 // Scale first column by 1/sqrt(A00)
2192 2066
2193 h00 = m[A00]; << 2067 h00 = m[A00];
2194 if(h00 <= 0) << 2068 if (h00 <= 0) { return; }
2195 { <<
2196 return; <<
2197 } <<
2198 h00 = 1.0 / std::sqrt(h00); 2069 h00 = 1.0 / std::sqrt(h00);
2199 2070
2200 g10 = m[A10] * h00; 2071 g10 = m[A10] * h00;
2201 g20 = m[A20] * h00; 2072 g20 = m[A20] * h00;
2202 g30 = m[A30] * h00; 2073 g30 = m[A30] * h00;
2203 g40 = m[A40] * h00; 2074 g40 = m[A40] * h00;
2204 g50 = m[A50] * h00; 2075 g50 = m[A50] * h00;
2205 2076
2206 // Form G11 (actually, just h11) 2077 // Form G11 (actually, just h11)
2207 2078
2208 h11 = m[A11] - (g10 * g10); 2079 h11 = m[A11] - (g10 * g10);
2209 if(h11 <= 0) << 2080 if (h11 <= 0) { return; }
2210 { <<
2211 return; <<
2212 } <<
2213 h11 = 1.0 / std::sqrt(h11); 2081 h11 = 1.0 / std::sqrt(h11);
2214 2082
2215 // Subtract inter-column column dot product 2083 // Subtract inter-column column dot products from rest of column 1 and
2216 // scale to get column 1 of G << 2084 // scale to get column 1 of G
2217 2085
2218 g21 = (m[A21] - (g10 * g20)) * h11; 2086 g21 = (m[A21] - (g10 * g20)) * h11;
2219 g31 = (m[A31] - (g10 * g30)) * h11; 2087 g31 = (m[A31] - (g10 * g30)) * h11;
2220 g41 = (m[A41] - (g10 * g40)) * h11; 2088 g41 = (m[A41] - (g10 * g40)) * h11;
2221 g51 = (m[A51] - (g10 * g50)) * h11; 2089 g51 = (m[A51] - (g10 * g50)) * h11;
2222 2090
2223 // Form G22 (actually, just h22) 2091 // Form G22 (actually, just h22)
2224 2092
2225 h22 = m[A22] - (g20 * g20) - (g21 * g21); 2093 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2226 if(h22 <= 0) << 2094 if (h22 <= 0) { return; }
2227 { <<
2228 return; <<
2229 } <<
2230 h22 = 1.0 / std::sqrt(h22); 2095 h22 = 1.0 / std::sqrt(h22);
2231 2096
2232 // Subtract inter-column column dot product 2097 // Subtract inter-column column dot products from rest of column 2 and
2233 // scale to get column 2 of G << 2098 // scale to get column 2 of G
2234 2099
2235 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 2100 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2236 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 2101 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2237 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) 2102 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2238 2103
2239 // Form G33 (actually, just h33) 2104 // Form G33 (actually, just h33)
2240 2105
2241 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 2106 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2242 if(h33 <= 0) << 2107 if (h33 <= 0) { return; }
2243 { <<
2244 return; <<
2245 } <<
2246 h33 = 1.0 / std::sqrt(h33); 2108 h33 = 1.0 / std::sqrt(h33);
2247 2109
2248 // Subtract inter-column column dot product 2110 // Subtract inter-column column dot products from rest of column 3 and
2249 // scale to get column 3 of G << 2111 // scale to get column 3 of G
2250 2112
2251 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 2113 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2252 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - 2114 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2253 2115
2254 // Form G44 (actually, just h44) 2116 // Form G44 (actually, just h44)
2255 2117
2256 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 2118 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2257 if(h44 <= 0) << 2119 if (h44 <= 0) { return; }
2258 { <<
2259 return; <<
2260 } <<
2261 h44 = 1.0 / std::sqrt(h44); 2120 h44 = 1.0 / std::sqrt(h44);
2262 2121
2263 // Subtract inter-column column dot product 2122 // Subtract inter-column column dot product from M54 and scale to get G54
2264 2123
2265 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - 2124 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2266 2125
2267 // Finally form h55 - if this is possible i 2126 // Finally form h55 - if this is possible inversion succeeds
2268 2127
2269 h55 = m[A55] - (g50 * g50) - (g51 * g51) - << 2128 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
2270 (g54 * g54); << 2129 if (h55 <= 0) { return; }
2271 if(h55 <= 0) <<
2272 { <<
2273 return; <<
2274 } <<
2275 h55 = 1.0 / std::sqrt(h55); 2130 h55 = 1.0 / std::sqrt(h55);
2276 2131
2277 // Form H = 1/G -- diagonal members of H ar 2132 // Form H = 1/G -- diagonal members of H are already correct
2278 //------------- 2133 //-------------
2279 2134
2280 // The order here is dictated by speed cons 2135 // The order here is dictated by speed considerations
2281 2136
2282 h54 = -h44 * g54 * h55; << 2137 h54 = -h44 * g54 * h55;
2283 h43 = -h33 * g43 * h44; << 2138 h43 = -h33 * g43 * h44;
2284 h53 = -h33 * (g43 * h54 + g53 * h55); 2139 h53 = -h33 * (g43 * h54 + g53 * h55);
2285 h32 = -h22 * g32 * h33; << 2140 h32 = -h22 * g32 * h33;
2286 h42 = -h22 * (g32 * h43 + g42 * h44); 2141 h42 = -h22 * (g32 * h43 + g42 * h44);
2287 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * 2142 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
2288 h21 = -h11 * g21 * h22; << 2143 h21 = -h11 * g21 * h22;
2289 h31 = -h11 * (g21 * h32 + g31 * h33); 2144 h31 = -h11 * (g21 * h32 + g31 * h33);
2290 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2145 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2291 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * 2146 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2292 h10 = -h00 * g10 * h11; << 2147 h10 = -h00 * g10 * h11;
2293 h20 = -h00 * (g10 * h21 + g20 * h22); 2148 h20 = -h00 * (g10 * h21 + g20 * h22);
2294 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2149 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2295 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2150 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2296 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * 2151 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2297 2152
2298 // Change this to its inverse = H^T*H 2153 // Change this to its inverse = H^T*H
2299 //------------------------------------ 2154 //------------------------------------
2300 2155
2301 m[A00] = << 2156 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 2157 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2304 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2158 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2305 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 2159 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2306 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 2160 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2307 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 2161 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
2308 m[A03] = h30 * h33 + h40 * h43 + h50 * h53; 2162 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2309 m[A13] = h31 * h33 + h41 * h43 + h51 * h53; 2163 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2310 m[A23] = h32 * h33 + h42 * h43 + h52 * h53; 2164 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2311 m[A33] = h33 * h33 + h43 * h43 + h53 * h53; 2165 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2312 m[A04] = h40 * h44 + h50 * h54; 2166 m[A04] = h40 * h44 + h50 * h54;
2313 m[A14] = h41 * h44 + h51 * h54; 2167 m[A14] = h41 * h44 + h51 * h54;
2314 m[A24] = h42 * h44 + h52 * h54; 2168 m[A24] = h42 * h44 + h52 * h54;
2315 m[A34] = h43 * h44 + h53 * h54; 2169 m[A34] = h43 * h44 + h53 * h54;
2316 m[A44] = h44 * h44 + h54 * h54; 2170 m[A44] = h44 * h44 + h54 * h54;
2317 m[A05] = h50 * h55; 2171 m[A05] = h50 * h55;
2318 m[A15] = h51 * h55; 2172 m[A15] = h51 * h55;
2319 m[A25] = h52 * h55; 2173 m[A25] = h52 * h55;
2320 m[A35] = h53 * h55; 2174 m[A35] = h53 * h55;
2321 m[A45] = h54 * h55; 2175 m[A45] = h54 * h55;
2322 m[A55] = h55 * h55; 2176 m[A55] = h55 * h55;
2323 2177
2324 ifail = 0; 2178 ifail = 0;
2325 return; 2179 return;
2326 } 2180 }
2327 2181
2328 void G4ErrorSymMatrix::invert4(G4int& ifail) << 2182 void G4ErrorSymMatrix::invert4 (G4int & ifail)
2329 { 2183 {
2330 ifail = 0; 2184 ifail = 0;
2331 2185
2332 // Find all NECESSARY 2x2 dets: (14 of the 2186 // Find all NECESSARY 2x2 dets: (14 of them)
2333 2187
2334 G4double Det2_12_01 = m[A10] * m[A21] - m[A << 2188 G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
2335 G4double Det2_12_02 = m[A10] * m[A22] - m[A << 2189 G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
2336 G4double Det2_12_12 = m[A11] * m[A22] - m[A << 2190 G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
2337 G4double Det2_13_01 = m[A10] * m[A31] - m[A << 2191 G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
2338 G4double Det2_13_02 = m[A10] * m[A32] - m[A << 2192 G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
2339 G4double Det2_13_03 = m[A10] * m[A33] - m[A << 2193 G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
2340 G4double Det2_13_12 = m[A11] * m[A32] - m[A << 2194 G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
2341 G4double Det2_13_13 = m[A11] * m[A33] - m[A << 2195 G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
2342 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 2196 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
2343 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 2197 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
2344 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 2198 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
2345 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 2199 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
2346 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 2200 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
2347 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 2201 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
2348 2202
2349 // Find all NECESSARY 3x3 dets: (10 of th 2203 // Find all NECESSARY 3x3 dets: (10 of them)
2350 2204
2351 G4double Det3_012_012 = << 2205 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 << 2206 + m[A02]*Det2_12_01;
2353 G4double Det3_013_012 = << 2207 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 << 2208 + m[A02]*Det2_13_01;
2355 G4double Det3_013_013 = << 2209 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 << 2210 + m[A03]*Det2_13_01;
2357 G4double Det3_023_012 = << 2211 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 << 2212 + m[A02]*Det2_23_01;
2359 G4double Det3_023_013 = << 2213 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 << 2214 + m[A03]*Det2_23_01;
2361 G4double Det3_023_023 = << 2215 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 << 2216 + m[A03]*Det2_23_02;
2363 G4double Det3_123_012 = << 2217 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 << 2218 + m[A12]*Det2_23_01;
2365 G4double Det3_123_013 = << 2219 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 << 2220 + m[A13]*Det2_23_01;
2367 G4double Det3_123_023 = << 2221 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 << 2222 + m[A13]*Det2_23_02;
2369 G4double Det3_123_123 = << 2223 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 << 2224 + m[A13]*Det2_23_12;
2371 2225
2372 // Find the 4x4 det: 2226 // Find the 4x4 det:
2373 2227
2374 G4double det = m[A00] * Det3_123_123 - m[A0 << 2228 G4double det = m[A00]*Det3_123_123
2375 m[A02] * Det3_123_013 - m[A0 << 2229 - m[A01]*Det3_123_023
>> 2230 + m[A02]*Det3_123_013
>> 2231 - m[A03]*Det3_123_012;
2376 2232
2377 if(det == 0) << 2233 if ( det == 0 )
2378 { << 2234 {
2379 ifail = 1; 2235 ifail = 1;
2380 return; 2236 return;
2381 } << 2237 }
>> 2238
>> 2239 G4double oneOverDet = 1.0/det;
>> 2240 G4double mn1OverDet = - oneOverDet;
2382 2241
2383 G4double oneOverDet = 1.0 / det; << 2242 m[A00] = Det3_123_123 * oneOverDet;
2384 G4double mn1OverDet = -oneOverDet; << 2243 m[A01] = Det3_123_023 * mn1OverDet;
>> 2244 m[A02] = Det3_123_013 * oneOverDet;
>> 2245 m[A03] = Det3_123_012 * mn1OverDet;
2385 2246
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 2247
2395 m[A22] = Det3_013_013 * oneOverDet; << 2248 m[A11] = Det3_023_023 * oneOverDet;
2396 m[A23] = Det3_013_012 * mn1OverDet; << 2249 m[A12] = Det3_023_013 * mn1OverDet;
>> 2250 m[A13] = Det3_023_012 * oneOverDet;
2397 2251
2398 m[A33] = Det3_012_012 * oneOverDet; << 2252 m[A22] = Det3_013_013 * oneOverDet;
>> 2253 m[A23] = Det3_013_012 * mn1OverDet;
>> 2254
>> 2255 m[A33] = Det3_012_012 * oneOverDet;
2399 2256
2400 return; 2257 return;
2401 } 2258 }
2402 2259
2403 void G4ErrorSymMatrix::invertHaywood4(G4int& << 2260 void G4ErrorSymMatrix::invertHaywood4 (G4int & ifail)
2404 { 2261 {
2405 invert4(ifail); // For the 4x4 case, the m << 2262 invert4(ifail); // For the 4x4 case, the method we use for invert is already
2406 // the Haywood method. << 2263 // the Haywood method.
2407 } 2264 }
>> 2265
2408 2266