//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id: G4ErrorSymMatrix.cc,v 1.3 2007/06/21 15:04:10 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
// ------------------------------------------------------------
//      GEANT 4 class implementation file
// ------------------------------------------------------------

#include "globals.hh"
#include <iostream>
#include <cmath>

#include "G4ErrorSymMatrix.hh"
#include "G4ErrorMatrix.hh"

// Simple operation for all elements

#define SIMPLE_UOP(OPER)                      \
    G4ErrorMatrixIter a=m.begin();            \
    G4ErrorMatrixIter e=m.begin()+num_size(); \
   for(;a<e; a++) (*a) OPER t;

#define SIMPLE_BOP(OPER)                           \
    G4ErrorMatrixIter a=m.begin();                 \
    G4ErrorMatrixConstIter b=m2.m.begin();         \
    G4ErrorMatrixConstIter e=m.begin()+num_size(); \
   for(;a<e; a++, b++) (*a) OPER (*b);

#define SIMPLE_TOP(OPER)                                 \
    G4ErrorMatrixConstIter a=m1.m.begin();               \
    G4ErrorMatrixConstIter b=m2.m.begin();               \
    G4ErrorMatrixIter t=mret.m.begin();                  \
    G4ErrorMatrixConstIter e=m1.m.begin()+m1.num_size(); \
   for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);

#define CHK_DIM_2(r1,r2,c1,c2,fun) \
   if (r1!=r2 || c1!=c2)  { \
    G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
   }

#define CHK_DIM_1(c1,r2,fun) \
   if (c1!=r2) { \
     G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
   }

// Constructors. (Default constructors are inlined and in .icc file)

G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
   : m(p*(p+1)/2), nrow(p)
{
   size = nrow * (nrow+1) / 2;
   m.assign(size,0);
}

G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4int init)
   : m(p*(p+1)/2), nrow(p)
{
   size = nrow * (nrow+1) / 2;

   m.assign(size,0);
   switch(init)
   {
   case 0:
      break;
      
   case 1:
      {
         G4ErrorMatrixIter a = m.begin();
         for(G4int i=1;i<=nrow;i++)
         {
            *a = 1.0;
            a += (i+1);
         }
         break;
      }
   default:
      G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
   }
}

//
// Destructor
//

G4ErrorSymMatrix::~G4ErrorSymMatrix()
{
}

G4ErrorSymMatrix::G4ErrorSymMatrix(const G4ErrorSymMatrix &m1)
   : m(m1.size), nrow(m1.nrow), size(m1.size)
{
   m = m1.m;
}

//
//
// Sub matrix
//
//

G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) const
{
  G4ErrorSymMatrix mret(max_row-min_row+1);
  if(max_row > num_row())
    { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
  G4ErrorMatrixIter a = mret.m.begin();
  G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
  for(G4int irow=1; irow<=mret.num_row(); irow++)
  {
    G4ErrorMatrixConstIter b = b1;
    for(G4int icol=1; icol<=irow; icol++)
    {
      *(a++) = *(b++);
    }
    b1 += irow+min_row-1;
  }
  return mret;
}

G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) 
{
  G4ErrorSymMatrix mret(max_row-min_row+1);
  if(max_row > num_row())
    { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
  G4ErrorMatrixIter a = mret.m.begin();
  G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
  for(G4int irow=1; irow<=mret.num_row(); irow++)
  {
    G4ErrorMatrixIter b = b1;
    for(G4int icol=1; icol<=irow; icol++)
    {
      *(a++) = *(b++);
    }
    b1 += irow+min_row-1;
  }
  return mret;
}

void G4ErrorSymMatrix::sub(G4int row,const G4ErrorSymMatrix &m1)
{
  if(row <1 || row+m1.num_row()-1 > num_row() )
    { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
  G4ErrorMatrixConstIter a = m1.m.begin();
  G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
  for(G4int irow=1; irow<=m1.num_row(); irow++)
  {
    G4ErrorMatrixIter b = b1;
    for(G4int icol=1; icol<=irow; icol++)
    {
      *(b++) = *(a++);
    }
    b1 += irow+row-1;
  }
}

//
// Direct sum of two matricies
//

G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &m1,
                      const G4ErrorSymMatrix &m2)
{
  G4ErrorSymMatrix mret(m1.num_row() + m2.num_row(), 0);
  mret.sub(1,m1);
  mret.sub(m1.num_row()+1,m2);
  return mret;
}

/* -----------------------------------------------------------------------
   This section contains support routines for matrix.h. This section contains
   The two argument functions +,-. They call the copy constructor and +=,-=.
   ----------------------------------------------------------------------- */

G4ErrorSymMatrix G4ErrorSymMatrix::operator- () const 
{
  G4ErrorSymMatrix m2(nrow);
  G4ErrorMatrixConstIter a=m.begin();
  G4ErrorMatrixIter b=m2.m.begin();
  G4ErrorMatrixConstIter e=m.begin()+num_size();
  for(;a<e; a++, b++) { (*b) = -(*a); }
  return m2;
}

G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
{
  G4ErrorMatrix mret(m1);
  CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+);
  mret += m2;
  return mret;
}

G4ErrorMatrix operator+(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
{
  G4ErrorMatrix mret(m2);
  CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),+);
  mret += m1;
  return mret;
}

G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1,
                           const G4ErrorSymMatrix &m2)
{
  G4ErrorSymMatrix mret(m1.nrow);
  CHK_DIM_1(m1.nrow, m2.nrow,+);
  SIMPLE_TOP(+)
  return mret;
}

//
// operator -
//

G4ErrorMatrix operator-(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
{
  G4ErrorMatrix mret(m1);
  CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
  mret -= m2;
  return mret;
}

G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
{
  G4ErrorMatrix mret(m1);
  CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
  mret -= m2;
  return mret;
}

G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1,
                           const G4ErrorSymMatrix &m2)
{
  G4ErrorSymMatrix mret(m1.num_row());
  CHK_DIM_1(m1.num_row(),m2.num_row(),-);
  SIMPLE_TOP(-)
  return mret;
}

/* -----------------------------------------------------------------------
   This section contains support routines for matrix.h. This file contains
   The two argument functions *,/. They call copy constructor and then /=,*=.
   ----------------------------------------------------------------------- */

G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1,G4double t)
{
  G4ErrorSymMatrix mret(m1);
  mret /= t;
  return mret;
}

G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &m1,G4double t)
{
  G4ErrorSymMatrix mret(m1);
  mret *= t;
  return mret;
}

G4ErrorSymMatrix operator*(G4double t,const G4ErrorSymMatrix &m1)
{
  G4ErrorSymMatrix mret(m1);
  mret *= t;
  return mret;
}

G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
{
  G4ErrorMatrix mret(m1.num_row(),m2.num_col());
  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
  G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
  G4double temp;
  G4ErrorMatrixIter mir=mret.m.begin();
  G4int step,stept;
  for(mit1=m1.m.begin();
      mit1<m1.m.begin()+m1.num_row()*m1.num_col(); mit1 = mit2)
  {
    for(step=1,snp=m2.m.begin();step<=m2.num_row();)
    {
      mit2=mit1;
      sp=snp;
      snp+=step;
      temp=0;
      while(sp<snp)
      {
        temp+=*(sp++)*(*(mit2++));
      }
      sp+=step-1;
      for(stept=++step;stept<=m2.num_row();stept++)
      {
        temp+=*sp*(*(mit2++));
        sp+=stept;
      }
      *(mir++)=temp;
    }
  }
  return mret;
}

