diff --git a/smatrix/inc/Math/Dinv.h b/smatrix/inc/Math/Dinv.h
index 6ea9243a27ef444c179f705e796eecc0d304900d..baf42e99c11217f6c3fa82218e77b53d3cbf0ef9 100644
--- a/smatrix/inc/Math/Dinv.h
+++ b/smatrix/inc/Math/Dinv.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: Dinv.h,v 1.3 2005/12/07 16:44:05 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: Dinv.h,v 1.4 2006/02/08 14:45:35 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Dinv
@@ -67,7 +67,7 @@ public:
}
#endif
-
+#ifdef OLD_IMPL
/* Initialized data */
static unsigned int work[n];
@@ -83,6 +83,23 @@ public:
return false;
}
return Dfinv<MatrixRep,n,idim>(rhs,work);
+#else
+
+ /* Initialized data */
+ static unsigned int work[n+1];
+ for(unsigned int i=0; i<n+1; ++i) work[i] = 0;
+
+ static typename MatrixRep::value_type det = 0;
+
+ if (DfactMatrix(rhs,det,work) != 0) {
+ std::cerr << "Dfact_matrix failed!!" << std::endl;
+ return false;
+ }
+
+ int ifail = DfinvMatrix(rhs,work);
+ if (ifail == 0) return true;
+ return false;
+#endif
} // Dinv
@@ -91,6 +108,7 @@ public:
static bool Dinv(MatRepSym<T,idim> & rhs) {
// not very efficient but need to re-do Dsinv for new storage of
// symmetric matrices
+#ifdef OLD_IMPL
MatRepStd<T,idim> tmp;
for (unsigned int i = 0; i< idim*idim; ++i)
tmp[i] = rhs[i];
@@ -100,9 +118,30 @@ public:
rhs[i] = tmp[i];
return true;
+#else
+ int ifail = 0;
+ InvertBunchKaufman(rhs,ifail);
+ if (ifail == 0) return true;
+ return false;
+#endif
}
+ /**
+ Bunch-Kaufman method for inversion of symmetric matrices
+ */
+ template <class T>
+ static int DfactMatrix(MatRepStd<T,idim,n> & rhs, T & det, unsigned int * work);
+ template <class T>
+ static int DfinvMatrix(MatRepStd<T,idim,n> & rhs, unsigned int * work);
+
+ /**
+ Bunch-Kaufman method for inversion of symmetric matrices
+ */
+ template <class T>
+ static void InvertBunchKaufman(MatRepSym<T,idim> & rhs, int &ifail);
+
+
}; // class Inverter
@@ -279,5 +318,6 @@ public:
#include "CramerInversion.icc"
#include "CramerInversionSym.icc"
+#include "MatrixInversion.icc"
#endif /* ROOT_Math_Dinv */
diff --git a/smatrix/inc/Math/MatrixFunctions.h b/smatrix/inc/Math/MatrixFunctions.h
index a190e2d4123f179190b1ce3404a476d8a2831a51..4001fe5e0cb0e89fdaa0f37909277cef6a1e91a7 100644
--- a/smatrix/inc/Math/MatrixFunctions.h
+++ b/smatrix/inc/Math/MatrixFunctions.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: MatrixFunctions.h,v 1.9 2006/03/20 17:11:44 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: MatrixFunctions.h,v 1.10 2006/04/25 13:54:01 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_MatrixFunctions
@@ -57,6 +57,10 @@ SVector<T,D1> operator*(const SMatrix<T,D1,D2,R>& rhs, const SVector<T,D2>& lhs)
#endif
+// matrix-vector product:
+// use apply(i) funciton for matrices. Tested (11/05/06) with using (i,j) but
+// performances are slightly worse (not clear why)
+
//==============================================================================
// meta_row_dot
//==============================================================================
@@ -736,6 +740,88 @@ inline SMatrix<T,D2,D2,MatRepSym<T,D2> > SimilarityT(const Expr<A,T,D1,D2,R>& lh
+//==============================================================================
+// TensorMulOp
+// tensor (or outer) product between two vectors giving a matrix
+//==============================================================================
+template <class Vector1, class Vector2>
+class TensorMulOp {
+public:
+ ///
+ TensorMulOp( const Vector1 & lhs, const Vector2 & rhs) :
+ lhs_(lhs),
+ rhs_(rhs) {}
+
+ ///
+ ~TensorMulOp() {}
+
+ /// Vector2::kSize is the number of columns in the resulting matrix
+ inline typename Vector1::value_type apply(unsigned int i) const {
+ return lhs_.apply( i/ Vector2::kSize) * rhs_.apply( i % Vector2::kSize );
+ }
+ inline typename Vector1::value_type operator() (unsigned int i, unsigned j) const {
+ return lhs_.apply(i) * rhs_.apply(j);
+ }
+
+ inline bool IsInUse (const typename Vector1::value_type * ) const {
+ return false;
+ }
+
+
+protected:
+
+ const Vector1 & lhs_;
+ const Vector2 & rhs_;
+
+};
+
+ // Tensor Prod (use default MatRepStd for the returned expression
+ // cannot make a symmetric matrix
+//==============================================================================
+// TensorProd (SVector x SVector)
+//==============================================================================
+template < class T, unsigned int D1, unsigned int D2>
+inline Expr<TensorMulOp<SVector<T,D1>, SVector<T,D2> >, T, D1, D2 >
+ TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
+ typedef TensorMulOp<SVector<T,D1>, SVector<T,D2> > TVMulOp;
+ return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// TensorProd (VecExpr x SVector)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> >, T, D1, D2 >
+ TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
+ typedef TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> > TVMulOp;
+ return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// TensorProd (SVector x VecExpr)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> >, T, D1, D2 >
+ TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
+ typedef TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> > TVMulOp;
+ return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
+}
+
+
+//==============================================================================
+// TensorProd (VecExpr x VecExpr)
+//==============================================================================
+template <class A, class B, class T, unsigned int D1, unsigned int D2>
+inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> >, T, D1, D2 >
+ TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
+ typedef TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> > TVMulOp;
+ return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
+}
+
+
+
+
+
} // namespace Math
} // namespace ROOT
diff --git a/smatrix/inc/Math/MatrixInversion.icc b/smatrix/inc/Math/MatrixInversion.icc
new file mode 100644
index 0000000000000000000000000000000000000000..8e62d1214a937d5d7a3ffe17c910027990dcef7b
--- /dev/null
+++ b/smatrix/inc/Math/MatrixInversion.icc
@@ -0,0 +1,659 @@
+// @(#)root/smatrix:$Name: $:$Id: MatrixInversion.icc,v 1.1 2006/02/08 14:45:35 moneta Exp $
+// Authors: CLHEP authors, L. Moneta 2006
+
+#include "Math/SVector.h"
+
+
+// inversion algorithms for matrices
+// taken from CLHEP (L. Moneta May 2006)
+
+namespace ROOT {
+
+ namespace Math {
+
+
+ /** Bunch-Kaufman diagonal pivoting method
+ It is decribed in J.R. Bunch, L. Kaufman (1977).
+ "Some Stable Methods for Calculating Inertia and Solving Symmetric
+ Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
+ /Charles F. van Loan, "Matrix Computations" (the second edition
+ has a bug.) and implemented in "lapack"
+ Mario Stanke, 09/97
+
+ */
+
+template <unsigned int idim, unsigned int N>
+template<class T>
+void Inverter<idim,N>::InvertBunchKaufman(MatRepSym<T,idim> & rhs, int &ifail) {
+
+
+
+
+ int i, j, k, s;
+ int pivrow;
+
+ const int nrow = N;
+
+ // Establish the two working-space arrays needed: x and piv are
+ // used as pointers to arrays of doubles and ints respectively, each
+ // of length nrow. We do not want to reallocate each time through
+ // unless the size needs to grow. We do not want to leak memory, even
+ // by having a new without a delete that is only done once.
+
+
+ static SVector<double, N> xvec;
+ static SVector<int, N> pivv;
+
+ typedef int* pivIter;
+ typedef T* mIter;
+
+
+ // Note - resize shuld do nothing if the size is already larger than nrow,
+ // but on VC++ there are indications that it does so we check.
+ // Note - the data elements in a vector are guaranteed to be contiguous,
+ // so x[i] and piv[i] are optimally fast.
+ mIter x = xvec.begin();
+ // x[i] is used as helper storage, needs to have at least size nrow.
+ pivIter piv = pivv.begin();
+ // piv[i] is used to store details of exchanges
+
+ double temp1, temp2;
+ mIter ip, mjj, iq;
+ double lambda, sigma;
+ const double alpha = .6404; // = (1+sqrt(17))/8
+ const double epsilon = 32*std::numeric_limits<T>::epsilon();
+ // whenever a sum of two doubles is below or equal to epsilon
+ // it is set to zero.
+ // this constant could be set to zero but then the algorithm
+ // doesn't neccessarily detect that a matrix is singular
+
+ for (i = 0; i < nrow; i++)
+ piv[i] = i+1;
+
+ ifail = 0;
+
+ // compute the factorization P*A*P^T = L * D * L^T
+ // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
+ // L and D^-1 are stored in A = *this, P is stored in piv[]
+
+ for (j=1; j < nrow; j+=s) // main loop over columns
+ {
+ mjj = rhs.Array() + j*(j-1)/2 + j-1;
+ lambda = 0; // compute lambda = max of A(j+1:n,j)
+ pivrow = j+1;
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (i=j+1; i <= nrow ; ip += i++)
+ if (std::fabs(*ip) > lambda)
+ {
+ lambda = std::fabs(*ip);
+ pivrow = i;
+ }
+
+ if (lambda == 0 )
+ {
+ if (*mjj == 0)
+ {
+ ifail = 1;
+ return;
+ }
+ s=1;
+ *mjj = 1./ *mjj;
+ }
+ else
+ {
+ if (std::fabs(*mjj) >= lambda*alpha)
+ {
+ s=1;
+ pivrow=j;
+ }
+ else
+ {
+ sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
+ ip = rhs.Array() + pivrow*(pivrow-1)/2+j-1;
+ for (k=j; k < pivrow; k++)
+ {
+ if (std::fabs(*ip) > sigma)
+ sigma = std::fabs(*ip);
+ ip++;
+ }
+ if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
+ {
+ s=1;
+ pivrow = j;
+ }
+ else if (std::fabs(*(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1))
+ >= alpha * sigma)
+ s=1;
+ else
+ s=2;
+ }
+ if (pivrow == j) // no permutation neccessary
+ {
+ piv[j-1] = pivrow;
+ if (*mjj == 0)
+ {
+ ifail=1;
+ return;
+ }
+ temp2 = *mjj = 1./ *mjj; // invert D(j,j)
+
+ // update A(j+1:n, j+1,n)
+ for (i=j+1; i <= nrow; i++)
+ {
+ temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
+ ip = rhs.Array()+i*(i-1)/2+j;
+ for (k=j+1; k<=i; k++)
+ {
+ *ip -= temp1 * *(rhs.Array() + k*(k-1)/2 + j-1);
+ if (std::fabs(*ip) <= epsilon)
+ *ip=0;
+ ip++;
+ }
+ }
+ // update L
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (i=j+1; i <= nrow; ip += i++)
+ *ip *= temp2;
+ }
+ else if (s==1) // 1x1 pivot
+ {
+ piv[j-1] = pivrow;
+
+ // interchange rows and columns j and pivrow in
+ // submatrix (j:n,j:n)
+ ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
+ for (i=j+1; i < pivrow; i++, ip++)
+ {
+ temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
+ *(rhs.Array() + i*(i-1)/2 + j-1)= *ip;
+ *ip = temp1;
+ }
+ temp1 = *mjj;
+ *mjj = *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1);
+ *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
+ ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1;
+ iq = ip + pivrow-j;
+ for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
+ {
+ temp1 = *iq;
+ *iq = *ip;
+ *ip = temp1;
+ }
+
+ if (*mjj == 0)
+ {
+ ifail = 1;
+ return;
+ }
+ temp2 = *mjj = 1./ *mjj; // invert D(j,j)
+
+ // update A(j+1:n, j+1:n)
+ for (i = j+1; i <= nrow; i++)
+ {
+ temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
+ ip = rhs.Array()+i*(i-1)/2+j;
+ for (k=j+1; k<=i; k++)
+ {
+ *ip -= temp1 * *(rhs.Array() + k*(k-1)/2 + j-1);
+ if (std::fabs(*ip) <= epsilon)
+ *ip=0;
+ ip++;
+ }
+ }
+ // update L
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (i=j+1; i<=nrow; ip += i++)
+ *ip *= temp2;
+ }
+ else // s=2, ie use a 2x2 pivot
+ {
+ piv[j-1] = -pivrow;
+ piv[j] = 0; // that means this is the second row of a 2x2 pivot
+
+ if (j+1 != pivrow)
+ {
+ // interchange rows and columns j+1 and pivrow in
+ // submatrix (j:n,j:n)
+ ip = rhs.Array() + pivrow*(pivrow-1)/2 + j+1;
+ for (i=j+2; i < pivrow; i++, ip++)
+ {
+ temp1 = *(rhs.Array() + i*(i-1)/2 + j);
+ *(rhs.Array() + i*(i-1)/2 + j) = *ip;
+ *ip = temp1;
+ }
+ temp1 = *(mjj + j + 1);
+ *(mjj + j + 1) =
+ *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
+ *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
+ temp1 = *(mjj + j);
+ *(mjj + j) = *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1);
+ *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1) = temp1;
+ ip = rhs.Array() + (pivrow+1)*pivrow/2 + j;
+ iq = ip + pivrow-(j+1);
+ for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
+ {
+ temp1 = *iq;
+ *iq = *ip;
+ *ip = temp1;
+ }
+ }
+ // invert D(j:j+1,j:j+1)
+ temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
+ if (temp2 == 0)
+ std::cerr
+ << "SymMatrix::bunch_invert: error in pivot choice"
+ << std::endl;
+ temp2 = 1. / temp2;
+ // this quotient is guaranteed to exist by the choice
+ // of the pivot
+ temp1 = *mjj;
+ *mjj = *(mjj + j + 1) * temp2;
+ *(mjj + j + 1) = temp1 * temp2;
+ *(mjj + j) = - *(mjj + j) * temp2;
+
+ if (j < nrow-1) // otherwise do nothing
+ {
+ // update A(j+2:n, j+2:n)
+ for (i=j+2; i <= nrow ; i++)
+ {
+ ip = rhs.Array() + i*(i-1)/2 + j-1;
+ temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
+ if (std::fabs(temp1 ) <= epsilon)
+ temp1 = 0;
+ temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
+ if (std::fabs(temp2 ) <= epsilon)
+ temp2 = 0;
+ for (k = j+2; k <= i ; k++)
+ {
+ ip = rhs.Array() + i*(i-1)/2 + k-1;
+ iq = rhs.Array() + k*(k-1)/2 + j-1;
+ *ip -= temp1 * *iq + temp2 * *(iq+1);
+ if (std::fabs(*ip) <= epsilon)
+ *ip = 0;
+ }
+ }
+ // update L
+ for (i=j+2; i <= nrow ; i++)
+ {
+ ip = rhs.Array() + i*(i-1)/2 + j-1;
+ temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
+ if (std::fabs(temp1) <= epsilon)
+ temp1 = 0;
+ *(ip+1) = *ip * *(mjj + j)
+ + *(ip+1) * *(mjj + j + 1);
+ if (std::fabs(*(ip+1)) <= epsilon)
+ *(ip+1) = 0;
+ *ip = temp1;
+ }
+ }
+ }
+ }
+ } // end of main loop over columns
+
+ if (j == nrow) // the the last pivot is 1x1
+ {
+ mjj = rhs.Array() + j*(j-1)/2 + j-1;
+ if (*mjj == 0)
+ {
+ ifail = 1;
+ return;
+ }
+ else
+ *mjj = 1. / *mjj;
+ } // end of last pivot code
+
+ // computing the inverse from the factorization
+
+ for (j = nrow ; j >= 1 ; j -= s) // loop over columns
+ {
+ mjj = rhs.Array() + j*(j-1)/2 + j-1;
+ if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
+ {
+ s = 1;
+ if (j < nrow)
+ {
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (i=0; i < nrow-j; ip += 1+j+i++)
+ x[i] = *ip;
+ for (i=j+1; i<=nrow ; i++)
+ {
+ temp2=0;
+ ip = rhs.Array() + i*(i-1)/2 + j;
+ for (k=0; k <= i-j-1; k++)
+ temp2 += *ip++ * x[k];
+ for (ip += i-1; k < nrow-j; ip += 1+j+k++)
+ temp2 += *ip * x[k];
+ *(rhs.Array()+ i*(i-1)/2 + j-1) = -temp2;
+ }
+ temp2 = 0;
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (k=0; k < nrow-j; ip += 1+j+k++)
+ temp2 += x[k] * *ip;
+ *mjj -= temp2;
+ }
+ }
+ else //2x2 pivot, compute columns j and j-1 of the inverse
+ {
+ if (piv[j-1] != 0)
+ std::cerr << "error in piv" << piv[j-1] << std::endl;
+ s=2;
+ if (j < nrow)
+ {
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (i=0; i < nrow-j; ip += 1+j+i++)
+ x[i] = *ip;
+ for (i=j+1; i<=nrow ; i++)
+ {
+ temp2 = 0;
+ ip = rhs.Array() + i*(i-1)/2 + j;
+ for (k=0; k <= i-j-1; k++)
+ temp2 += *ip++ * x[k];
+ for (ip += i-1; k < nrow-j; ip += 1+j+k++)
+ temp2 += *ip * x[k];
+ *(rhs.Array()+ i*(i-1)/2 + j-1) = -temp2;
+ }
+ temp2 = 0;
+ ip = rhs.Array() + (j+1)*j/2 + j-1;
+ for (k=0; k < nrow-j; ip += 1+j+k++)
+ temp2 += x[k] * *ip;
+ *mjj -= temp2;
+ temp2 = 0;
+ ip = rhs.Array() + (j+1)*j/2 + j-2;
+ for (i=j+1; i <= nrow; ip += i++)
+ temp2 += *ip * *(ip+1);
+ *(mjj-1) -= temp2;
+ ip = rhs.Array() + (j+1)*j/2 + j-2;
+ for (i=0; i < nrow-j; ip += 1+j+i++)
+ x[i] = *ip;
+ for (i=j+1; i <= nrow ; i++)
+ {
+ temp2 = 0;
+ ip = rhs.Array() + i*(i-1)/2 + j;
+ for (k=0; k <= i-j-1; k++)
+ temp2 += *ip++ * x[k];
+ for (ip += i-1; k < nrow-j; ip += 1+j+k++)
+ temp2 += *ip * x[k];
+ *(rhs.Array()+ i*(i-1)/2 + j-2)= -temp2;
+ }
+ temp2 = 0;
+ ip = rhs.Array() + (j+1)*j/2 + j-2;
+ for (k=0; k < nrow-j; ip += 1+j+k++)
+ temp2 += x[k] * *ip;
+ *(mjj-j) -= temp2;
+ }
+ }
+
+ // interchange rows and columns j and piv[j-1]
+ // or rows and columns j and -piv[j-2]
+
+ pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
+ ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
+ for (i=j+1;i < pivrow; i++, ip++)
+ {
+ temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
+ *(rhs.Array() + i*(i-1)/2 + j-1) = *ip;
+ *ip = temp1;
+ }
+ temp1 = *mjj;
+ *mjj = *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
+ *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
+ if (s==2)
+ {
+ temp1 = *(mjj-1);
+ *(mjj-1) = *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2);
+ *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2) = temp1;
+ }
+
+ ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
+ iq = ip + pivrow-j;
+ for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
+ {
+ temp1 = *iq;
+ *iq = *ip;
+ *ip = temp1;
+ }
+ } // end of loop over columns (in computing inverse from factorization)
+
+ return; // inversion successful
+
+}
+
+
+
+
+ /**
+ LU factorization : code originally from CERNLIB dfact routine and ported in C++ for CLHEP
+ */
+
+template <unsigned int idim, unsigned int n>
+template<class T>
+int Inverter<idim,n>::DfactMatrix(MatRepStd<T,idim,n> & rhs, T &det, unsigned int *ir) {
+
+ if (idim != n) return -1;
+
+ int ifail, jfail;
+
+ typedef T* mIter;
+
+
+ double tf;
+ double g1 = 1.0e-19, g2 = 1.0e19;
+
+ double p, q, t;
+ double s11, s12;
+
+ double epsilon = 8*std::numeric_limits<double>::epsilon();
+ // could be set to zero (like it was before)
+ // but then the algorithm often doesn't detect
+ // that a matrix is singular
+
+ int normal = 0, imposs = -1;
+ int jrange = 0, jover = 1, junder = -1;
+ ifail = normal;
+ jfail = jrange;
+ int nxch = 0;
+ det = 1.0;
+ mIter mj = rhs.Array();
+ mIter mjj = mj;
+ for (unsigned int j=1;j<=n;j++) {
+ unsigned int k = j;
+ p = (std::fabs(*mjj));
+ if (j!=n) {
+ mIter mij = mj + n + j - 1;
+ for (unsigned int i=j+1;i<=n;i++) {
+ q = (std::fabs(*(mij)));
+ if (q > p) {
+ k = i;
+ p = q;
+ }
+ mij += n;
+ }
+ if (k==j) {
+ if (p <= epsilon) {
+ det = 0;
+ ifail = imposs;
+ jfail = jrange;
+ return ifail;
+ }
+ det = -det; // in this case the sign of the determinant
+ // must not change. So I change it twice.
+ }
+ mIter mjl = mj;
+ mIter mkl = rhs.Array() + (k-1)*n;
+ for (unsigned int l=1;l<=n;l++) {
+ tf = *mjl;
+ *(mjl++) = *mkl;
+ *(mkl++) = tf;
+ }
+ nxch = nxch + 1; // this makes the determinant change its sign
+ ir[nxch] = (((j)<<12)+(k));
+ } else {
+ if (p <= epsilon) {
+ det = 0.0;
+ ifail = imposs;
+ jfail = jrange;
+ return ifail;
+ }
+ }
+ det *= *mjj;
+ *mjj = 1.0 / *mjj;
+ t = (std::fabs(det));
+ if (t < g1) {
+ det = 0.0;
+ if (jfail == jrange) jfail = junder;
+ } else if (t > g2) {
+ det = 1.0;
+ if (jfail==jrange) jfail = jover;
+ }
+ if (j!=n) {
+ mIter mk = mj + n;
+ mIter mkjp = mk + j;
+ mIter mjk = mj + j;
+ for (unsigned k=j+1;k<=n;k++) {
+ s11 = - (*mjk);
+ s12 = - (*mkjp);
+ if (j!=1) {
+ mIter mik = rhs.Array() + k - 1;
+ mIter mijp = rhs.Array() + j;
+ mIter mki = mk;
+ mIter mji = mj;
+ for (unsigned int i=1;i<j;i++) {
+ s11 += (*mik) * (*(mji++));
+ s12 += (*mijp) * (*(mki++));
+ mik += n;
+ mijp += n;
+ }
+ }
+ *(mjk++) = -s11 * (*mjj);
+ *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
+ mk += n;
+ mkjp += n;
+ }
+ }
+ mj += n;
+ mjj += (n+1);
+ }
+ if (nxch%2==1) det = -det;
+ if (jfail !=jrange) det = 0.0;
+ ir[n] = nxch;
+ return 0;
+}
+
+
+
+ /**
+ Inversion for General square matrices. Code from dfinv routine from CERNLIB
+ Assumed first the LU decomposition via DfactMatrix function
+
+ taken from CLHEP : L. Moneta May 2006
+ */
+
+template <unsigned int idim, unsigned int n>
+template<class T>
+int Inverter<idim,n>::DfinvMatrix(MatRepStd<T,idim,n> & rhs,unsigned int * ir) {
+
+
+ typedef T* mIter;
+
+ if (idim != n) return -1;
+
+
+ double s31, s32;
+ register double s33, s34;
+
+ mIter m11 = rhs.Array();
+ mIter m12 = m11 + 1;
+ mIter m21 = m11 + n;
+ mIter m22 = m12 + n;
+ *m21 = -(*m22) * (*m11) * (*m21);
+ *m12 = -(*m12);
+ if (n>2) {
+ mIter mi = rhs.Array() + 2 * n;
+ mIter mii= rhs.Array() + 2 * n + 2;
+ mIter mimim = rhs.Array() + n + 1;
+ for (unsigned int i=3;i<=n;i++) {
+ unsigned int im2 = i - 2;
+ mIter mj = rhs.Array();
+ mIter mji = mj + i - 1;
+ mIter mij = mi;
+ for (unsigned int j=1;j<=im2;j++) {
+ s31 = 0.0;
+ s32 = *mji;
+ mIter mkj = mj + j - 1;
+ mIter mik = mi + j - 1;
+ mIter mjkp = mj + j;
+ mIter mkpi = mj + n + i - 1;
+ for (unsigned int k=j;k<=im2;k++) {
+ s31 += (*mkj) * (*(mik++));
+ s32 += (*(mjkp++)) * (*mkpi);
+ mkj += n;
+ mkpi += n;
+ }
+ *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
+ *mji = -s32;
+ mj += n;
+ mji += n;
+ mij++;
+ }
+ *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
+ *(mimim+1) = -(*(mimim+1));
+ mi += n;
+ mimim += (n+1);
+ mii += (n+1);
+ }
+ }
+ mIter mi = rhs.Array();
+ mIter mii = rhs.Array();
+ for (unsigned int i=1;i<n;i++) {
+ unsigned int ni = n - i;
+ mIter mij = mi;
+ //int j;
+ for (unsigned j=1; j<=i;j++) {
+ s33 = *mij;
+ register mIter mikj = mi + n + j - 1;
+ register mIter miik = mii + 1;
+ mIter min_end = mi + n;
+ for (;miik<min_end;) {
+ s33 += (*mikj) * (*(miik++));
+ mikj += n;
+ }
+ *(mij++) = s33;
+ }
+ for (unsigned j=1;j<=ni;j++) {
+ s34 = 0.0;
+ mIter miik = mii + j;
+ mIter mikij = mii + j * n + j;
+ for (unsigned int k=j;k<=ni;k++) {
+ s34 += *mikij * (*(miik++));
+ mikij += n;
+ }
+ *(mii+j) = s34;
+ }
+ mi += n;
+ mii += (n+1);
+ }
+ unsigned int nxch = ir[n];
+ if (nxch==0) return 0;
+ for (unsigned int mm=1;mm<=nxch;mm++) {
+ unsigned int k = nxch - mm + 1;
+ int ij = ir[k];
+ int i = ij >> 12;
+ int j = ij%4096;
+ mIter mki = rhs.Array() + i - 1;
+ mIter mkj = rhs.Array() + j - 1;
+ for (k=1; k<=n;k++) {
+ // 2/24/05 David Sachs fix of improper swap bug that was present
+ // for many years:
+ double ti = *mki; // 2/24/05
+ *mki = *mkj;
+ *mkj = ti; // 2/24/05
+ mki += n;
+ mkj += n;
+ }
+ }
+ return 0;
+}
+
+
+ } // end namespace Math
+} // end namespace ROOT
diff --git a/smatrix/inc/Math/SMatrix.h b/smatrix/inc/Math/SMatrix.h
index f28e456da0fc8f5482c3597af1480fe9ba21b807..bfc5eaf8eb1ba5f62c3002dd3de68a0f96b25c14 100644
--- a/smatrix/inc/Math/SMatrix.h
+++ b/smatrix/inc/Math/SMatrix.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SMatrix.h,v 1.19 2006/04/20 13:13:21 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.h,v 1.20 2006/04/25 13:54:01 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SMatrix
@@ -337,6 +337,7 @@ public:
#endif
+#ifdef OLD_IMPL
/** @name --- Expert functions --- */
/**
@@ -362,19 +363,23 @@ public:
this method will preserve the contents of the Matrix!
*/
bool Sdet2(T& det) const;
+#endif
/**
- invert square Matrix via Dinv.
- This method change the current matrix
+ invert square Matrix ( this method change the current matrix)
+ The method used for general square matrices is the LU factorization taken from Dinv routine
+ from the CERNLIB (written in C++ from CLHEP authors)
+ In case of symmetric matrices Bunch-Kaufman diagonal pivoting method is used
+ (The implementation is the one written by the CLHEP authors)
*/
bool Invert();
/**
- invert square Matrix via Dinv.
- This method returns a new matrix. In case the inversion fails
- the current matrix is returned
- Return ifail = 0 when successfull
+ invert a square Matrix and returns a new matrix. In case the inversion fails
+ the current matrix is returned.
+ Return ifail = 0 when successfull.
+ See ROOT::Math::SMatrix::Invert for the inversion algorithm
*/
SMatrix<T,D1,D2,R> Inverse(int & ifail ) const;
diff --git a/smatrix/inc/Math/SMatrix.icc b/smatrix/inc/Math/SMatrix.icc
index 512986519417b8eb647c3001b72b1aa889c83124..7efd8041916a35a3f3a861322388c5822fd06cb2 100644
--- a/smatrix/inc/Math/SMatrix.icc
+++ b/smatrix/inc/Math/SMatrix.icc
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SMatrix.icc,v 1.19 2006/03/30 16:18:05 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.icc,v 1.20 2006/04/25 13:54:01 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SMatrix_icc
@@ -431,6 +431,7 @@ bool SMatrix<T,D1,D2,R>::operator<(const Expr<A,T,D1,D2,R2>& rhs) const {
return rc;
}
+#ifdef OLD_IMPL
//==============================================================================
// sinvert
//==============================================================================
@@ -451,7 +452,6 @@ inline SMatrix<T,D1,D2,R> SMatrix<T,D1,D2,R>::Sinverse(int & ifail) const {
}
-
//==============================================================================
// sdet
//==============================================================================
@@ -469,6 +469,7 @@ inline bool SMatrix<T,D1,D2,R>::Sdet2(T& det) const {
return tmp.Sdet(det);
}
+#endif
//==============================================================================
// invert
@@ -476,7 +477,7 @@ inline bool SMatrix<T,D1,D2,R>::Sdet2(T& det) const {
template <class T, unsigned int D1, unsigned int D2, class R>
inline bool SMatrix<T,D1,D2,R>::Invert() {
STATIC_CHECK( D1 == D2,SMatrix_not_square);
- return Inverter<D1>::Dinv((*this).fRep);
+ return Inverter<D1,D2>::Dinv((*this).fRep);
}
template <class T, unsigned int D1, unsigned int D2, class R>
diff --git a/smatrix/inc/Math/SVector.h b/smatrix/inc/Math/SVector.h
index 82174e70a0627c6418516656f9853635d85a8e6d..22239f94654a82d2e3f1176f2098f10287ae78f6 100644
--- a/smatrix/inc/Math/SVector.h
+++ b/smatrix/inc/Math/SVector.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SVector.h,v 1.7 2006/02/27 18:41:58 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SVector.h,v 1.8 2006/04/20 13:13:21 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SVector
@@ -119,7 +119,12 @@ public:
// if you use iterator this is not necessary
/// fill from array, len must be equal to D!
- SVector(const T* a, unsigned int len);
+ SVector( const T * a, unsigned int len);
+
+ /// fill from iterators
+ //(iterator is T* to skip ambiguities)
+ SVector(const_iterator begin, const_iterator end);
+
#endif
///
SVector(const T& rhs);
diff --git a/smatrix/inc/Math/SVector.icc b/smatrix/inc/Math/SVector.icc
index 427e01a50f5299738db25c9aa2c19196e3ca130c..10f36322733fce3e2fa019b731594cd75dec263a 100644
--- a/smatrix/inc/Math/SVector.icc
+++ b/smatrix/inc/Math/SVector.icc
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SVector.icc,v 1.6 2006/02/08 14:45:35 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SVector.icc,v 1.7 2006/02/27 18:41:58 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SVector_icc
@@ -114,6 +114,13 @@ SVector<T,D>::SVector(const T* a, unsigned int len) {
fArray[i] = a[i];
}
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const_iterator begin, const_iterator end) {
+ assert( begin + D == end);
+ std::copy(begin, end, fArray);
+}
+
+
#endif
template <class T, unsigned int D>
diff --git a/smatrix/test/Makefile b/smatrix/test/Makefile
index a16ed9d748d435fa78a412de40552999f1c737b2..2ecde26e1fd1f8864f1d299dd4043e5fff235c93 100644
--- a/smatrix/test/Makefile
+++ b/smatrix/test/Makefile
@@ -12,7 +12,7 @@ include ../../test/Makefile.arch
ifeq ($(PLATFORM),macosx)
#unroll loop better on gcc > 4
-CXXFLAGS+= -g -funroll-loops
+CXXFLAGS+= -O3 -g -funroll-loops
endif
# if have clhep
diff --git a/smatrix/test/stressKalman.cxx b/smatrix/test/stressKalman.cxx
index d53a4a69de9767f5f47c4147c9bdfb22351afd09..db93ffbf7fea98d3ec92744304b0928426cd7b87 100644
--- a/smatrix/test/stressKalman.cxx
+++ b/smatrix/test/stressKalman.cxx
@@ -19,7 +19,7 @@
#include <map>
-
+//#undef HAVE_CLHEP
#ifdef HAVE_CLHEP
#include "CLHEP/Matrix/SymMatrix.h"
@@ -58,7 +58,7 @@ public:
test_smatrix_kalman();
test_tmatrix_kalman();
#ifdef HAVE_CLHEP
- test_clhep_matrix();
+ test_clhep_kalman();
#endif
return;
@@ -188,8 +188,8 @@ void printTime(TStopwatch & time, std::string s) {
}
// reference times for sizes <=6 and > 6 on Linux slc3 P4 3Ghz ("SMatrix","SMatrix_sym","TMatrix")
-double refTime1[3] = { 40.49, 53.75, 83.21 };
-double refTime2[3] = { 393.81, 462.16, 785.50 };
+double refTime1[4] = { 40.49, 53.75, 83.21,1000 };
+double refTime2[4] = { 393.81, 462.16, 785.50,10000 };
#define NMAX1 9 // matrix storese results from 2 to 10
#define NMAX2 7 // results from 2 to 8
@@ -807,8 +807,8 @@ int TestRunner<NDIM1,NDIM2>::test_clhep_kalman() {
fillRandomMat(r,K0,second,first,1);
fillRandomSym(r,Cp,second,1);
fillRandomSym(r,V,first,1);
- fillRandomVec(r,m,first,1);
- fillRandomVec(r,xp,second,1);
+ fillRandomVec(r,m,first);
+ fillRandomVec(r,xp,second);
MnSymMatrix I(second,1);//Identity matrix
@@ -843,12 +843,11 @@ int TestRunner<NDIM1,NDIM2>::test_clhep_kalman() {
x2= RinvSym.similarity(vtmp1);
if(ifail!=0) { std::cout << "Error inverting Rinv" << std::endl; break; }
}
- // std::cout << k << " chi2 " << x2 << std::endl;
x2sum += x2;
c2 = 0;
- for (int i=1; i<=NDIM2; ++i)
- for (int j=1; j<=NDIM2; ++j)
+ for (unsigned int i=1; i<=NDIM2; ++i)
+ for (unsigned int j=1; j<=NDIM2; ++j)
c2 += C(i,j);
c2sum += c2;
}
diff --git a/smatrix/test/testKalman.cxx b/smatrix/test/testKalman.cxx
index bbc408697915257fdffa219f444375078016d076..c9ee5ce65f92523d3237ceb617139034a050a86a 100644
--- a/smatrix/test/testKalman.cxx
+++ b/smatrix/test/testKalman.cxx
@@ -489,8 +489,8 @@ int test_clhep_kalman() {
fillRandomMat(r,K0,second,first,1);
fillRandomSym(r,Cp,second,1);
fillRandomSym(r,V,first,1);
- fillRandomVec(r,m,first,1);
- fillRandomVec(r,xp,second,1);
+ fillRandomVec(r,m,first);
+ fillRandomVec(r,xp,second);
MnSymMatrix I(second,1);//Identity matrix
diff --git a/smatrix/test/testOperations.cxx b/smatrix/test/testOperations.cxx
index be886229a3f340d69f3bb1f78f2eb63a85648370..d1deeca31b1653f9a5896d54de65dec1db6600cd 100644
--- a/smatrix/test/testOperations.cxx
+++ b/smatrix/test/testOperations.cxx
@@ -796,8 +796,9 @@ int main(int argc , char *argv[] ) {
NLOOP = 1000*NLOOP_MIN
initValues();
-
TEST(5)
+// NLOOP = 50*NLOOP_MIN;
+// TEST(30);
#else
NLOOP = 5000*NLOOP_MIN;
initValues();
diff --git a/smatrix/test/testSMatrix.cxx b/smatrix/test/testSMatrix.cxx
index 70b0b7fe0e1d1525ecb1d7aad12b9dd6583da1b8..50a45c741dfb4742bd7767026f1ab8a9bd2b5f8e 100644
--- a/smatrix/test/testSMatrix.cxx
+++ b/smatrix/test/testSMatrix.cxx
@@ -79,11 +79,11 @@ int test1() {
//std::vector<float> vp(sp.begin(), sp.end() );
std::vector<float> vm(sm.begin(), sm.end() );
- //SVector<float, 2> sp2(vp.begin(),vp.end());
- //SVector<float, 2> sp2(vp.begin(),vp.size());
+ SVector<float, 4> sp1(m,4);
+ SVector<float, 4> sp2(sm.begin(),sm.end());
SMatrix<float, 2,2> sm2(vm.begin(),vm.end());
- //if ( sp2 != sp) { cout << "Test STL interface for SVector failed" << endl; return -1; }
+ if ( sp1 != sp2) { cout << "Test STL interface for SVector failed" << endl; return -1; }
if ( sm2 != sm) { cout << "Test STL interface for SMatrix failed" << endl; return -1; }
@@ -137,7 +137,7 @@ int test3() {
cout << "A: " << endl << A << endl;
double det = 0.;
- A.Sdet(det);
+ A.Det(det);
cout << "Determinant: " << det << endl;
// WARNING: A has changed!!
cout << "A again: " << endl << A << endl;
@@ -446,6 +446,24 @@ int test12() {
double a = 3;
std::cout << "S\n" << a * S << std::endl;
+ // test inversion
+ int ifail1 = 0;
+ //int ifail2 = 0;
+ SMatrix<double,2,2,MatRepSym<double,2> > Sinv1 = S.Inverse (ifail1);
+ //SMatrix<double,2,2,MatRepSym<double,2> > Sinv2 = S.Sinverse(ifail2);
+ std::cout << "Inverse: S-1 " << Sinv1 << "\nifail=" << ifail1 << std::endl;
+ //std::cout << "Sinverse: S-1" << Sinv2 << "\nifail=" << ifail2 << std::endl;
+
+ SMatrix<double,2> IS1 = S*Sinv1;
+ //SMatrix<double,2> IS2 = S*Sinv2;
+ double d1 = std::sqrt( IS1(1,0)*IS1(1,0) + IS1(0,1)*IS1(0,1) );
+ //double d2 = std::sqrt( IS2(1,0)*IS2(1,0) + IS2(0,1)*IS2(0,1) );
+
+
+ iret |= compare( d1 < 1E-6,true,"inversion1" );
+ //iret |= compare( d2 < 1E-6,true,"inversion2" );
+
+
SMatrix<double,2,3 > M1 = SMatrixIdentity();
M1(0,1) = 1; M1(1,0) = 1; M1(0,2) = 1; M1(1,2) = 1;
SMatrix<double,2,3 > M2 = SMatrixIdentity();
@@ -881,6 +899,42 @@ int test16() {
}
+int test17() {
+ int iret =0;
+ // tets tensor product
+ SVector<double,2> v1(1,2);
+ SVector<double,3> v2(1,2,3);
+
+ SMatrix<double,2,3> m = TensorProd(v1,v2);
+ std::cout << "Tensor Product \n" << m << std::endl;
+
+ SVector<double,4> a1(2,-1,3,4);
+ SVector<double,4> a2(5,6,1,-2);
+
+ SMatrix<double,4> A = TensorProd(a1,a2);
+ double r1 = Dot(a1, A * a2 );
+ double r2 = Dot(a1, a1) * Dot(a2,a2 );
+ iret |= compare(r1,r2,"tensor prod");
+
+ A = TensorProd(2.*a1,a2);
+ r1 = Dot(a1, A * a2 )/2;
+ r2 = Dot(a1, a1) * Dot(a2,a2 );
+ iret |= compare(r1,r2,"tensor prod");
+
+
+ A = TensorProd(a1,2*a2);
+ r1 = Dot(a1, A * a2 )/2;
+ r2 = Dot(a1, a1) * Dot(a2,a2 );
+ iret |= compare(r1,r2,"tensor prod");
+
+ A = TensorProd(0.5*a1,2*a2);
+ r1 = Dot(a1, A * a2 );
+ r2 = Dot(a1, a1) * Dot(a2,a2 );
+ iret |= compare(r1,r2,"tensor prod");
+
+ return iret;
+}
+
#define TEST(N) \
itest = N; \
@@ -909,6 +963,7 @@ int main() {
TEST(14);
TEST(15);
TEST(16);
+ TEST(17);