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