G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
{
  G4ErrorMatrix mret(m1.num_row(),m2.num_col());
  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
  G4int step,stept;
  G4ErrorMatrixConstIter mit1,mit2,sp,snp;
  G4double temp;
  G4ErrorMatrixIter mir=mret.m.begin();
  for(step=1,snp=m1.m.begin();step<=m1.num_row();snp+=step++)
  {
    for(mit1=m2.m.begin();mit1<m2.m.begin()+m2.num_col();mit1++)
    {
      mit2=mit1;
      sp=snp;
      temp=0;
      while(sp<snp+step)
      {
        temp+=*mit2*(*(sp++));
        mit2+=m2.num_col();
      }
      sp+=step-1;
      for(stept=step+1;stept<=m1.num_row();stept++)
      {
        temp+=*mit2*(*sp);
        mit2+=m2.num_col();
        sp+=stept;
      }
      *(mir++)=temp;
    }
  }
  return mret;
}

G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
{
  G4ErrorMatrix mret(m1.num_row(),m1.num_row());
  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
  G4int step1,stept1,step2,stept2;
  G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
  G4double temp;
  G4ErrorMatrixIter mr = mret.m.begin();
  for(step1=1,snp1=m1.m.begin();step1<=m1.num_row();snp1+=step1++)
  {
    for(step2=1,snp2=m2.m.begin();step2<=m2.num_row();)
    {
      sp1=snp1;
      sp2=snp2;
      snp2+=step2;
      temp=0;
      if(step1<step2)
      {
        while(sp1<snp1+step1)
          { temp+=(*(sp1++))*(*(sp2++)); }
        sp1+=step1-1;
        for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
          { temp+=(*sp1)*(*(sp2++)); }
        sp2+=step2-1;
        for(stept2=++step2;stept2<=m2.num_row();sp1+=stept1++,sp2+=stept2++)
          { temp+=(*sp1)*(*sp2); }
      }
      else
      {
        while(sp2<snp2)
          { temp+=(*(sp1++))*(*(sp2++)); }
        sp2+=step2-1;
        for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
          { temp+=(*(sp1++))*(*sp2); }
        sp1+=step1-1;
        for(stept1=step1+1;stept1<=m1.num_row();sp1+=stept1++,sp2+=stept2++)
          { temp+=(*sp1)*(*sp2); }
      }
      *(mr++)=temp;
    }
  }
  return mret;
}

/* -----------------------------------------------------------------------
   This section contains the assignment and inplace operators =,+=,-=,*=,/=.
   ----------------------------------------------------------------------- */

G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorSymMatrix &m2)
{
  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
  G4int n = num_col();
  G4ErrorMatrixConstIter sjk = m2.m.begin();
  G4ErrorMatrixIter m1j = m.begin();
  G4ErrorMatrixIter mj = m.begin();
  // j >= k
  for(G4int j=1;j<=num_row();j++)
  {
    G4ErrorMatrixIter mjk = mj;
    G4ErrorMatrixIter mkj = m1j;
    for(G4int k=1;k<=j;k++)
    {
      *(mjk++) += *sjk;
      if(j!=k) *mkj += *sjk;
      sjk++;
      mkj += n;
    }
    mj += n;
    m1j++;
  }
  return (*this);
}

G4ErrorSymMatrix & G4ErrorSymMatrix::operator+=(const G4ErrorSymMatrix &m2)
{
  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
  SIMPLE_BOP(+=)
  return (*this);
}

G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorSymMatrix &m2)
{
  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
  G4int n = num_col();
  G4ErrorMatrixConstIter sjk = m2.m.begin();
  G4ErrorMatrixIter m1j = m.begin();
  G4ErrorMatrixIter mj = m.begin();
  // j >= k
  for(G4int j=1;j<=num_row();j++)
  {
    G4ErrorMatrixIter mjk = mj;
    G4ErrorMatrixIter mkj = m1j;
    for(G4int k=1;k<=j;k++)
    {
      *(mjk++) -= *sjk;
      if(j!=k) *mkj -= *sjk;
      sjk++;
      mkj += n;
    }
    mj += n;
    m1j++;
  }
  return (*this);
}

G4ErrorSymMatrix & G4ErrorSymMatrix::operator-=(const G4ErrorSymMatrix &m2)
{
  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
  SIMPLE_BOP(-=)
  return (*this);
}

G4ErrorSymMatrix & G4ErrorSymMatrix::operator/=(G4double t)
{
  SIMPLE_UOP(/=)
  return (*this);
}

G4ErrorSymMatrix & G4ErrorSymMatrix::operator*=(G4double t)
{
  SIMPLE_UOP(*=)
  return (*this);
}

G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorSymMatrix &m1)
{
   if(m1.nrow*m1.nrow != size)
   {
      size = m1.nrow * m1.nrow;
      m.resize(size);
   }
   nrow = m1.nrow;
   ncol = m1.nrow;
   G4int n = ncol;
   G4ErrorMatrixConstIter sjk = m1.m.begin();
   G4ErrorMatrixIter m1j = m.begin();
   G4ErrorMatrixIter mj = m.begin();
   // j >= k
   for(G4int j=1;j<=num_row();j++)
   {
      G4ErrorMatrixIter mjk = mj;
      G4ErrorMatrixIter mkj = m1j;
      for(G4int k=1;k<=j;k++)
      {
         *(mjk++) = *sjk;
         if(j!=k) *mkj = *sjk;
         sjk++;
         mkj += n;
      }
      mj += n;
      m1j++;
   }
   return (*this);
}

G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &m1)
{
   if(m1.nrow != nrow)
   {
      nrow = m1.nrow;
      size = m1.size;
      m.resize(size);
   }
   m = m1.m;
   return (*this);
}

// Print the Matrix.

std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q)
{
  s << G4endl;

  // Fixed format needs 3 extra characters for field,
  // while scientific needs 7

  G4int width;
  if(s.flags() & std::ios::fixed)
  {
    width = s.precision()+3;
  }
  else
  {
    width = s.precision()+7;
  }
  for(G4int irow = 1; irow<= q.num_row(); irow++)
  {
    for(G4int icol = 1; icol <= q.num_col(); icol++)
    {
      s.width(width);
      s << q(irow,icol) << " ";
    }
    s << G4endl;
  }
  return s;
}

G4ErrorSymMatrix G4ErrorSymMatrix::
apply(G4double (*f)(G4double, G4int, G4int)) const
{
  G4ErrorSymMatrix mret(num_row());
  G4ErrorMatrixConstIter a = m.begin();
  G4ErrorMatrixIter b = mret.m.begin();
  for(G4int ir=1;ir<=num_row();ir++)
  {
    for(G4int ic=1;ic<=ir;ic++)
    {
      *(b++) = (*f)(*(a++), ir, ic);
    }
  }
  return mret;
}

void G4ErrorSymMatrix::assign (const G4ErrorMatrix &m1)
{
   if(m1.nrow != nrow)
   {
      nrow = m1.nrow;
      size = nrow * (nrow+1) / 2;
      m.resize(size);
   }
   G4ErrorMatrixConstIter a = m1.m.begin();
   G4ErrorMatrixIter b = m.begin();
   for(G4int r=1;r<=nrow;r++)
   {
      G4ErrorMatrixConstIter d = a;
      for(G4int c=1;c<=r;c++)
      {
         *(b++) = *(d++);
      }
      a += nrow;
   }
}

