Re: [ROOT] matrix inversion class; comments please

From: Valeri Fine (Faine) (fine@bnl.gov)
Date: Tue Dec 12 2000 - 17:48:55 MET


> > Hi George,
> > 
> > The efficiency issue (TArrayD vs. Double_t *) will be resolved quickly later
> > this week by some profiling. I do remember from similar checks that
> > the difference was noticable, in particular the boundary checking. Anyhow,
> > I will try to get more quantitative.
> 
> If you use operator[] or At() in the inner loop, you could definitely 
> see something, particularly in a less-optimized compilation mode (a 
> good optimizer will probably be able to observe that the value of size 
> is not changing in the loop and do the check only once or twice at the 
> edge conditions). I note that while ROOT does not provide the 
> equivalent of an UncheckedAt() function for the TArrayD (would probably 
> be a nice addition), it's not all THAT hard to evade the checks, since 
> TArrayD does give you unfettered access to its internals with 
> GetArray() (1). So if you need efficiency in your inner loop, just grab 
> the array afterwards and index off that.
> 
> (1) Could ROOT at least change the overloading of this operator so that 
> GetArray() const returns a const Double_t *, not a Double_t *, as well 
> as having Double_t * GetArray()?
> 

 I do noit understanbd why we need the index check at all for example for 
the object "s" ( see below)

TInverse::GetSubMat( 
. . .
 TArrayD s(k*l);
  for (Int_t i = 0; i < k; i++)
  {
    Int_t off1 = i*l;
    Int_t off2 = (m0+i)*n+n0;
    for (Int_t j = 0; j < l; j++)
      s[off1+j] = a[off2+j];
  }

  One may do this at once and "by hand" within the "end-user" methdos:

  static Int_t   PseudoInv       (TArrayD &a,const Int_t m,const Int_t n,const Double_t frac,
                                  Int_t &nrSingular);
  static Int_t   SvDeComp        (TArrayD &a,TArrayD &b,TArrayD &x,TArrayD &r,const Int_t m,
                                  const Int_t n,const Double_t frac,Int_t &nrSingular);
  static Int_t   SvDeCompBound   (TArrayD &a,TArrayD &b,TArrayD &e,TArrayD &d,TArrayD &x,
                                  TArrayD &r,const Int_t m,const Int_t n,const Int_t l,
                                  const Double_t frac,Int_t &nrSingular);
  
 And I think one should apply one and the same  memory allocation approach across of the entire 
package.
 Myself likes
        TArrayD s(k*l);
  from  TInverse::GetSubMat( 
 and do not like:

Int_t TInverse::Inverse(TArrayD &a)
{
  Int_t    n   = TMath::Sqrt(Long_t(a.GetSize()));
  Double_t *pa = a.GetArray();

   Int_t *index = new Int_t[n]; 

   Double_t d = 0.0;
   Int_t fail = LUDecompose(pa,n,index,d);
   if (!fail) LUBacksub(pa,n,index);

   delete [] index;
   return (fail);
}

==============

I think the solution below looks more elegant: no mess with "new" / "delete"
and no performance penalty as well.

Int_t TInverse::Inverse(TArrayD &a)
{
  Int_t    n   = TMath::Sqrt(Long_t(a.GetSize()));
  TArrayI index(n)
  Double_t *pa = a.fArray;

   Double_t d = 0.0;
   Int_t fail = LUDecompose(pa,n,index.fArray,d);
   if (!fail) LUBacksub(pa,n,index);
   return (fail);
}



This archive was generated by hypermail 2b29 : Tue Jan 02 2001 - 11:50:39 MET