This file is indexed.

/usr/include/root/TDecompSparse.h is in libroot-math-matrix-dev 5.34.14-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
// @(#)root/matrix:$Id$
// Authors: Fons Rademakers, Eddy Offermann   Apr 2004

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT_TDecompSparse
#define ROOT_TDecompSparse

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Sparse Decomposition class                                            //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#ifndef ROOT_TDecompBase
#include "TDecompBase.h"
#endif
#ifndef ROOT_TMatrixDSparse
#include "TMatrixDSparse.h"
#endif    
#ifndef ROOT_TArrayD
#include "TArrayD.h"
#endif
#ifndef ROOT_TArrayI
#include "TArrayI.h"
#endif

// globals
const Double_t kInitTreatAsZero         = 1.0e-12;
const Double_t kInitThresholdPivoting   = 1.0e-8;
const Double_t kInitPrecision           = 1.0e-7;

// the Threshold Pivoting parameter may need to be increased during
// the algorithm if poor precision is obtained from the linear
// solves.  kThresholdPivoting indicates the largest value we are
// willing to tolerate.

const Double_t kThresholdPivotingMax    = 1.0e-2;

// the factor in the range (1,inf) by which kThresholdPivoting is
// increased when it is found to be inadequate.

const Double_t kThresholdPivotingFactor = 10.0;

class TDecompSparse : public TDecompBase
{
protected :

   Int_t     fVerbose;

   Int_t     fIcntl[31]; // integer control numbers
   Double_t  fCntl[6];   // float control numbers
   Int_t     fInfo[21];  // array used for communication between programs

   Double_t  fPrecision; // precision we demand from the linear system solver. If it isn't
                         // attained on the first solve, we use iterative refinement and
                         // possibly refactorization with a higher value of kThresholdPivoting.

   TArrayI   fIkeep;     // pivot sequence and temporary storage information
   TArrayI   fIw;
   TArrayI   fIw1;
   TArrayI   fIw2;
   Int_t     fNsteps;
   Int_t     fMaxfrt;
   TArrayD   fW;         // temporary storage for the factorization

   Double_t  fIPessimism; // amounts by which to increase allocated factorization space when
   Double_t  fRPessimism; // inadequate space is detected. fIPessimism is for array "fIw",
                          // fRPessimism is for the array "fact".

   TMatrixDSparse fA; // original matrix; needed for the iterative solving procedure
   Int_t          fNrows;
   Int_t          fNnonZeros;
   TArrayD        fFact; // size of fFact array; may be increased during the numerical factorization
                         // if the estimate obtained during the symbolic factorization proves to be inadequate.
   TArrayI        fRowFact;
   TArrayI        fColFact;

   static Int_t NonZerosUpperTriang(const TMatrixDSparse &a);
   static void  CopyUpperTriang    (const TMatrixDSparse &a,Double_t *b);

          void  InitParam();
   static void  InitPivot (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,
                           TArrayI &Aiw,TArrayI &Aikeep,TArrayI &Aiw1,Int_t &nsteps,
                           const Int_t iflag,Int_t *icntl,Double_t *cntl,Int_t *info,Double_t &ops);
   static void   Factor   (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,TArrayD &Aa,
                           TArrayI &Aiw,TArrayI &Aikeep,const Int_t nsteps,Int_t &maxfrt,
                           TArrayI &Aiw1,Int_t *icntl,Double_t *cntl,Int_t *info);
   static void   Solve    (const Int_t n,TArrayD &Aa,TArrayI &Aiw,TArrayD &Aw,const Int_t maxfrt,
                           TVectorD &b,TArrayI &Aiw1,const Int_t nsteps,Int_t *icntl,Int_t *info);
 