G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorMatrix &m1) const
{
  G4ErrorSymMatrix mret(m1.num_row());
  G4ErrorMatrix temp = m1*(*this);

// If m1*(*this) has correct dimensions, then so will the m1.T multiplication.
// So there is no need to check dimensions again.

  G4int n = m1.num_col();
  G4ErrorMatrixIter mr = mret.m.begin();
  G4ErrorMatrixIter tempr1 = temp.m.begin();
  for(G4int r=1;r<=mret.num_row();r++)
  {
    G4ErrorMatrixConstIter m1c1 = m1.m.begin();
    for(G4int c=1;c<=r;c++)
    {
      G4double tmp = 0.0;
      G4ErrorMatrixIter tempri = tempr1;
      G4ErrorMatrixConstIter m1ci = m1c1;
      for(G4int i=1;i<=m1.num_col();i++)
      {
        tmp+=(*(tempri++))*(*(m1ci++));
      }
      *(mr++) = tmp;
      m1c1 += n;
    }
    tempr1 += n;
  }
  return mret;
}

G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorSymMatrix &m1) const
{
  G4ErrorSymMatrix mret(m1.num_row());
  G4ErrorMatrix temp = m1*(*this);
  G4int n = m1.num_col();
  G4ErrorMatrixIter mr = mret.m.begin();
  G4ErrorMatrixIter tempr1 = temp.m.begin();
  for(G4int r=1;r<=mret.num_row();r++)
  {
    G4ErrorMatrixConstIter m1c1 = m1.m.begin();
    G4int c;
    for(c=1;c<=r;c++)
    {
      G4double tmp = 0.0;
      G4ErrorMatrixIter tempri = tempr1;
      G4ErrorMatrixConstIter m1ci = m1c1;
      G4int i;
      for(i=1;i<c;i++)
      {
        tmp+=(*(tempri++))*(*(m1ci++));
      }
      for(i=c;i<=m1.num_col();i++)
      {
        tmp+=(*(tempri++))*(*(m1ci));
        m1ci += i;
      }
      *(mr++) = tmp;
      m1c1 += c;
    }
    tempr1 += n;
  }
  return mret;
}

G4ErrorSymMatrix G4ErrorSymMatrix::similarityT(const G4ErrorMatrix &m1) const
{
  G4ErrorSymMatrix mret(m1.num_col());
  G4ErrorMatrix temp = (*this)*m1;
  G4int n = m1.num_col();
  G4ErrorMatrixIter mrc = mret.m.begin();
  G4ErrorMatrixIter temp1r = temp.m.begin();
  for(G4int r=1;r<=mret.num_row();r++)
  {
    G4ErrorMatrixConstIter m11c = m1.m.begin();
    for(G4int c=1;c<=r;c++)
    {
      G4double tmp = 0.0;
      G4ErrorMatrixIter tempir = temp1r;
      G4ErrorMatrixConstIter m1ic = m11c;
      for(G4int i=1;i<=m1.num_row();i++)
      {
        tmp+=(*(tempir))*(*(m1ic));
        tempir += n;
        m1ic += n;
      }
      *(mrc++) = tmp;
      m11c++;
    }
    temp1r++;
  }
  return mret;
}

void G4ErrorSymMatrix::invert(G4int &ifail)
{  
  ifail = 0;

  switch(nrow)
  {
  case 3:
    {
      G4double det, temp;
      G4double t1, t2, t3;
      G4double c11,c12,c13,c22,c23,c33;
      c11 = (*(m.begin()+2)) * (*(m.begin()+5))
          - (*(m.begin()+4)) * (*(m.begin()+4));
      c12 = (*(m.begin()+4)) * (*(m.begin()+3))
          - (*(m.begin()+1)) * (*(m.begin()+5));
      c13 = (*(m.begin()+1)) * (*(m.begin()+4))
          - (*(m.begin()+2)) * (*(m.begin()+3));
      c22 = (*(m.begin()+5)) * (*m.begin())
          - (*(m.begin()+3)) * (*(m.begin()+3));
      c23 = (*(m.begin()+3)) * (*(m.begin()+1))
          - (*(m.begin()+4)) * (*m.begin());
      c33 = (*m.begin()) * (*(m.begin()+2))
          - (*(m.begin()+1)) * (*(m.begin()+1));
      t1 = std::fabs(*m.begin());
      t2 = std::fabs(*(m.begin()+1));
      t3 = std::fabs(*(m.begin()+3));
      if (t1 >= t2)
      {
        if (t3 >= t1)
        {
          temp = *(m.begin()+3);
          det = c23*c12-c22*c13;
        }
        else
        {
          temp = *m.begin();
          det = c22*c33-c23*c23;
        }
      }
      else if (t3 >= t2)
      {
        temp = *(m.begin()+3);
        det = c23*c12-c22*c13;
      }
      else
      {
        temp = *(m.begin()+1);
        det = c13*c23-c12*c33;
      }
      if (det==0)
      {
        ifail = 1;
        return;
      }
      {
        G4double s = temp/det;
        G4ErrorMatrixIter mm = m.begin();
        *(mm++) = s*c11;
        *(mm++) = s*c12;
        *(mm++) = s*c22;
        *(mm++) = s*c13;
        *(mm++) = s*c23;
        *(mm) = s*c33;
      }
    }
    break;
 case 2:
    {
      G4double det, temp, s;
      det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
      if (det==0)
      {
        ifail = 1;
        return;
      }
      s = 1.0/det;
      *(m.begin()+1) *= -s;
      temp = s*(*(m.begin()+2));
      *(m.begin()+2) = s*(*m.begin());
      *m.begin() = temp;
      break;
    }
 case 1:
    {
      if ((*m.begin())==0)
      {
        ifail = 1;
        return;
      }
      *m.begin() = 1.0/(*m.begin());
      break;
    }
 case 5:
    {
      invert5(ifail);
      return;
    }
 case 6:
    {
      invert6(ifail);
      return;
    }
 case 4:
    {
      invert4(ifail);
      return;
    }
 default:
    {
      invertBunchKaufman(ifail);
      return;
    }
  }
  return; // inversion successful
}

G4double G4ErrorSymMatrix::determinant() const
{
  static const G4int max_array = 20;

  // ir must point to an array which is ***1 longer than*** nrow

  static std::vector<G4int> ir_vec (max_array+1); 
  if (ir_vec.size() <= static_cast<unsigned int>(nrow))
  {
    ir_vec.resize(nrow+1);
  }
  G4int * ir = &ir_vec[0];   

  G4double det;
  G4ErrorMatrix mt(*this);
  G4int i = mt.dfact_matrix(det, ir);
  if(i==0) { return det; }
  return 0.0;
}

G4double G4ErrorSymMatrix::trace() const
{
   G4double t = 0.0;
   for (G4int i=0; i<nrow; i++) 
     { t += *(m.begin() + (i+3)*i/2); }
   return t;
}

