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