From 82f3a0a318ebbc6ab05d99b9853e7b3e5acc16e2 Mon Sep 17 00:00:00 2001
From: Rene Brun <Rene.Brun@cern.ch>
Date: Wed, 4 Feb 2004 17:12:44 +0000
Subject: [PATCH] From Eddy Offermann:
1. LinkDef.h : Added NormalEqn friend functions from
TDecompChol class
2. TMatrixDBase: GetDecompMatrix() returns a TMatrixD
instead of TMatrixDBase
Added Solve and TransSolve functions
that return the solution (others
replace
the input argument which is
inconvinient
when output vector size is smaller
than
input (in ncase of fit).
Also added Invert that returns matrix
instead of replacing the input
argument.
This is convinient for a generalized
(pseudo) inverse: in (m x n) out (n x
m)
3. TDecompChol: Solve/TrnasSolve had incorrect
forward/
backward substitution function ! Tests
have been added to stressLinear to
check
Added NormlEqn solver which can handle
weights in data
4. TDecompSVD: Bug in Solve function ! index was
stop-
ping one too early . Was not noticed
by
stressLinear because the linear
equation
test wasnumerically to difficult
causing accuracy tolerances to be not
strict enough (was 0.005 , is now
DBL_EPSILON on a diagonal dominant
matrix)
git-svn-id: http://root.cern.ch/svn/root/trunk@8118 27541ba8-7e3a-0410-8455-c3a389f83636
---
matrix/inc/LinkDef.h | 7 +-
matrix/inc/TDecompBase.h | 67 ++++++++--------
matrix/inc/TDecompChol.h | 20 +++--
matrix/inc/TDecompLU.h | 17 ++--
matrix/inc/TDecompQRH.h | 18 +++--
matrix/inc/TDecompSVD.h | 20 ++---
matrix/src/TDecompBase.cxx | 28 +++++--
matrix/src/TDecompChol.cxx | 154 ++++++++++++++++++++++++++-----------
matrix/src/TDecompLU.cxx | 45 ++++++++---
matrix/src/TDecompQRH.cxx | 56 ++++++++++----
matrix/src/TDecompSVD.cxx | 57 ++++++++++----
11 files changed, 332 insertions(+), 157 deletions(-)
diff --git a/matrix/inc/LinkDef.h b/matrix/inc/LinkDef.h
index e1e751e4020..31f275707eb 100644
--- a/matrix/inc/LinkDef.h
+++ b/matrix/inc/LinkDef.h
@@ -1,4 +1,4 @@
-/* @(#)root/matrix:$Name: $:$Id: LinkDef.h,v 1.9 2004/01/25 20:33:32 brun Exp $ */
+/* @(#)root/matrix:$Name: $:$Id: LinkDef.h,v 1.10 2004/01/25 23:28:44 rdm Exp $ */
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
@@ -167,4 +167,9 @@
#pragma link C++ function VerifyMatrixValue (const TMatrixDBase &,Double_t ,Int_t,Double_t);
#pragma link C++ function VerifyMatrixIdentity(const TMatrixDBase &,const TMatrixDBase &,Int_t,Double_t);
+#pragma link C++ function NormalEqn (const TMatrixD &,const TVectorD & );
+#pragma link C++ function NormalEqn (const TMatrixD &,const TVectorD &,const TVectorD &);
+#pragma link C++ function NormalEqn (const TMatrixD &,const TMatrixD & );
+#pragma link C++ function NormalEqn (const TMatrixD &,const TMatrixD &,const TVectorD &);
+
#endif
diff --git a/matrix/inc/TDecompBase.h b/matrix/inc/TDecompBase.h
index 623c6a0d29d..53112f9bbcf 100644
--- a/matrix/inc/TDecompBase.h
+++ b/matrix/inc/TDecompBase.h
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompBase.h,v 1.1 2004/01/25 20:33:32 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompBase.h,v 1.2 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -23,8 +23,8 @@
//#include <limits>
-#ifndef ROOT_TMatrixDBase
-#include "TMatrixDBase.h"
+#ifndef ROOT_TMatrixD
+#include "TMatrixD.h"
#endif
class TDecompBase : public TObject
@@ -38,7 +38,7 @@ protected :
Int_t Hager(Double_t& est,Int_t iter=5);
- virtual const TMatrixDBase &GetDecompMatrix() const = 0;
+ virtual const TMatrixD &GetDecompMatrix() const = 0;
public :
enum EMatrixDecompStat { kInit=0,kDecomposed=1,kDetermined=2,kCondition=4,kSingular=8 };
@@ -48,27 +48,30 @@ public :
TDecompBase(const TDecompBase &another);
virtual ~TDecompBase() {};
- inline Int_t GetStatus () const { return fStatus; }
- inline Double_t GetTol () const { return fTol; }
- inline Double_t GetDet1 () const { return fDet1; }
- inline Double_t GetDet2 () const { return fDet2; }
- inline Double_t GetCondition () const { return fCondition; }
- virtual Int_t GetNrows () const = 0;
- virtual Int_t GetNcols () const = 0;
- inline Double_t SetTol (Double_t tol);
-
- virtual Double_t Condition ();
- virtual void Det (Double_t &d1,Double_t &d2);
- virtual Int_t Decompose (const TMatrixDBase &a) = 0;
- virtual Bool_t Solve (TVectorD& b) = 0;
- virtual Bool_t TransSolve(TVectorD& b) = 0;
- virtual Bool_t Solve (TMatrixDColumn& b) = 0;
- virtual Bool_t TransSolve(TMatrixDColumn& b) = 0;
-
- virtual Bool_t MultiSolve(TMatrixDBase& B);
- virtual void Invert (TMatrixDBase& inv);
-
- static void DiagProd (const TVectorD &diag,Double_t tol,Double_t &d1,Double_t &d2);
+ inline Int_t GetStatus () const { return fStatus; }
+ inline Double_t GetTol () const { return fTol; }
+ inline Double_t GetDet1 () const { return fDet1; }
+ inline Double_t GetDet2 () const { return fDet2; }
+ inline Double_t GetCondition () const { return fCondition; }
+ virtual Int_t GetNrows () const = 0;
+ virtual Int_t GetNcols () const = 0;
+ inline Double_t SetTol (Double_t tol);
+
+ virtual Double_t Condition ();
+ virtual void Det (Double_t &d1,Double_t &d2);
+ virtual Int_t Decompose (const TMatrixDBase &a) = 0;
+ virtual Bool_t Solve ( TVectorD &b) = 0;
+ virtual TVectorD Solve (const TVectorD& b,Bool_t &ok) = 0;
+ virtual Bool_t Solve ( TMatrixDColumn& b) = 0;
+ virtual Bool_t TransSolve ( TVectorD &b) = 0;
+ virtual TVectorD TransSolve (const TVectorD &b,Bool_t &ok) = 0;
+ virtual Bool_t TransSolve ( TMatrixDColumn& b) = 0;
+
+ virtual Bool_t MultiSolve (TMatrixD& B);
+ virtual void Invert (TMatrixD& inv);
+ virtual TMatrixD Invert ();
+
+ static void DiagProd (const TVectorD &diag,Double_t tol,Double_t &d1,Double_t &d2);
TDecompBase &operator= (const TDecompBase &source);
@@ -83,12 +86,12 @@ Double_t TDecompBase::SetTol(Double_t newTol)
return oldTol;
}
-Bool_t DefHouseHolder (const TVectorD &vc,Int_t lp,Int_t l,Double_t &up,Double_t &b,Double_t tol=0.0);
-void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TMatrixDRow &cr);
-void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TMatrixDColumn &cc);
-void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TVectorD &cv);
-void DefGivens (Double_t v1,Double_t v2,Double_t &c,Double_t &s);
-void DefAplGivens (Double_t &v1,Double_t &v2,Double_t &c,Double_t &s);
-void ApplyGivens (Double_t &z1,Double_t &z2,Double_t c,Double_t s);
+Bool_t DefHouseHolder (const TVectorD &vc,Int_t lp,Int_t l,Double_t &up,Double_t &b,Double_t tol=0.0);
+void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TMatrixDRow &cr);
+void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TMatrixDColumn &cc);
+void ApplyHouseHolder(const TVectorD &vc,Double_t up,Double_t b,Int_t lp,Int_t l,TVectorD &cv);
+void DefGivens ( Double_t v1,Double_t v2,Double_t &c,Double_t &s);
+void DefAplGivens ( Double_t &v1,Double_t &v2,Double_t &c,Double_t &s);
+void ApplyGivens ( Double_t &z1,Double_t &z2,Double_t c,Double_t s);
#endif
diff --git a/matrix/inc/TDecompChol.h b/matrix/inc/TDecompChol.h
index 5f2cc0332fb..69ca5bc4c22 100644
--- a/matrix/inc/TDecompChol.h
+++ b/matrix/inc/TDecompChol.h
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompChol.h,v 1.1 2004/01/25 20:33:32 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompChol.h,v 1.2 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -28,7 +28,7 @@ protected :
TMatrixD fU; // decomposed matrix fU so that a = fU^T fU
- virtual const TMatrixDBase &GetDecompMatrix() const { return fU; }
+ virtual const TMatrixD &GetDecompMatrix() const { return fU; }
public :
@@ -43,12 +43,14 @@ public :
virtual Int_t GetNcols () const { return fU.GetNcols(); }
const TMatrixD &GetU () const { return fU; }
- virtual Int_t Decompose (const TMatrixDBase &a);
- virtual Bool_t Solve (TVectorD &b);
- virtual Bool_t Solve (TMatrixDColumn &b);
- virtual Bool_t TransSolve(TVectorD &b) { return Solve(b); }
- virtual Bool_t TransSolve(TMatrixDColumn &b) { return Solve(b); }
- virtual void Det (Double_t &d1,Double_t &d2);
+ virtual Int_t Decompose (const TMatrixDBase &a);
+ virtual Bool_t Solve ( TVectorD &b);
+ virtual TVectorD Solve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t Solve ( TMatrixDColumn &b);
+ 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) { return Solve(b); }
+ virtual void Det (Double_t &d1,Double_t &d2);
TDecompChol &operator= (const TDecompChol &source);
@@ -56,6 +58,8 @@ public :
};
TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b);
+TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b,const TVectorD &std);
TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &b);
+TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &B,const TVectorD &std);
#endif
diff --git a/matrix/inc/TDecompLU.h b/matrix/inc/TDecompLU.h
index a6731cba1cd..bccb3e6fdb4 100644
--- a/matrix/inc/TDecompLU.h
+++ b/matrix/inc/TDecompLU.h
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompLU.h,v 1.2 2004/01/28 07:39:18 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompLU.h,v 1.3 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -32,7 +32,7 @@ protected :
TMatrixD fLU; // decomposed matrix so that a = l u where
// l is stored lower left and u upper right side
- virtual const TMatrixDBase &GetDecompMatrix() const { return fLU; }
+ virtual const TMatrixD &GetDecompMatrix() const { return fLU; }
public :
@@ -46,11 +46,14 @@ public :
virtual Int_t GetNcols () const { return fLU.GetNcols(); }
const TMatrixD &GetLU () const { return fLU; }
- virtual Int_t Decompose (const TMatrixDBase &a);
- virtual Bool_t Solve (TVectorD &b);
- virtual Bool_t Solve (TMatrixDColumn &b);
- virtual Bool_t TransSolve(TVectorD &b);
- virtual Bool_t TransSolve(TMatrixDColumn &b);
+ virtual Int_t Decompose (const TMatrixDBase &a);
+ virtual Bool_t Solve ( TVectorD &b);
+ virtual TVectorD Solve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t Solve ( TMatrixDColumn &b);
+ virtual Bool_t TransSolve ( TVectorD &b);
+ virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t TransSolve ( TMatrixDColumn &b);
+
virtual Double_t Condition ();
virtual void Det (Double_t &d1,Double_t &d2);
diff --git a/matrix/inc/TDecompQRH.h b/matrix/inc/TDecompQRH.h
index 9c69acf0c0b..0d4a33bbcbb 100644
--- a/matrix/inc/TDecompQRH.h
+++ b/matrix/inc/TDecompQRH.h
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompQRH.h,v 1.1 2004/01/25 20:33:32 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompQRH.h,v 1.2 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -34,7 +34,7 @@ protected :
static Bool_t QRH(TMatrixD &q,TVectorD &diagR,TVectorD &up,TVectorD &w,Double_t tol);
- virtual const TMatrixDBase &GetDecompMatrix() const { return fR; }
+ virtual const TMatrixD &GetDecompMatrix() const { return fR; }
public :
@@ -50,12 +50,14 @@ public :
virtual const TVectorD &GetUp () const { return fUp; }
virtual const TVectorD &GetW () const { return fW; }
- virtual Int_t Decompose (const TMatrixDBase &a);
- virtual Bool_t Solve (TVectorD &b);
- virtual Bool_t Solve (TMatrixDColumn &b);
- virtual Bool_t TransSolve(TVectorD &b);
- virtual Bool_t TransSolve(TMatrixDColumn &b);
- virtual void Det (Double_t &d1,Double_t &d2);
+ virtual Int_t Decompose (const TMatrixDBase &a);
+ virtual Bool_t Solve ( TVectorD &b);
+ virtual TVectorD Solve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t Solve ( TMatrixDColumn &b);
+ virtual Bool_t TransSolve ( TVectorD &b);
+ virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t TransSolve ( TMatrixDColumn &b);
+ virtual void Det (Double_t &d1,Double_t &d2);
TDecompQRH &operator= (const TDecompQRH &source);
diff --git a/matrix/inc/TDecompSVD.h b/matrix/inc/TDecompSVD.h
index 86f25b2d5ae..20ed36c9a2b 100644
--- a/matrix/inc/TDecompSVD.h
+++ b/matrix/inc/TDecompSVD.h
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompSVD.h,v 1.1 2004/01/25 20:33:32 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompSVD.h,v 1.2 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -38,7 +38,7 @@ protected :
static void Diag_3 (TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l);
static void SortSingular (TMatrixD &v,TMatrixD &u,TVectorD &sDiag);
- virtual const TMatrixDBase &GetDecompMatrix() const { return fU; }
+ virtual const TMatrixD &GetDecompMatrix() const { return fU; }
public :
@@ -54,13 +54,15 @@ public :
const TMatrixD &GetV () const { return fV; }
const TVectorD &GetSig () const { return fSig; }
- virtual Int_t Decompose (const TMatrixDBase &a);
- virtual Bool_t Solve (TVectorD &b);
- virtual Bool_t Solve (TMatrixDColumn &b);
- virtual Bool_t TransSolve (TVectorD &b);
- virtual Bool_t TransSolve (TMatrixDColumn &b);
- virtual Double_t Condition ();
- virtual void Det (Double_t &d1,Double_t &d2);
+ virtual Int_t Decompose (const TMatrixDBase &a);
+ virtual Bool_t Solve ( TVectorD &b);
+ virtual TVectorD Solve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t Solve ( TMatrixDColumn &b);
+ virtual Bool_t TransSolve ( TVectorD &b);
+ virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok);
+ virtual Bool_t TransSolve ( TMatrixDColumn &b);
+ virtual Double_t Condition ();
+ virtual void Det (Double_t &d1,Double_t &d2);
TDecompSVD &operator= (const TDecompSVD &source);
diff --git a/matrix/src/TDecompBase.cxx b/matrix/src/TDecompBase.cxx
index a6472db8aeb..b8aea09fab7 100644
--- a/matrix/src/TDecompBase.cxx
+++ b/matrix/src/TDecompBase.cxx
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompBase.cxx,v 1.4 2004/01/27 08:12:26 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompBase.cxx,v 1.5 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -157,7 +157,7 @@ Int_t TDecompBase::Hager(Double_t &est,Int_t iter)
// This routine uses Hager's Convex Optimisation Algorithm.
// See Applied Numerical Linear Algebra, p139 & SIAM J Sci Stat Comp 1984 pp 311-16
- const TMatrixDBase &m = GetDecompMatrix();
+ const TMatrixD &m = GetDecompMatrix();
Assert(m.IsValid());
const Int_t n = m.GetNrows();
@@ -261,11 +261,11 @@ Double_t TDecompBase::Condition()
}
//______________________________________________________________________________
-Bool_t TDecompBase::MultiSolve(TMatrixDBase &B)
+Bool_t TDecompBase::MultiSolve(TMatrixD &B)
{
// Solve set of equations with RHS in columns of B
- const TMatrixDBase &m = GetDecompMatrix();
+ const TMatrixD &m = GetDecompMatrix();
Assert(m.IsValid() && B.IsValid());
const Int_t colLwb = B.GetColLwb();
@@ -283,9 +283,9 @@ Bool_t TDecompBase::MultiSolve(TMatrixDBase &B)
}
//______________________________________________________________________________
-void TDecompBase::Invert(TMatrixDBase &inv)
+void TDecompBase::Invert(TMatrixD &inv)
{
- // For a matrix A(m,n) the inverse is so that A * A_inv = A_inv * A = unit
+ // For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
// The user should always supply a matrix of size (m x m) !
// If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
// should be used .
@@ -301,6 +301,20 @@ void TDecompBase::Invert(TMatrixDBase &inv)
MultiSolve(inv);
}
+//______________________________________________________________________________
+TMatrixD TDecompBase::Invert()
+{
+ // For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
+ // (n x m) Ainv is returned .
+
+ TMatrixD inv(GetNrows(),GetNrows());
+ inv.UnitMatrix();
+ MultiSolve(inv);
+ inv.ResizeTo(GetNcols(),GetNrows());
+
+ return inv;
+}
+
//______________________________________________________________________________
void TDecompBase::Det(Double_t &d1,Double_t &d2)
{
@@ -309,7 +323,7 @@ void TDecompBase::Det(Double_t &d1,Double_t &d2)
fDet1 = 0.0;
fDet2 = 0.0;
} else {
- const TMatrixDBase &m = GetDecompMatrix();
+ const TMatrixD &m = GetDecompMatrix();
Assert(m.IsValid());
const TVectorD diagv = TMatrixDDiag_const(m);
DiagProd(diagv,fTol,fDet1,fDet2);
diff --git a/matrix/src/TDecompChol.cxx b/matrix/src/TDecompChol.cxx
index 5d7ff578d7e..8a099969aaf 100644
--- a/matrix/src/TDecompChol.cxx
+++ b/matrix/src/TDecompChol.cxx
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompChol.cxx,v 1.2 2004/01/27 08:12:26 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompChol.cxx,v 1.3 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -86,6 +86,7 @@ Int_t TDecompChol::Decompose(const TMatrixDBase &a)
if (irow == icol)
{
if (pU[rowOff+irow] <= 0) {
+ Error("Decompose(const TMatrixDBase &","matrix not positive definite");
fU.Invalidate();
return kFALSE;
}
@@ -112,8 +113,8 @@ const TMatrixD TDecompChol::GetMatrix() const
Bool_t TDecompChol::Solve(TVectorD &b)
{
// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
-// assumed to be in upper triang of fU. The real fTol is used to determine if
-// diagonal element is zero. The solution is returned in b.
+// assumed to be in upper triang of fU. fTol is used to determine if diagonal
+// element is zero. The solution is returned in b.
Assert(b.IsValid());
Assert(fStatus & kDecomposed);
@@ -127,83 +128,96 @@ Bool_t TDecompChol::Solve(TVectorD &b)
const Int_t n = fU.GetNrows();
- const Double_t *pA = fU.GetMatrixArray();
+ const Double_t *pU = fU.GetMatrixArray();
Double_t *pb = b.GetMatrixArray();
Int_t i;
- // step 1: Forward substitution
- for (i = n-1; i >= 0; i--) {
+ // step 1: Forward substitution on U^T
+ for (i = 0; i < n; i++) {
const Int_t off_i = i*n;
- if (pA[off_i+i] < fTol)
+ if (pU[off_i+i] < fTol)
{
- Error("Solve(TVectorD &b)","u[%d,%d]=%.4e < %.4e",i,i,pA[off_i+i],fTol);
+ Error("Solve(TVectorD &b)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
return kFALSE;
}
Double_t r = pb[i];
- for (Int_t j = n-1; j >= i; j--) {
+ for (Int_t j = 0; j < i; j++) {
const Int_t off_j = j*n;
- r -= pA[off_j+i]*pb[j];
+ r -= pU[off_j+i]*pb[j];
}
- pb[i] = r/pA[off_i+i];
+ pb[i] = r/pU[off_i+i];
}
- // step 2: Backward substitution
- for (i = 0; i < n; i++) {
+ // step 2: Backward substitution on U
+ for (i = n-1; i >= 0; i--) {
const Int_t off_i = i*n;
Double_t r = pb[i];
- for (Int_t j = 0; j < i-1; j++)
- r -= pA[off_i+j]*pb[j];
- pb[i] = r/pA[off_i+i];
+ for (Int_t j = i+1; j < n; j++)
+ r -= pU[off_i+j]*pb[j];
+ pb[i] = r/pU[off_i+i];
}
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompChol::Solve(const TVectorD &b,Bool_t &ok)
+{
+// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
+// assumed to be in upper triang of fU. fTol is used to determine if diagonal
+// element is zero.
+
+ TVectorD x = b;
+ ok = Solve(x);
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompChol::Solve(TMatrixDColumn &cb)
-{
+{
const TMatrixDBase *b = cb.GetMatrix();
Assert(b->IsValid());
Assert(fStatus & kDecomposed);
-
+
if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb())
- {
+ {
Error("Solve(TMatrixDColumn &cb","vector and matrix incompatible");
- return kFALSE;
+ return kFALSE;
}
-
+
const Int_t n = fU.GetNrows();
-
- const Double_t *pA = fU.GetMatrixArray();
+
+ const Double_t *pU = fU.GetMatrixArray();
Double_t *pcb = cb.GetPtr();
- const Int_t inc = cb.GetInc();
-
+ const Int_t inc = cb.GetInc();
+
Int_t i;
- // step 1: Forward substitution
- for (i = n-1; i >= 0; i--) {
+ // step 1: Forward substitution U^T
+ for (i = 0; i < n; i++) {
const Int_t off_i = i*n;
const Int_t off_i2 = i*inc;
- if (pA[off_i+i] < fTol)
+ if (pU[off_i+i] < fTol)
{
- Error("Solve(TMatrixDColumn &cb)","u[%d,%d]=%.4e < %.4e",i,i,pA[off_i+i],fTol);
+ Error("Solve(TMatrixDColumn &cb)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
return kFALSE;
}
Double_t r = pcb[off_i2];
- for (Int_t j = n-1; j >= i; j--) {
+ for (Int_t j = 0; j < i; j++) {
const Int_t off_j = j*n;
- r -= pA[off_j+i]*pcb[j*inc];
+ r -= pU[off_j+i]*pcb[j*inc];
}
- pcb[off_i2] = r/pA[off_i+i];
+ pcb[off_i2] = r/pU[off_i+i];
}
-
- // step 2: Backward substitution
- for (i = 0; i < n; i++) {
+
+ // step 2: Backward substitution U
+ for (i = n-1; i >= 0; i--) {
const Int_t off_i = i*n;
const Int_t off_i2 = i*inc;
Double_t r = pcb[off_i2];
- for (Int_t j = 0; j < i-1; j++)
- r -= pA[off_i+j]*pcb[j*inc];
- pcb[off_i2] = r/pA[off_i+i];
+ for (Int_t j = i+1; j < n; j++)
+ r -= pU[off_i+j]*pcb[j*inc];
+ pcb[off_i2] = r/pU[off_i+i];
}
return kTRUE;
@@ -243,27 +257,77 @@ TDecompChol &TDecompChol::operator=(const TDecompChol &source)
//______________________________________________________________________________
TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b)
{
- // Solve A . x = b for vector x where
+ // Solve min {(A . x - b)^T (A . x - b)} for vector x where
// A : (m x n) matrix, m >= n
// b : (m) vector
// x : (n) vector
- TDecompChol ch(TMatrixD(TMatrixDBase::kAtA,A));
- TVectorD x = TMatrixD(TMatrixDBase::kTransposed,A)*b;
- ch.Solve(x);
- return x;
+ TDecompChol ch(TMatrixDSym(TMatrixDBase::kAtA,A));
+ Bool_t ok;
+ return ch.Solve(TMatrixD(TMatrixDBase::kTransposed,A)*b,ok);
+}
+
+//______________________________________________________________________________
+TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b,const TVectorD &std)
+{
+ // Solve min {(A . x - b)^T W (A . x - b)} for vector x where
+ // A : (m x n) matrix, m >= n
+ // b : (m) vector
+ // x : (n) vector
+ // W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
+ // = 0 fir i != j
+
+ if (!AreCompatible(b,std)) {
+ ::Error("NormalEqn","vectors b and std are not compatible");
+ return TVectorD();
+ }
+
+ TMatrixD Aw = A;
+ TVectorD bw = b;
+ for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
+ TMatrixDRow(Aw,irow) *= 1/std(irow);
+ bw(irow) /= std(irow);
+ }
+ TDecompChol ch(TMatrixDSym(TMatrixDBase::kAtA,Aw));
+ Bool_t ok;
+ return ch.Solve(TMatrixD(TMatrixDBase::kTransposed,Aw)*bw,ok);
}
//______________________________________________________________________________
TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &B)
{
- // Solve A . x = B for matrix x where
+ // Solve min {(A . X_j - B_j)^T (A . X_j - B_j)} for each column j in
+ // B and X
// A : (m x n ) matrix, m >= n
// B : (m x nb) matrix, nb >= 1
- // x : (n x nb) matrix
+ // X : (n x nb) matrix
TDecompChol ch(TMatrixDSym(TMatrixDBase::kAtA,A));
TMatrixD X(A,TMatrixDBase::kTransposeMult,B);
ch.MultiSolve(X);
return X;
}
+
+//______________________________________________________________________________
+TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &B,const TVectorD &std)
+{
+ // Solve min {(A . X_j - B_j)^T W (A . X_j - B_j)} for each column j in
+ // B and X
+ // A : (m x n ) matrix, m >= n
+ // B : (m x nb) matrix, nb >= 1
+ // X : (n x nb) matrix
+ // W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
+ // = 0 fir i != j
+
+ TMatrixD Aw = A;
+ TMatrixD Bw = B;
+ for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
+ TMatrixDRow(Aw,irow) *= 1/std(irow);
+ TMatrixDRow(Bw,irow) *= 1/std(irow);
+ }
+
+ TDecompChol ch(TMatrixDSym(TMatrixDBase::kAtA,Aw));
+ TMatrixD X(Aw,TMatrixDBase::kTransposeMult,Bw);
+ ch.MultiSolve(X);
+ return X;
+}
diff --git a/matrix/src/TDecompLU.cxx b/matrix/src/TDecompLU.cxx
index 05e9bf51d89..25c652c0c23 100644
--- a/matrix/src/TDecompLU.cxx
+++ b/matrix/src/TDecompLU.cxx
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompLU.cxx,v 1.3 2004/01/28 07:39:18 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompLU.cxx,v 1.4 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -182,20 +182,32 @@ Bool_t TDecompLU::Solve(TVectorD &b)
}
//______________________________________________________________________________
-Bool_t TDecompLU::Solve(TMatrixDColumn &cb)
+TVectorD TDecompLU::Solve(const TVectorD &b,Bool_t &ok)
{
// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
-// been transformed. Solution returned in b.
+// been transformed.
+ TVectorD x = b;
+ ok = Solve(x);
+
+ return x;
+}
+
+//______________________________________________________________________________
+Bool_t TDecompLU::Solve(TMatrixDColumn &cb)
+{
+// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
+// been transformed. Solution returned in b.
+
const TMatrixDBase *b = cb.GetMatrix();
Assert(b->IsValid());
Assert(fStatus & kDecomposed);
-
- if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
+
+ if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
return kFALSE;
- }
-
+ }
+
const Int_t n = fLU.GetNrows();
const Double_t *pLU = fLU.GetMatrixArray();
@@ -243,6 +255,7 @@ Bool_t TDecompLU::Solve(TMatrixDColumn &cb)
return kTRUE;
}
+
//______________________________________________________________________________
Bool_t TDecompLU::TransSolve(TVectorD &b)
{
@@ -304,6 +317,18 @@ Bool_t TDecompLU::TransSolve(TVectorD &b)
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompLU::TransSolve(const TVectorD &b,Bool_t &ok)
+{
+// Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has *not*
+// been transformed.
+
+ TVectorD x = b;
+ ok = TransSolve(x);
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompLU::TransSolve(TMatrixDColumn &cb)
{
@@ -314,10 +339,10 @@ Bool_t TDecompLU::TransSolve(TMatrixDColumn &cb)
Assert(b->IsValid());
Assert(fStatus & kDecomposed);
- if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
+ if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
- return kFALSE;
- }
+ return kFALSE;
+ }
const Int_t n = fLU.GetNrows();
const Int_t lwb = fLU.GetRowLwb();
diff --git a/matrix/src/TDecompQRH.cxx b/matrix/src/TDecompQRH.cxx
index 88bb2b732eb..dc8cb2589a9 100644
--- a/matrix/src/TDecompQRH.cxx
+++ b/matrix/src/TDecompQRH.cxx
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompQRH.cxx,v 1.2 2004/01/27 08:12:26 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompQRH.cxx,v 1.3 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -177,31 +177,45 @@ Bool_t TDecompQRH::Solve(TVectorD &b)
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompQRH::Solve(const TVectorD &b,Bool_t &ok)
+{
+// Solve A x=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
+// has *not* been transformed.
+
+ TVectorD x = b;
+ ok = Solve(x);
+ if (GetNrows() != GetNcols())
+ x.ResizeTo(GetNcols());
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompQRH::Solve(TMatrixDColumn &cb)
-{
+{
const TMatrixDBase *b = cb.GetMatrix();
Assert(b->IsValid());
Assert(fStatus & kDecomposed);
if (fQ.GetNrows() != b->GetNrows() || fQ.GetRowLwb() != b->GetRowLwb())
- {
+ {
Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
return kFALSE;
- }
-
+ }
+
const Int_t nQRow = fQ.GetNrows();
const Int_t nQCol = fQ.GetNcols();
-
+
// Calculate Q^T.b
const Int_t nQ = (nQRow <= nQCol) ? nQRow-1 : nQCol;
for (Int_t k = 0; k < nQ; k++) {
const TVectorD qc_k = TMatrixDColumn_const(fQ,k);
ApplyHouseHolder(qc_k,fUp(k),fW(k),k,k+1,cb);
- }
-
+ }
+
const Int_t nRCol = fR.GetNcols();
-
+
const Double_t *pR = fR.GetMatrixArray();
Double_t *pcb = cb.GetPtr();
const Int_t inc = cb.GetInc();
@@ -277,6 +291,18 @@ Bool_t TDecompQRH::TransSolve(TVectorD &b)
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompQRH::TransSolve(const TVectorD &b,Bool_t &ok)
+{
+// Solve A^T x=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
+// has *not* been transformed.
+
+ TVectorD x = b;
+ ok = TransSolve(x);
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompQRH::TransSolve(TMatrixDColumn &cb)
{
@@ -288,18 +314,18 @@ Bool_t TDecompQRH::TransSolve(TMatrixDColumn &cb)
Error("TransSolve(TMatrixDColumn &","matrix should be square");
return kFALSE;
}
-
- if (fR.GetNrows() != b->GetNrows() || fR.GetRowLwb() != b->GetRowLwb()) {
+
+ if (fR.GetNrows() != b->GetNrows() || fR.GetRowLwb() != b->GetRowLwb()) {
Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
return kFALSE;
- }
-
+ }
+
const Double_t *pR = fR.GetMatrixArray();
Double_t *pcb = cb.GetPtr();
const Int_t inc = cb.GetInc();
-
+
const Int_t nRCol = fR.GetNcols();
-
+
// Backward substitution
for (Int_t i = 0; i < nRCol; i++) {
const Int_t off_i = i*nRCol;
diff --git a/matrix/src/TDecompSVD.cxx b/matrix/src/TDecompSVD.cxx
index 1053cbd9cb4..2ecc85a55e0 100644
--- a/matrix/src/TDecompSVD.cxx
+++ b/matrix/src/TDecompSVD.cxx
@@ -1,4 +1,4 @@
-// @(#)root/matrix:$Name: $:$Id: TDecompSVD.cxx,v 1.5 2004/01/29 16:08:54 brun Exp $
+// @(#)root/matrix:$Name: $:$Id: TDecompSVD.cxx,v 1.6 2004/02/03 16:50:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Dec 2003
/*************************************************************************
@@ -502,7 +502,7 @@ Bool_t TDecompSVD::Solve(TVectorD &b)
const Double_t threshold = fSig(0)*fTol;
TVectorD tmp(lwb,upb);
- for (Int_t irow = lwb; irow < upb; irow++) {
+ for (Int_t irow = lwb; irow <= upb; irow++) {
Double_t r = 0.0;
if (fSig(irow-lwb) > threshold) {
const TVectorD uc_i = TMatrixDColumn(fU,irow);
@@ -522,30 +522,46 @@ Bool_t TDecompSVD::Solve(TVectorD &b)
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompSVD::Solve(const TVectorD &b,Bool_t &ok)
+{
+// Solve Ax=b assuming the SVD form of A is stored .
+// If A is of size (m x n), the solution will be of size (n) .
+//
+// For m > n , x is the least-squares solution of min(A . x - b)
+
+ TVectorD x = b;
+ ok = Solve(x);
+ if (GetNrows() != GetNcols())
+ x.ResizeTo(GetNcols());
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompSVD::Solve(TMatrixDColumn &cb)
-{
+{
const TMatrixDBase *b = cb.GetMatrix();
- Assert(b->IsValid());
+ Assert(b->IsValid());
Assert(fStatus & kDecomposed);
-
+
if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb())
{
Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
return kFALSE;
}
-
+
// We start with fU fSig fV^T x = b, and turn it into fV^T x = fSig^-1 fU^T b
// Form tmp = fSig^-1 fU^T b but ignore diagonal elements in
// fSig(i) < fTol * max(fSig)
-
+
const Int_t lwb = fU.GetColLwb();
const Int_t upb = lwb+fV.GetNcols()-1;
const Double_t threshold = fSig(0)*fTol;
-
+
TVectorD tmp(lwb,upb);
const TVectorD vb = cb;
- for (Int_t irow = lwb; irow < upb; irow++) {
+ for (Int_t irow = lwb; irow <= upb; irow++) {
Double_t r = 0.0;
if (fSig(irow-lwb) > threshold) {
const TVectorD uc_i = TMatrixDColumn(fU,irow);
@@ -609,6 +625,17 @@ Bool_t TDecompSVD::TransSolve(TVectorD &b)
return kTRUE;
}
+//______________________________________________________________________________
+TVectorD TDecompSVD::TransSolve(const TVectorD &b,Bool_t &ok)
+{
+// Solve A^T x=b assuming the SVD form of A is stored .
+
+ TVectorD x = b;
+ ok = TransSolve(x);
+
+ return x;
+}
+
//______________________________________________________________________________
Bool_t TDecompSVD::TransSolve(TMatrixDColumn &cb)
{
@@ -619,21 +646,21 @@ Bool_t TDecompSVD::TransSolve(TMatrixDColumn &cb)
if (fU.GetNcols() != fU.GetNrows()) {
Error("TransSolve(TMatrixDColumn &","matrix should be square");
return kFALSE;
- }
-
+ }
+
if (fV.GetNrows() != b->GetNrows() || fV.GetRowLwb() != b->GetRowLwb())
- {
+ {
Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
return kFALSE;
}
-
+
// We start with fV fSig fU^T x = b, and turn it into fU^T x = fSig^-1 fV^T b
// Form tmp = fSig^-1 fV^T b but ignore diagonal elements in
// fSig(i) < fTol * max(fSig)
-
+
const Int_t nCol_v = fV.GetNcols();
const Double_t threshold = fSig(0)*fTol;
-
+
const TVectorD vb = cb;
TVectorD tmp(nCol_v);
for (Int_t i = 0; i < nCol_v; i++) {
--
GitLab