   static void   InitPivot_sub1 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *iw,Int_t *ipe,
                                 Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
   static void   InitPivot_sub2 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *nv,
                                 Int_t *nxt,Int_t *lst,Int_t *ipd,Int_t *flag,const Int_t iovflo,Int_t &ncmpa,
                                 const Double_t fratio);
   static void   InitPivot_sub2a(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t &ncmpa);
   static void   InitPivot_sub3 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *iw,
                                 Int_t *ipe,Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
   static void   InitPivot_sub4 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *ips,
                                 Int_t *ipv,Int_t *nv,Int_t *flag,Int_t &ncmpa);
   static void   InitPivot_sub5 (const Int_t n,Int_t *ipe,Int_t *nv,Int_t *ips,Int_t *ne,Int_t *na,Int_t *nd,
                                 Int_t &nsteps,const Int_t nemin);
   static void   InitPivot_sub6 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *na,
                                 Int_t *ne,Int_t *nd,const Int_t nsteps,Int_t *lstki,Int_t *lstkr,Int_t *iw,
                                 Int_t *info,Double_t &ops);

   static void   Factor_sub1    (const Int_t n,const Int_t nz,Int_t &nz1,Double_t *a,const Int_t la,
                                 Int_t *irn,Int_t *icn,Int_t *iw,const Int_t liw,Int_t *perm,Int_t *iw2,
                                 Int_t *icntl,Int_t *info);
   static void   Factor_sub2    (const Int_t n,const Int_t nz,Double_t *a,const Int_t la,Int_t *iw,
                                 const Int_t liw,Int_t *perm,Int_t *nstk,const Int_t nsteps,Int_t &maxfrt,
                                 Int_t *nelim,Int_t *iw2,Int_t *icntl,Double_t *cntl,Int_t *info);
   static void   Factor_sub3    (Double_t *a,Int_t *iw,Int_t &j1,Int_t &j2,const Int_t itop,const Int_t ireal,
                                 Int_t &ncmpbr,Int_t &ncmpbi);

   static void   Solve_sub1     (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
                                 const Int_t nblk,Int_t &latop,Int_t *icntl);
   static void   Solve_sub2     (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
                                 const Int_t nblk,const Int_t latop,Int_t *icntl);
   static Int_t  IDiag          (Int_t ix,Int_t iy) { return ((iy-1)*(2*ix-iy+2))/2; }

   inline Int_t IError          () { return fInfo[2]; }
   inline Int_t MinRealWorkspace() { return fInfo[5]; }
   inline Int_t MinIntWorkspace () { return fInfo[6]; }
   inline Int_t ErrorFlag       () { return fInfo[1]; }

   // Takes values in the range [0,1]. Larger values enforce greater stability in
   // the factorization as they insist on larger pivots. Smaller values preserve
   // sparsity at the cost of using smaller pivots.

   inline Double_t GetThresholdPivoting() { return fCntl[1]; }
   inline Double_t GetTreatAsZero      () { return fCntl[3]; }

   // The factorization will not accept a pivot whose absolute value is less than fCntl[3] as
   // a 1x1 pivot or as the off-diagonal in a 2x2 pivot.

   inline void     SetThresholdPivoting(Double_t piv) { fCntl[1] = piv; }
   inline void     SetTreatAsZero      (Double_t tol) { fCntl[3] = tol; }

   virtual const TMatrixDBase &GetDecompMatrix() const { MayNotUse("GetDecompMatrix()"); return fA; }

public :

   TDecompSparse();
   TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose);
   TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose);
   TDecompSparse(const TMatrixDSparse &a,Int_t verbose);
   TDecompSparse(const TDecompSparse &another);
   virtual ~TDecompSparse() {}

   inline  void     SetVerbose (Int_t v) { fVerbose = (v) ? 1 : 0;
                                            if (fVerbose) { fIcntl[1] = fIcntl[2] = 1; fIcntl[3] = 2; }
                                            else          { fIcntl[1] = fIcntl[2] = fIcntl[3] = 0; }
                                         }
   virtual Int_t    GetNrows   () const { return fA.GetNrows(); }
   virtual Int_t    GetNcols   () const { return fA.GetNcols(); }

   virtual void     SetMatrix  (const TMatrixDSparse &a);

   virtual Bool_t   Decompose  ();
   virtual Bool_t   Solve      (      TVectorD &b);
   virtual TVectorD Solve      (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
   virtual Bool_t   Solve      (      TMatrixDColumn & /*b*/)
                               { MayNotUse("Solve(TMatrixDColumn &)"); return kFALSE; }
   virtual Bool_t   TransSolve (      TVectorD &b)            { return Solve(b); }
   virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
   virtual Bool_t   TransSolve (      TMatrixDColumn & /*b*/)
                               { MayNotUse("TransSolve(TMatrixDColumn &)"); return kFALSE; }

   virtual void     Det        (Double_t &/*d1*/,Double_t &/*d2*/)
                                { MayNotUse("Det(Double_t&,Double_t&)"); }

   void Print(Option_t *opt ="") const; // *MENU*

   TDecompSparse &operator= (const TDecompSparse &source);

   ClassDef(TDecompSparse,1) // Matrix Decompositition LU
};

#endif