void G4ErrorSymMatrix::invertBunchKaufman(G4int &ifail)
{
  // Bunch-Kaufman diagonal pivoting method
  // It is decribed in J.R. Bunch, L. Kaufman (1977). 
  // "Some Stable Methods for Calculating Inertia and Solving Symmetric 
  // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub, 
  // Charles F. van Loan, "Matrix Computations" (the second edition 
  // has a bug.) and implemented in "lapack"
  // Mario Stanke, 09/97

  G4int i, j, k, s;
  G4int pivrow;

  // Establish the two working-space arrays needed:  x and piv are
  // used as pointers to arrays of doubles and ints respectively, each
  // of length nrow.  We do not want to reallocate each time through
  // unless the size needs to grow.  We do not want to leak memory, even
  // by having a new without a delete that is only done once.
  
  static const G4int max_array = 25;
  static std::vector<G4double> xvec (max_array);
  static std::vector<G4int>    pivv (max_array);
  typedef std::vector<G4int>::iterator pivIter; 
  if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
  if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
    // Note - resize should do  nothing if the size is already larger than nrow,
    //        but on VC++ there are indications that it does so we check.
    // Note - the data elements in a vector are guaranteed to be contiguous,
    //        so x[i] and piv[i] are optimally fast.
  G4ErrorMatrixIter   x   = xvec.begin();
    // x[i] is used as helper storage, needs to have at least size nrow.
  pivIter piv = pivv.begin();
    // piv[i] is used to store details of exchanges
      
  G4double temp1, temp2;
  G4ErrorMatrixIter ip, mjj, iq;
  G4double lambda, sigma;
  const G4double alpha = .6404; // = (1+sqrt(17))/8
  const G4double epsilon = 32*DBL_EPSILON;
    // whenever a sum of two doubles is below or equal to epsilon
    // it is set to zero.
    // this constant could be set to zero but then the algorithm
    // doesn't neccessarily detect that a matrix is singular
  
  for (i = 0; i < nrow; i++)
  {
    piv[i] = i+1;
  }
      
  ifail = 0;
      
  // compute the factorization P*A*P^T = L * D * L^T 
  // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
  // L and D^-1 are stored in A = *this, P is stored in piv[]
        
  for (j=1; j < nrow; j+=s)  // main loop over columns
  {
          mjj = m.begin() + j*(j-1)/2 + j-1;
          lambda = 0;           // compute lambda = max of A(j+1:n,j)
          pivrow = j+1;
          ip = m.begin() + (j+1)*j/2 + j-1;
          for (i=j+1; i <= nrow ; ip += i++)
          {
            if (std::fabs(*ip) > lambda)
              {
                lambda = std::fabs(*ip);
                pivrow = i;
              }
          }
          if (lambda == 0 )
            {
              if (*mjj == 0)
                {
                  ifail = 1;
                  return;
                }
              s=1;
              *mjj = 1./ *mjj;
            }
          else
            {
              if (std::fabs(*mjj) >= lambda*alpha)
                {
                  s=1;
                  pivrow=j;
                }
              else
                {
                  sigma = 0;  // compute sigma = max A(pivrow, j:pivrow-1)
                  ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
                  for (k=j; k < pivrow; k++)
                    {
                      if (std::fabs(*ip) > sigma)
                        sigma = std::fabs(*ip);
                      ip++;
                    }
                  if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
                    {
                      s=1;
                      pivrow = j;
                    }
                  else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1)) 
                                >= alpha * sigma)
                    { s=1; }
                  else
                    { s=2; }
                }
              if (pivrow == j)  // no permutation neccessary
                {
                  piv[j-1] = pivrow;
                  if (*mjj == 0)
                    {
                      ifail=1;
                      return;
                    }
                  temp2 = *mjj = 1./ *mjj; // invert D(j,j)
                  
                  // update A(j+1:n, j+1,n)
                  for (i=j+1; i <= nrow; i++)
                    {
                      temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
                      ip = m.begin()+i*(i-1)/2+j;
                      for (k=j+1; k<=i; k++)
                        {
                          *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
                          if (std::fabs(*ip) <= epsilon)
                            { *ip=0; }
                          ip++;
                        }
                    }
                  // update L 
                  ip = m.begin() + (j+1)*j/2 + j-1; 
                  for (i=j+1; i <= nrow; ip += i++)
                  {
                    *ip *= temp2;
                  }
                }
              else if (s==1) // 1x1 pivot 
                {
                  piv[j-1] = pivrow;
                  
                  // interchange rows and columns j and pivrow in
                  // submatrix (j:n,j:n)
                  ip = m.begin() + pivrow*(pivrow-1)/2 + j;
                  for (i=j+1; i < pivrow; i++, ip++)
                    {
                      temp1 = *(m.begin() + i*(i-1)/2 + j-1);
                      *(m.begin() + i*(i-1)/2 + j-1)= *ip;
                      *ip = temp1;
                    }
                  temp1 = *mjj;
                  *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
                  *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
                  ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
                  iq = ip + pivrow-j;
                  for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
                    {
                      temp1 = *iq;
                      *iq = *ip;
                      *ip = temp1;
                    }
                  
                  if (*mjj == 0)
                    {
                      ifail = 1;
                      return;
                    }
                  temp2 = *mjj = 1./ *mjj; // invert D(j,j)
                  
                  // update A(j+1:n, j+1:n)
                  for (i = j+1; i <= nrow; i++)
                    {
                      temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
                      ip = m.begin()+i*(i-1)/2+j;
                      for (k=j+1; k<=i; k++)
                        {
                          *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
                          if (std::fabs(*ip) <= epsilon)
                            { *ip=0; }
                          ip++;
                        }
                    }
                  // update L
                  ip = m.begin() + (j+1)*j/2 + j-1;
                  for (i=j+1; i<=nrow; ip += i++)
                  {
                    *ip *= temp2;
                  }
                }
              else // s=2, ie use a 2x2 pivot
                {
                  piv[j-1] = -pivrow;
                  piv[j] = 0; // that means this is the second row of a 2x2 pivot
                  
                  if (j+1 != pivrow) 
                    {
                      // interchange rows and columns j+1 and pivrow in
                      // submatrix (j:n,j:n) 
                      ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
                      for (i=j+2; i < pivrow; i++, ip++)
                        {
                          temp1 = *(m.begin() + i*(i-1)/2 + j);
                          *(m.begin() + i*(i-1)/2 + j) = *ip;
                          *ip = temp1;
                        }
                      temp1 = *(mjj + j + 1);
                      *(mjj + j + 1) = 
                        *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
                      *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
                      temp1 = *(mjj + j);
                      *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
                      *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
                      ip = m.begin() + (pivrow+1)*pivrow/2 + j;
                      iq = ip + pivrow-(j+1);
                      for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
                        {
                          temp1 = *iq;
                          *iq = *ip;
                          *ip = temp1;
                        }
                    } 
                  // invert D(j:j+1,j:j+1)
                  temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j); 
                  if (temp2 == 0)
                  {
                    G4cerr
                      << "G4ErrorSymMatrix::bunch_invert: error in pivot choice"
                      << G4endl;
                  }
                  temp2 = 1. / temp2;

                  // this quotient is guaranteed to exist by the choice 
                  // of the pivot

                  temp1 = *mjj;
                  *mjj = *(mjj + j + 1) * temp2;
                  *(mjj + j + 1) = temp1 * temp2;
                  *(mjj + j) = - *(mjj + j) * temp2;
                  
                  if (j < nrow-1) // otherwise do nothing
                    {
                      // update A(j+2:n, j+2:n)
                      for (i=j+2; i <= nrow ; i++)
                        {
                          ip = m.begin() + i*(i-1)/2 + j-1;
                          temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
                          if (std::fabs(temp1 ) <= epsilon)
                            { temp1 = 0; }
                          temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
                          if (std::fabs(temp2 ) <= epsilon)
                            { temp2 = 0; }
                          for (k = j+2; k <= i ; k++)
                            {
                              ip = m.begin() + i*(i-1)/2 + k-1;
                              iq = m.begin() + k*(k-1)/2 + j-1;
                              *ip -= temp1 * *iq + temp2 * *(iq+1);
                              if (std::fabs(*ip) <= epsilon)
                                { *ip = 0; }
                            }
                        }
                      // update L
                      for (i=j+2; i <= nrow ; i++)
                        {
                          ip = m.begin() + i*(i-1)/2 + j-1;
                          temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
                          if (std::fabs(temp1) <= epsilon)
                            { temp1 = 0; }
                          *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
                          if (std::fabs(*(ip+1)) <= epsilon)
                            { *(ip+1) = 0; }
                          *ip = temp1;
                        }
                    }
                }
            }
  } // end of main loop over columns

  if (j == nrow) // the the last pivot is 1x1
  {
    mjj = m.begin() + j*(j-1)/2 + j-1;
    if (*mjj == 0)
    {
      ifail = 1;
      return;
    }
    else
    {
      *mjj = 1. / *mjj;
    }
  } // end of last pivot code

  // computing the inverse from the factorization
         
  for (j = nrow ; j >= 1 ; j -= s) // loop over columns
  {
          mjj = m.begin() + j*(j-1)/2 + j-1;
          if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
            {
              s = 1; 
              if (j < nrow)
                {
                  ip = m.begin() + (j+1)*j/2 + j-1;
                  for (i=0; i < nrow-j; ip += 1+j+i++)
                  {
                    x[i] = *ip;
                  }
                  for (i=j+1; i<=nrow ; i++)
                  {
                    temp2=0;
                    ip = m.begin() + i*(i-1)/2 + j;
                    for (k=0; k <= i-j-1; k++)
                      { temp2 += *ip++ * x[k]; }
                    for (ip += i-1; k < nrow-j; ip += 1+j+k++) 
                      { temp2 += *ip * x[k]; }
                    *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
                  }
                  temp2 = 0;
                  ip = m.begin() + (j+1)*j/2 + j-1;
                  for (k=0; k < nrow-j; ip += 1+j+k++)
                    { temp2 += x[k] * *ip; }
                  *mjj -= temp2;
                }
            }
          else //2x2 pivot, compute columns j and j-1 of the inverse
            {
              if (piv[j-1] != 0)
                { G4cerr << "error in piv" << piv[j-1] << G4endl; }
              s=2; 
              if (j < nrow)
                {
                  ip = m.begin() + (j+1)*j/2 + j-1;
                  for (i=0; i < nrow-j; ip += 1+j+i++)
                  {
                    x[i] = *ip;
                  }
                  for (i=j+1; i<=nrow ; i++)
                  {
                    temp2 = 0;
                    ip = m.begin() + i*(i-1)/2 + j;
                    for (k=0; k <= i-j-1; k++)
                      { temp2 += *ip++ * x[k]; }
                    for (ip += i-1; k < nrow-j; ip += 1+j+k++)
                      { temp2 += *ip * x[k]; }
                    *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
                  }    
                  temp2 = 0;
                  ip = m.begin() + (j+1)*j/2 + j-1;
                  for (k=0; k < nrow-j; ip += 1+j+k++)
                    { temp2 += x[k] * *ip; }
                  *mjj -= temp2;
                  temp2 = 0;
                  ip = m.begin() + (j+1)*j/2 + j-2;
                  for (i=j+1; i <= nrow; ip += i++)
                    { temp2 += *ip * *(ip+1); }
                  *(mjj-1) -= temp2;
                  ip = m.begin() + (j+1)*j/2 + j-2;
                  for (i=0; i < nrow-j; ip += 1+j+i++)
                  {
                    x[i] = *ip;
                  }
                  for (i=j+1; i <= nrow ; i++)
                  {
                    temp2 = 0;
                    ip = m.begin() + i*(i-1)/2 + j;
                    for (k=0; k <= i-j-1; k++)
                      { temp2 += *ip++ * x[k]; }
                    for (ip += i-1; k < nrow-j; ip += 1+j+k++)
                      { temp2 += *ip * x[k]; }
                    *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
                  }
                  temp2 = 0;
                  ip = m.begin() + (j+1)*j/2 + j-2;
                  for (k=0; k < nrow-j; ip += 1+j+k++)
                    { temp2 += x[k] * *ip; }
                  *(mjj-j) -= temp2;
                }
            }  
          
          // interchange rows and columns j and piv[j-1] 
          // or rows and columns j and -piv[j-2]
          
          pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
          ip = m.begin() + pivrow*(pivrow-1)/2 + j;
          for (i=j+1;i < pivrow; i++, ip++)
            {
              temp1 = *(m.begin() + i*(i-1)/2 + j-1);
              *(m.begin() + i*(i-1)/2 + j-1) = *ip;
              *ip = temp1;
            }
          temp1 = *mjj;
          *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
          *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
          if (s==2)
            {
              temp1 = *(mjj-1);
              *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
              *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
            }
          
          ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;  // &A(i,j)
          iq = ip + pivrow-j;
          for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
            {
              temp1 = *iq;
              *iq = *ip;
              *ip = temp1;
            } 
  } // end of loop over columns (in computing inverse from factorization)

  return; // inversion successful
}

