<commit_before>#include "Riostream.h"
#include "TMatrixD.h"
#include "TMatrixDLazy.h"
#include "TVectorD.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"

// This macro shows several ways to invert a matrix . Each  method
// is a trade-off between accuracy of the inversion and speed.
// Which method to chose depends on "how well-behaved" the matrix is.
// This is best checked through a call to Condition(), available in each
// decomposition class. A second possibilty (less preferred) would be to
// check the determinant
//
// You can run this script with different matrix sizes

void invertMatrix(Int_t msize=6)
{
  if (msize < 2 || msize > 10) {
    cout << "2 <= msize <= 10" <<endl;
    return;
  }
  cout << "--------------------------------------------------------" <<endl;
  cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<endl;
  cout << "For each inversion procedure we check the maxmimum size " <<endl;
  cout << "of the off-diagonal elements of Inv(A) * A              " <<endl;
  cout << "--------------------------------------------------------" <<endl;

  TMatrixD H_square = THilbertMatrixD(msize,msize);

//  1. InvertFast(Double_t *det=0)
//   This method is only available for the matrix class TMatrixD/F, because an
//   symmetric matrix is not necessarily symmetric after inversion .
//   It is identical to Invert() for sizes > 6 x 6 but for smaller sizes, the
//   inversion is performed according to Cramer's rule by explicitly calculating
//   all Jacobi's sub-determinants . For instance for a 6 x 6 matrix this means:
//    # of 5 x 5 determinant : 36 
//    # of 4 x 4 determinant : 75 
//    # of 3 x 3 determinant : 80 
//    # of 2 x 2 determinant : 45    (see TMatrixD/FCramerInv.cxx)
//
//    The only "quality" control in this process is to check whether the 6 x 6
//    determinant is unequal 0 . But speed gains are significant compared to Invert() ,
//    upto an order of magnitude for sizes <= 4 x 4
//
//    The inversion is done "in place", so the original matrix will be overwritten
//    If a pointer to a Double_t is supplied the determinant is calculated
//

  cout << "1. Use .InvertFast(&det)" <<endl;
  if (msize > 6)
    cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<endl;

  Double_t det1;
  TMatrixD H1 = H_square;
  H1.InvertFast(&det1);

  // Get the maximum off-diagonal matrix value . One way to do this is to set the
  // diagonal to zero .

  TMatrixD U1(H1,TMatrixD::kMult,H_square);
  TMatrixDDiag diag1(U1); diag1 = 0.0;
  const Double_t U1_max_offdiag = (U1.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U1_max_offdiag << endl;
  cout << "  Determinant          = " << det1 <<endl;

// 2. Invert(Double_t *det=0)
//   Also only available for TMatrixD/F . Again the inversion is performed in place .
//   It consists out of a sequence of calls related to the LU decomposition:
//    - The matrix is decomposed using a scheme according to Crout which involves
//      "implicit partial pivoting", see for instance Num. Recip. (we have also available
//      a decomposition scheme that does not the scaling and is therefore even slightly
//      faster but less stable)
//      With each decomposition, a tolerance has to be specified . If this tolerance
//      requirement is not met, the matrix is regarded as being singular. The value
//      passed to this decomposition, is the data member fTol of the matrix . Its
//      default value is DBL_EPSILON, which is defined as the smallest nuber so that
//      1+DBL_EPSILON > 1
//    - The last step is a standard forward/backward substitution .<commit_msg>From <NAME>: After the introduction of the Bnch-Kaufman decomposition for symmetric matrices, we have to update some comments in the invertMatrix.C tutorial,<commit_after>#include "Riostream.h"
#include "TMatrixD.h"
#include "TMatrixDLazy.h"
#include "TVectorD.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"

// This macro shows several ways to invert a matrix . Each  method
// is a trade-off between accuracy of the inversion and speed.
// Which method to chose depends on "how well-behaved" the matrix is.
// This is best checked through a call to Condition(), available in each
// decomposition class. A second possibilty (less preferred) would be to
// check the determinant
//
// You can run this script with different matrix sizes

void invertMatrix(Int_t msize=6)
{
  if (msize < 2 || msize > 10) {
    cout << "2 <= msize <= 10" <<endl;
    return;
  }
  cout << "--------------------------------------------------------" <<endl;
  cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<endl;
  cout << "For each inversion procedure we check the maxmimum size " <<endl;
  cout << "of the off-diagonal elements of Inv(A) * A              " <<endl;
  cout << "--------------------------------------------------------" <<endl;

  TMatrixD H_square = THilbertMatrixD(msize,msize);

//  1. InvertFast(Double_t *det=0)
//   It is identical to Invert() for sizes > 6 x 6 but for smaller sizes, the
//   inversion is performed according to Cramer's rule by explicitly calculating
//   all Jacobi's sub-determinants . For instance for a 6 x 6 matrix this means:
//    # of 5 x 5 determinant : 36 
//    # of 4 x 4 determinant : 75 
//    # of 3 x 3 determinant : 80 
//    # of 2 x 2 determinant : 45    (see TMatrixD/FCramerInv.cxx)
//
//    The only "quality" control in this process is to check whether the 6 x 6
//    determinant is unequal 0 . But speed gains are significant compared to Invert() ,
//    upto an order of magnitude for sizes <= 4 x 4
//
//    The inversion is done "in place", so the original matrix will be overwritten
//    If a pointer to a Double_t is supplied the determinant is calculated
//

  cout << "1. Use .InvertFast(&det)" <<endl;
  if (msize > 6)
    cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<endl;

  Double_t det1;
  TMatrixD H1 = H_square;
  H1.InvertFast(&det1);

  // Get the maximum off-diagonal matrix value . One way to do this is to set the
  // diagonal to zero .

  TMatrixD U1(H1,TMatrixD::kMult,H_square);
  TMatrixDDiag diag1(U1); diag1 = 0.0;
  const Double_t U1_max_offdiag = (U1.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U1_max_offdiag << endl;
  cout << "  Determinant          = " << det1 <<endl;

// 2. Invert(Double_t *det=0)
//   Again the inversion is performed in place .
//   It consists out of a sequence of calls to the decomposition classes . For instance
//   for the general dense matrix TMatrixD the LU decomposition is invoked:
//    - The matrix is decomposed using a scheme according to Crout which involves
//      "implicit partial pivoting", see for instance Num. Recip. (we have also available
//      a decomposition scheme that does not the scaling and is therefore even slightly
//      faster but less stable)
//      With each decomposition, a tolerance has to be specified . If this tolerance
//      requirement is not met, the matrix is regarded as being singular. The value
//      passed to this decomposition, is the data member fTol of the matrix . Its
//      default value is DBL_EPSILON, which is defined as the smallest nuber so that
//      1+DBL_EPSILON > 1
//    - The last step is a standard forward/backward substitution .
{'fraction_non_alphanumeric': 0.13183060109289618, 'fraction_numerical': 0.015027322404371584, 'mean_word_length': 2.7012133468149644, 'pattern_counts': {'":': 0, '<': 132, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7187763', 'n_tokens_mistral': 2257, 'n_tokens_neox': 2116, 'n_words': 940}