mnvert.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id: mnvert.cxx 20880 2007-11-19 11:23:41Z rdm $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005  

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
 *                                                                    *
 **********************************************************************/

#include "Minuit2/MnMatrix.h"

#include <cmath>

namespace ROOT {

   namespace Minuit2 {


/** Inverts a symmetric matrix. Matrix is first scaled to have all ones on 
    the diagonal (equivalent to change of units) but no pivoting is done 
    since matrix is positive-definite. 
 */

int mnvert(MnAlgebraicSymMatrix& a) {
  
   unsigned int nrow = a.Nrow();
   MnAlgebraicVector s(nrow);
   MnAlgebraicVector q(nrow);
   MnAlgebraicVector pp(nrow);
   
   for(unsigned int i = 0; i < nrow; i++) {
      double si = a(i,i);
      if (si < 0.) return 1;
      s(i) = 1./std::sqrt(si);
   }
   
   for(unsigned int i = 0; i < nrow; i++)
      for(unsigned int j = i; j < nrow; j++)
         a(i,j) *= (s(i)*s(j));
   
   for(unsigned i = 0; i < nrow; i++) {
      unsigned int k = i;
      if(a(k,k) == 0.) return 1;
      q(k) = 1./a(k,k);
      pp(k) = 1.;
      a(k,k) = 0.;
      unsigned int kp1 = k + 1;
      if(k != 0) {
         for(unsigned int j = 0; j < k; j++) {
            pp(j) = a(j,k);
            q(j) = a(j,k)*q(k);
            a(j,k) = 0.;
         }
      }
      if (k != nrow-1) {
         for(unsigned int j = kp1; j < nrow; j++) {
            pp(j) = a(k,j);
            q(j) = -a(k,j)*q(k);
            a(k,j) = 0.;
         }
      }
      for(unsigned int j = 0; j < nrow; j++)
         for(k = j; k < nrow; k++)
            a(j,k) += (pp(j)*q(k)); 
   }
   
   for(unsigned int j = 0; j < nrow; j++)
      for(unsigned int k = j; k < nrow; k++)
         a(j,k) *= (s(j)*s(k)); 
   
   return 0;
}

   }  // namespace Minuit2

}  // namespace ROOT

---