G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;

// Aij are indices for a 6x6 symmetric matrix.
//     The indices for 5x5 or 4x4 symmetric matrices are the same, 
//     ignoring all combinations with an index which is inapplicable.

#define A00 0
#define A01 1
#define A02 3
#define A03 6
#define A04 10
#define A05 15

#define A10 1
#define A11 2
#define A12 4
#define A13 7
#define A14 11
#define A15 16

#define A20 3
#define A21 4
#define A22 5
#define A23 8
#define A24 12
#define A25 17

#define A30 6
#define A31 7
#define A32 8
#define A33 9
#define A34 13
#define A35 18

#define A40 10
#define A41 11
#define A42 12
#define A43 13
#define A44 14
#define A45 19

#define A50 15
#define A51 16
#define A52 17
#define A53 18
#define A54 19
#define A55 20

void G4ErrorSymMatrix::invert5(G4int & ifail)
{ 
      if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
      {
        invertCholesky5(ifail);
        posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
        if (ifail!=0)    // Cholesky failed -- invert using Haywood
        {
          invertHaywood5(ifail);      
        }
      }
      else
      {
        if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
        {
          invertCholesky5(ifail);
          posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
          if (ifail!=0)    // Cholesky failed -- invert using Haywood
          {
            invertHaywood5(ifail);      
            adjustment5x5 = 0;
          }
        }
        else
        {
          invertHaywood5(ifail);      
          adjustment5x5 += CHOLESKY_CREEP_5x5;
        }
      }
      return;
}

void G4ErrorSymMatrix::invert6(G4int & ifail)
{
      if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
      {
        invertCholesky6(ifail);
        posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
        if (ifail!=0)    // Cholesky failed -- invert using Haywood
        {
          invertHaywood6(ifail);      
        }
      }
      else
      {
        if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
        {
          invertCholesky6(ifail);
          posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
          if (ifail!=0)    // Cholesky failed -- invert using Haywood
          {
            invertHaywood6(ifail);      
            adjustment6x6 = 0;
          }
        }
        else
        {
          invertHaywood6(ifail);      
          adjustment6x6 += CHOLESKY_CREEP_6x6;
        }
      }
      return;
}

void G4ErrorSymMatrix::invertHaywood5  (G4int & ifail)
{
  ifail = 0;

  // Find all NECESSARY 2x2 dets:  (25 of them)

  G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
  G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
  G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
  G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
  G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
  G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
  G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
  G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
  G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
  G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
  G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
  G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
  G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
  G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
  G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
  G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
  G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
  G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
  G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
  G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
  G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
  G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
  G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
  G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
  G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];

  // Find all NECESSARY 3x3 dets:   (30 of them)

  G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 
                                + m[A12]*Det2_23_01;
  G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 
                                + m[A13]*Det2_23_01;
  G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 
                                + m[A13]*Det2_23_02;
  G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 
                                + m[A13]*Det2_23_12;
  G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02 
                                + m[A12]*Det2_24_01;
  G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03 
                                + m[A13]*Det2_24_01;
  G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04 
                                + m[A14]*Det2_24_01;
  G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03 
                                + m[A13]*Det2_24_02;
  G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04 
                                + m[A14]*Det2_24_02;
  G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13 
                                + m[A13]*Det2_24_12;
  G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14 
                                + m[A14]*Det2_24_12;
  G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02 
                                + m[A12]*Det2_34_01;
  G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03 
                                + m[A13]*Det2_34_01;
  G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04 
                                + m[A14]*Det2_34_01;
  G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03 
                                + m[A13]*Det2_34_02;
  G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04 
                                + m[A14]*Det2_34_02;
  G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04 
                                + m[A14]*Det2_34_03;
  G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13 
                                + m[A13]*Det2_34_12;
  G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14 
                                + m[A14]*Det2_34_12;
  G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14 
                                + m[A14]*Det2_34_13;
  G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
                                + m[A22]*Det2_34_01;
  G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
                                + m[A23]*Det2_34_01;
  G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
                                + m[A24]*Det2_34_01;
  G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
                                + m[A23]*Det2_34_02;
  G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
                                + m[A24]*Det2_34_02;
  G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
                                + m[A24]*Det2_34_03;
  G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
                                + m[A23]*Det2_34_12;
  G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
                                + m[A24]*Det2_34_12;
  G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
                                + m[A24]*Det2_34_13;
  G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
                                + m[A24]*Det2_34_23;

  // Find all NECESSARY 4x4 dets:   (15 of them)

  G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023 
                                + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
  G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023 
                                + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
  G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024 
                                + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
  G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023 
                                + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
  G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024 
                                + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
  G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034 
                                + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
  G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023 
                                + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
  G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024 
                                + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
  G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034 
                                + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
  G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034 
                                + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
  G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
                                + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
  G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
                                + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
  G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
                                + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
  G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
                                + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
  G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
                                + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;

  // Find the 5x5 det:

  G4double det =  m[A00]*Det4_1234_1234 
                - m[A01]*Det4_1234_0234 
                + m[A02]*Det4_1234_0134 
                - m[A03]*Det4_1234_0124 
                + m[A04]*Det4_1234_0123;

  if ( det == 0 )
  {  
    ifail = 1;
    return;
  } 

  G4double oneOverDet = 1.0/det;
  G4double mn1OverDet = - oneOverDet;

  m[A00] =  Det4_1234_1234 * oneOverDet;
  m[A01] =  Det4_1234_0234 * mn1OverDet;
  m[A02] =  Det4_1234_0134 * oneOverDet;
  m[A03] =  Det4_1234_0124 * mn1OverDet;
  m[A04] =  Det4_1234_0123 * oneOverDet;

  m[A11] =  Det4_0234_0234 * oneOverDet;
  m[A12] =  Det4_0234_0134 * mn1OverDet;
  m[A13] =  Det4_0234_0124 * oneOverDet;
  m[A14] =  Det4_0234_0123 * mn1OverDet;

  m[A22] =  Det4_0134_0134 * oneOverDet;
  m[A23] =  Det4_0134_0124 * mn1OverDet;
  m[A24] =  Det4_0134_0123 * oneOverDet;

  m[A33] =  Det4_0124_0124 * oneOverDet;
  m[A34] =  Det4_0124_0123 * mn1OverDet;

  m[A44] =  Det4_0123_0123 * oneOverDet;

  return;
}

void G4ErrorSymMatrix::invertHaywood6  (G4int & ifail)
{
  ifail = 0;

  // Find all NECESSARY 2x2 dets:  (39 of them)

  G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
  G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
  G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
  G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
  G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
  G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
  G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
  G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
  G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
  G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
  G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
  G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
  G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
  G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
  G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
  G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
  G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
  G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
  G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
  G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
  G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
  G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
  G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
  G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
  G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
  G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
  G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
  G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
  G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
  G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
  G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
  G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
  G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
  G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
  G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
  G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
  G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
  G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
  G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];

  // Find all NECESSARY 3x3 dets:  (65 of them)

  G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
                                                + m[A22]*Det2_34_01;
  G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
                                                + m[A23]*Det2_34_01;
  G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
                                                + m[A24]*Det2_34_01;
  G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
                                                + m[A23]*Det2_34_02;
  G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
                                                + m[A24]*Det2_34_02;
  G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
                                                + m[A24]*Det2_34_03;
  G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
                                                + m[A23]*Det2_34_12;
  G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
                                                + m[A24]*Det2_34_12;
  G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
                                                + m[A24]*Det2_34_13;
  G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
                                                + m[A24]*Det2_34_23;
  G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 
                                                + m[A22]*Det2_35_01;
  G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 
                                                + m[A23]*Det2_35_01;
  G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 
                                                + m[A24]*Det2_35_01;
  G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 
                                                + m[A25]*Det2_35_01;
  G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 
                                                + m[A23]*Det2_35_02;
  G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 
                                                + m[A24]*Det2_35_02;
  G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 
                                                + m[A25]*Det2_35_02;
  G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 
                                                + m[A24]*Det2_35_03;
  G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 
                                                + m[A25]*Det2_35_03;
  G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 
                                                + m[A23]*Det2_35_12;
  G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 
                                                + m[A24]*Det2_35_12;
  G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 
                                                + m[A25]*Det2_35_12;
  G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 
                                                + m[A24]*Det2_35_13;
  G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 
                                                + m[A25]*Det2_35_13;
  G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 
                                                + m[A24]*Det2_35_23;
  G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 
                                                + m[A25]*Det2_35_23;
  G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 
                                                + m[A22]*Det2_45_01;
  G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 
                                                + m[A23]*Det2_45_01;
  G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 
                                                + m[A24]*Det2_45_01;
  G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 
                                                + m[A25]*Det2_45_01;
  G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 
                                                + m[A23]*Det2_45_02;
  G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 
                                                + m[A24]*Det2_45_02;
  G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 
                                                + m[A25]*Det2_45_02;
  G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 
                                                + m[A24]*Det2_45_03;
  G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 
                                                + m[A25]*Det2_45_03;
  G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 
                                                + m[A25]*Det2_45_04;
  G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 
                                                + m[A23]*Det2_45_12;
  G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 
                                                + m[A24]*Det2_45_12;
  G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 
                                                + m[A25]*Det2_45_12;
  G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 
                                                + m[A24]*Det2_45_13;
  G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 
                                                + m[A25]*Det2_45_13;
  G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 
                                                + m[A25]*Det2_45_14;
  G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 
                                                + m[A24]*Det2_45_23;
  G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 
                                                + m[A25]*Det2_45_23;
  G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 
                                                + m[A25]*Det2_45_24;
  G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 
                                                + m[A32]*Det2_45_01;
  G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 
                                                + m[A33]*Det2_45_01;
  G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 
                                                + m[A34]*Det2_45_01;
  G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 
                                                + m[A35]*Det2_45_01;
  G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 
                                                + m[A33]*Det2_45_02;
  G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 
                                                + m[A34]*Det2_45_02;
  G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 
                                                + m[A35]*Det2_45_02;
  G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 
                                                + m[A34]*Det2_45_03;
  G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 
                                                + m[A35]*Det2_45_03;
  G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 
                                                + m[A35]*Det2_45_04;
  G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 
                                                + m[A33]*Det2_45_12;
  G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 
                                                + m[A34]*Det2_45_12;
  G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 
                                                + m[A35]*Det2_45_12;
  G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 
                                                + m[A34]*Det2_45_13;
  G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 
                                                + m[A35]*Det2_45_13;
  G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 
                                                + m[A35]*Det2_45_14;
  G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 
                                                + m[A34]*Det2_45_23;
  G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 
                                                + m[A35]*Det2_45_23;
  G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 
                                                + m[A35]*Det2_45_24;
  G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 
                                                + m[A35]*Det2_45_34;

  // Find all NECESSARY 4x4 dets:  (55 of them)

  G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
                        + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
  G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
                        + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
  G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
                        + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
  G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
                        + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
  G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
                        + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
  G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 
                        + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
  G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 
                        + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
  G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 
                        + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
  G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 
                        + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
  G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 
                        + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
  G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 
                        + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
  G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 
                        + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
  G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 
                        + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
  G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 
                        + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
  G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 
                        + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
  G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 
                        + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
  G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 
                        + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
  G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 
                        + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
  G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 
                        + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
  G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 
                        + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
  G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 
                        + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
  G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 
                        + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
  G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 
                        + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
  G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 
                        + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
  G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 
                        + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
  G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 
                        + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
  G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 
                        + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
  G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 
                        + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
  G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 
                        + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
  G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 
                        + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
  G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 
                        + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
  G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 
                        + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
  G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 
                        + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
  G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 
                        + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
  G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 
                        + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
  G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 
                        + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
  G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 
                        + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
  G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 
                        + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
  G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 
                        + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
  G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 
                        + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
  G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 
                        + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
  G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 
                        + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
  G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 
                        + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
  G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 
                        + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
  G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 
                        + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
  G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 
                        + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
  G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 
                        + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
  G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 
                        + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
  G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 
                        + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
  G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 
                        + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
  G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 
                        + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
  G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 
                        + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
  G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 
                        + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
  G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 
                        + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
  G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 
                        + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;

  // Find all NECESSARY 5x5 dets:  (19 of them)

  G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 
    + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
  G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 
    + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
  G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 
    + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
  G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 
    + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
  G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 
    + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
  G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 
    + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
  G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 
    + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
  G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 
    + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
  G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 
    + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
  G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 
    + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
  G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 
    + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
  G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 
    + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
  G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 
    + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
  G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 
    + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
  G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 
    + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
  G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 
    + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
  G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 
    + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
  G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 
    + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
  G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 
    + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
  G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 
    + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
  G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 
    + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;

  // Find the determinant 

  G4double det =  m[A00]*Det5_12345_12345 
                - m[A01]*Det5_12345_02345 
                + m[A02]*Det5_12345_01345 
                - m[A03]*Det5_12345_01245 
                + m[A04]*Det5_12345_01235 
                - m[A05]*Det5_12345_01234;

  if ( det == 0 )
  {  
    ifail = 1;
    return;
  } 

  G4double oneOverDet = 1.0/det;
  G4double mn1OverDet = - oneOverDet;

  m[A00] =  Det5_12345_12345*oneOverDet;
  m[A01] =  Det5_12345_02345*mn1OverDet;
  m[A02] =  Det5_12345_01345*oneOverDet;
  m[A03] =  Det5_12345_01245*mn1OverDet;
  m[A04] =  Det5_12345_01235*oneOverDet;
  m[A05] =  Det5_12345_01234*mn1OverDet;

  m[A11] =  Det5_02345_02345*oneOverDet;
  m[A12] =  Det5_02345_01345*mn1OverDet;
  m[A13] =  Det5_02345_01245*oneOverDet;
  m[A14] =  Det5_02345_01235*mn1OverDet;
  m[A15] =  Det5_02345_01234*oneOverDet;

  m[A22] =  Det5_01345_01345*oneOverDet;
  m[A23] =  Det5_01345_01245*mn1OverDet;
  m[A24] =  Det5_01345_01235*oneOverDet;
  m[A25] =  Det5_01345_01234*mn1OverDet;

  m[A33] =  Det5_01245_01245*oneOverDet;
  m[A34] =  Det5_01245_01235*mn1OverDet;
  m[A35] =  Det5_01245_01234*oneOverDet;

  m[A44] =  Det5_01235_01235*oneOverDet;
  m[A45] =  Det5_01235_01234*mn1OverDet;

  m[A55] =  Det5_01234_01234*oneOverDet;

  return;
}

void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
{

  // Invert by 
  //
  // a) decomposing M = G*G^T with G lower triangular
  //    (if M is not positive definite this will fail, leaving this unchanged)
  // b) inverting G to form H
  // c) multiplying H^T * H to get M^-1.
  //
  // If the matrix is pos. def. it is inverted and 1 is returned.
  // If the matrix is not pos. def. it remains unaltered and 0 is returned.

  G4double h10;                           // below-diagonal elements of H
  G4double h20, h21;
  G4double h30, h31, h32;
  G4double h40, h41, h42, h43;

  G4double h00, h11, h22, h33, h44;       // 1/diagonal elements of G = 
                                          // diagonal elements of H

  G4double g10;                           // below-diagonal elements of G
  G4double g20, g21;
  G4double g30, g31, g32;
  G4double g40, g41, g42, g43;

  ifail = 1;  // We start by assuing we won't succeed...

  // Form G -- compute diagonal members of H directly rather than of G
  //-------

  // Scale first column by 1/sqrt(A00)

  h00 = m[A00]; 
  if (h00 <= 0) { return; }
  h00 = 1.0 / std::sqrt(h00);

  g10 = m[A10] * h00;
  g20 = m[A20] * h00;
  g30 = m[A30] * h00;
  g40 = m[A40] * h00;

  // Form G11 (actually, just h11)

  h11 = m[A11] - (g10 * g10);
  if (h11 <= 0) { return; }
  h11 = 1.0 / std::sqrt(h11);

  // Subtract inter-column column dot products from rest of column 1 and
  // scale to get column 1 of G 

  g21 = (m[A21] - (g10 * g20)) * h11;
  g31 = (m[A31] - (g10 * g30)) * h11;
  g41 = (m[A41] - (g10 * g40)) * h11;

  // Form G22 (actually, just h22)

  h22 = m[A22] - (g20 * g20) - (g21 * g21);
  if (h22 <= 0) { return; }
  h22 = 1.0 / std::sqrt(h22);

  // Subtract inter-column column dot products from rest of column 2 and
  // scale to get column 2 of G 

  g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
  g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;

  // Form G33 (actually, just h33)

  h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
  if (h33 <= 0) { return; }
  h33 = 1.0 / std::sqrt(h33);

  // Subtract inter-column column dot product from A43 and scale to get G43

  g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;

  // Finally form h44 - if this is possible inversion succeeds

  h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
  if (h44 <= 0) { return; }
  h44 = 1.0 / std::sqrt(h44);

  // Form H = 1/G -- diagonal members of H are already correct
  //-------------

  // The order here is dictated by speed considerations

  h43 = -h33 *  g43 * h44;
  h32 = -h22 *  g32 * h33;
  h42 = -h22 * (g32 * h43 + g42 * h44);
  h21 = -h11 *  g21 * h22;
  h31 = -h11 * (g21 * h32 + g31 * h33);
  h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
  h10 = -h00 *  g10 * h11;
  h20 = -h00 * (g10 * h21 + g20 * h22);
  h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
  h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);

  // Change this to its inverse = H^T*H
  //------------------------------------

  m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
  m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
  m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
  m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
  m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
  m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
  m[A03] = h30 * h33 + h40 * h43;
  m[A13] = h31 * h33 + h41 * h43;
  m[A23] = h32 * h33 + h42 * h43;
  m[A33] = h33 * h33 + h43 * h43;
  m[A04] = h40 * h44;
  m[A14] = h41 * h44;
  m[A24] = h42 * h44;
  m[A34] = h43 * h44;
  m[A44] = h44 * h44;

  ifail = 0;
  return;
}

void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
{
  // Invert by 
  //
  // a) decomposing M = G*G^T with G lower triangular
  //    (if M is not positive definite this will fail, leaving this unchanged)
  // b) inverting G to form H
  // c) multiplying H^T * H to get M^-1.
  //
  // If the matrix is pos. def. it is inverted and 1 is returned.
  // If the matrix is not pos. def. it remains unaltered and 0 is returned.

  G4double h10;                           // below-diagonal elements of H
  G4double h20, h21;
  G4double h30, h31, h32;
  G4double h40, h41, h42, h43;
  G4double h50, h51, h52, h53, h54;

  G4double h00, h11, h22, h33, h44, h55;  // 1/diagonal elements of G = 
                                          // diagonal elements of H

  G4double g10;                           // below-diagonal elements of G
  G4double g20, g21;
  G4double g30, g31, g32;
  G4double g40, g41, g42, g43;
  G4double g50, g51, g52, g53, g54;

  ifail = 1;  // We start by assuing we won't succeed...

  // Form G -- compute diagonal members of H directly rather than of G
  //-------

  // Scale first column by 1/sqrt(A00)

  h00 = m[A00]; 
  if (h00 <= 0) { return; }
  h00 = 1.0 / std::sqrt(h00);

  g10 = m[A10] * h00;
  g20 = m[A20] * h00;
  g30 = m[A30] * h00;
  g40 = m[A40] * h00;
  g50 = m[A50] * h00;

  // Form G11 (actually, just h11)

  h11 = m[A11] - (g10 * g10);
  if (h11 <= 0) { return; }
  h11 = 1.0 / std::sqrt(h11);

  // Subtract inter-column column dot products from rest of column 1 and
  // scale to get column 1 of G 

  g21 = (m[A21] - (g10 * g20)) * h11;
  g31 = (m[A31] - (g10 * g30)) * h11;
  g41 = (m[A41] - (g10 * g40)) * h11;
  g51 = (m[A51] - (g10 * g50)) * h11;

  // Form G22 (actually, just h22)

  h22 = m[A22] - (g20 * g20) - (g21 * g21);
  if (h22 <= 0) { return; }
  h22 = 1.0 / std::sqrt(h22);

  // Subtract inter-column column dot products from rest of column 2 and
  // scale to get column 2 of G 

  g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
  g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
  g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;

  // Form G33 (actually, just h33)

  h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
  if (h33 <= 0) { return; }
  h33 = 1.0 / std::sqrt(h33);

  // Subtract inter-column column dot products from rest of column 3 and
  // scale to get column 3 of G 

  g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
  g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;

  // Form G44 (actually, just h44)

  h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
  if (h44 <= 0) { return; }
  h44 = 1.0 / std::sqrt(h44);

  // Subtract inter-column column dot product from M54 and scale to get G54

  g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;

  // Finally form h55 - if this is possible inversion succeeds

  h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
  if (h55 <= 0) { return; }
  h55 = 1.0 / std::sqrt(h55);

  // Form H = 1/G -- diagonal members of H are already correct
  //-------------

  // The order here is dictated by speed considerations

  h54 = -h44 *  g54 * h55;
  h43 = -h33 *  g43 * h44;
  h53 = -h33 * (g43 * h54 + g53 * h55);
  h32 = -h22 *  g32 * h33;
  h42 = -h22 * (g32 * h43 + g42 * h44);
  h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
  h21 = -h11 *  g21 * h22;
  h31 = -h11 * (g21 * h32 + g31 * h33);
  h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
  h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
  h10 = -h00 *  g10 * h11;
  h20 = -h00 * (g10 * h21 + g20 * h22);
  h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
  h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
  h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);

  // Change this to its inverse = H^T*H
  //------------------------------------

  m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
  m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
  m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
  m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
  m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
  m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
  m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
  m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
  m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
  m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
  m[A04] = h40 * h44 + h50 * h54;
  m[A14] = h41 * h44 + h51 * h54;
  m[A24] = h42 * h44 + h52 * h54;
  m[A34] = h43 * h44 + h53 * h54;
  m[A44] = h44 * h44 + h54 * h54;
  m[A05] = h50 * h55;
  m[A15] = h51 * h55;
  m[A25] = h52 * h55;
  m[A35] = h53 * h55;
  m[A45] = h54 * h55;
  m[A55] = h55 * h55;

  ifail = 0;
  return;
}

void G4ErrorSymMatrix::invert4  (G4int & ifail)
{
  ifail = 0;

  // Find all NECESSARY 2x2 dets:  (14 of them)

  G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
  G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
  G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
  G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
  G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
  G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
  G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
  G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
  G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
  G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
  G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
  G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
  G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
  G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];

  // Find all NECESSARY 3x3 dets:   (10 of them)

  G4double Det3_012_012 = m[A00]*Det2_12_12 - m[A01]*Det2_12_02 
                                + m[A02]*Det2_12_01;
  G4double Det3_013_012 = m[A00]*Det2_13_12 - m[A01]*Det2_13_02 
                                + m[A02]*Det2_13_01;
  G4double Det3_013_013 = m[A00]*Det2_13_13 - m[A01]*Det2_13_03
                                + m[A03]*Det2_13_01;
  G4double Det3_023_012 = m[A00]*Det2_23_12 - m[A01]*Det2_23_02 
                                + m[A02]*Det2_23_01;
  G4double Det3_023_013 = m[A00]*Det2_23_13 - m[A01]*Det2_23_03
                                + m[A03]*Det2_23_01;
  G4double Det3_023_023 = m[A00]*Det2_23_23 - m[A02]*Det2_23_03
                                + m[A03]*Det2_23_02;
  G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 
                                + m[A12]*Det2_23_01;
  G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 
                                + m[A13]*Det2_23_01;
  G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 
                                + m[A13]*Det2_23_02;
  G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 
                                + m[A13]*Det2_23_12;

  // Find the 4x4 det:

  G4double det =  m[A00]*Det3_123_123 
                - m[A01]*Det3_123_023 
                + m[A02]*Det3_123_013 
                - m[A03]*Det3_123_012;

  if ( det == 0 )
  {  
    ifail = 1;
    return;
  } 

  G4double oneOverDet = 1.0/det;
  G4double mn1OverDet = - oneOverDet;

  m[A00] =  Det3_123_123 * oneOverDet;
  m[A01] =  Det3_123_023 * mn1OverDet;
  m[A02] =  Det3_123_013 * oneOverDet;
  m[A03] =  Det3_123_012 * mn1OverDet;


  m[A11] =  Det3_023_023 * oneOverDet;
  m[A12] =  Det3_023_013 * mn1OverDet;
  m[A13] =  Det3_023_012 * oneOverDet;

  m[A22] =  Det3_013_013 * oneOverDet;
  m[A23] =  Det3_013_012 * mn1OverDet;

  m[A33] =  Det3_012_012 * oneOverDet;

  return;
}

void G4ErrorSymMatrix::invertHaywood4  (G4int & ifail)
{
  invert4(ifail); // For the 4x4 case, the method we use for invert is already
                  // the Haywood method.
}