diff --git a/smatrix/CreateBinaryOp.sh b/smatrix/CreateBinaryOp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..2c2edd753411cb8bffb52e9bea42e950069019b9
--- /dev/null
+++ b/smatrix/CreateBinaryOp.sh
@@ -0,0 +1,208 @@
+#!/bin/bash
+#
+# Script to generate binary operators for SVector & SMatrix class.
+# The operator is applied for each element:
+# C(i,j) = A(i,j) OP B(i,j)
+#
+# author: T. Glebe
+#
+# MPI fuer Kernphysik
+# Saupfercheckweg 1
+# D-69117 Heidelberg
+#
+# 19. Mar 2001 (TG) created
+#
+##########################################################################
+
+OUTPUTFILE="BinaryOperators.hh"
+
+# Format: Class name, Function name, operator
+OPLIST="
+AddOp,operator+,+
+MinOp,operator-,-
+MulOp,operator*,*
+DivOp,operator/,/
+"
+
+
+# generate code:
+(
+echo "\
+#ifndef __BINARYOPERATORS_HH
+#define __BINARYOPERATORS_HH
+//======================================================
+//
+// ATTENTION: This file was automatically generated,
+// do not edit!
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+//======================================================
+
+template <class T, unsigned int D> class SVector;
+template <class T, unsigned int D1, unsigned int D2> class SMatrix;
+"
+
+ for i in $OPLIST; do
+ CNAM=`echo $i | cut -d, -f1`
+ FNAM=`echo $i | cut -d, -f2`
+ OPER=`echo $i | cut -d, -f3`
+
+ echo "
+//==============================================================================
+// $CNAM
+//==============================================================================
+template <class T>
+class $CNAM {
+public:
+ static inline T apply(const T& lhs, const T& rhs) {
+ return lhs ${OPER} rhs;
+ }
+};
+"
+
+# =========== SVector / Expr ===============
+CMBLIST="
+ +SVector<T,D>++SVector<T,D>
+ class_A,+Expr<A,T,D>++SVector<T,D>
+ +SVector<T,D>+class_A,+Expr<A,T,D>
+ class_A,+Expr<A,T,D>+class_B,+Expr<B,T,D>
+"
+ for j in $CMBLIST; do
+ CL1=`echo $j | cut -d+ -f1 | tr '_' ' '`
+ TY1=`echo $j | cut -d+ -f2 | tr '_' ' '`
+ CL2=`echo $j | cut -d+ -f3 | tr '_' ' '`
+ TY2=`echo $j | cut -d+ -f4 | tr '_' ' '`
+
+echo "
+//==============================================================================
+// $FNAM (SVector, binary)
+//==============================================================================
+template <${CL1} ${CL2} class T, unsigned int D>
+inline Expr<BinaryOp<${CNAM}<T>, ${TY1}, ${TY2}, T>, T, D>
+ ${FNAM}(const ${TY1}& lhs, const ${TY2}& rhs) {
+ typedef BinaryOp<${CNAM}<T>, ${TY1}, ${TY2}, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D>(${CNAM}BinOp(${CNAM}<T>(),lhs,rhs));
+}
+"
+ done
+
+# =========== SVector Constant ===============
+CNSTLIST="
+ +SVector<T,D>
+ class_B,+Expr<B,T,D>
+"
+
+ for j in $CNSTLIST; do
+ CL1=`echo $j | cut -d+ -f1 | tr '_' ' '`
+ TY1=`echo $j | cut -d+ -f2 | tr '_' ' '`
+
+echo "
+//==============================================================================
+// $FNAM (SVector, binary, Constant)
+//==============================================================================
+template <class A, ${CL1} class T, unsigned int D>
+inline Expr<BinaryOp<${CNAM}<T>, ${TY1}, Constant<A>, T>, T, D>
+ ${FNAM}(const ${TY1}& lhs, const A& rhs) {
+ typedef BinaryOp<${CNAM}<T>, ${TY1}, Constant<A>, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D>(${CNAM}BinOp(${CNAM}<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// $FNAM (SVector, binary, Constant)
+//==============================================================================
+template <class A, ${CL1} class T, unsigned int D>
+inline Expr<BinaryOp<${CNAM}<T>, Constant<A>, ${TY1}, T>, T, D>
+ ${FNAM}(const A& lhs, const ${TY1}& rhs) {
+ typedef BinaryOp<${CNAM}<T>, Constant<A>, ${TY1}, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D>(${CNAM}BinOp(${CNAM}<T>(),Constant<A>(lhs),rhs));
+}
+"
+ done
+
+
+# =========== SMatrix / Expr ===============
+CMBLIST="
+ +SMatrix<T,D,D2>++SMatrix<T,D,D2>
+ class_A,+Expr<A,T,D,D2>++SMatrix<T,D,D2>
+ +SMatrix<T,D,D2>+class_A,+Expr<A,T,D,D2>
+ class_A,+Expr<A,T,D,D2>+class_B,+Expr<B,T,D,D2>
+"
+
+# component wise multiplication should not occupy operator*()
+ if [ "$OPER" == "*" ]; then
+ MFNAM="times"
+ else
+ MFNAM=$FNAM
+ fi
+
+ for j in $CMBLIST; do
+ CL1=`echo $j | cut -d+ -f1 | tr '_' ' '`
+ TY1=`echo $j | cut -d+ -f2 | tr '_' ' '`
+ CL2=`echo $j | cut -d+ -f3 | tr '_' ' '`
+ TY2=`echo $j | cut -d+ -f4 | tr '_' ' '`
+
+echo "
+//==============================================================================
+// $MFNAM (SMatrix, binary)
+//==============================================================================
+template <${CL1} ${CL2} class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<${CNAM}<T>, ${TY1}, ${TY2}, T>, T, D, D2>
+ ${MFNAM}(const ${TY1}& lhs, const ${TY2}& rhs) {
+ typedef BinaryOp<${CNAM}<T>, ${TY1}, ${TY2}, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D,D2>(${CNAM}BinOp(${CNAM}<T>(),lhs,rhs));
+}
+"
+ done
+
+
+# =========== SMatrix Constant ===============
+CNSTLIST="
+ +SMatrix<T,D,D2>
+ class_B,+Expr<B,T,D,D2>
+"
+
+ for j in $CNSTLIST; do
+ CL1=`echo $j | cut -d+ -f1 | tr '_' ' '`
+ TY1=`echo $j | cut -d+ -f2 | tr '_' ' '`
+
+echo "
+//==============================================================================
+// $FNAM (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, ${CL1} class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<${CNAM}<T>, ${TY1}, Constant<A>, T>, T, D, D2>
+ ${FNAM}(const ${TY1}& lhs, const A& rhs) {
+ typedef BinaryOp<${CNAM}<T>, ${TY1}, Constant<A>, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D,D2>(${CNAM}BinOp(${CNAM}<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// $FNAM (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, ${CL1} class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<${CNAM}<T>, Constant<A>, ${TY1}, T>, T, D, D2>
+ ${FNAM}(const A& lhs, const ${TY1}& rhs) {
+ typedef BinaryOp<${CNAM}<T>, Constant<A>, ${TY1}, T> ${CNAM}BinOp;
+
+ return Expr<${CNAM}BinOp,T,D,D2>(${CNAM}BinOp(${CNAM}<T>(),Constant<A>(lhs),rhs));
+}
+"
+ done
+done
+
+echo "#endif"
+) > $OUTPUTFILE
+
+echo "$OUTPUTFILE generated"
\ No newline at end of file
diff --git a/smatrix/CreateUnaryOp.sh b/smatrix/CreateUnaryOp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..4c2866f1093cdd17bbf73f3e8bb69af24e28c04e
--- /dev/null
+++ b/smatrix/CreateUnaryOp.sh
@@ -0,0 +1,120 @@
+#!/bin/bash
+#
+# Script to generate unary operators for SVector and SMatrix class.
+# The operator is applied for each element:
+# C(i,j) = OP( B(i,j) )
+#
+# author: T. Glebe
+#
+# MPI fuer Kernphysik
+# Saupfercheckweg 1
+# D-69117 Heidelberg
+#
+# 19. Mar 2001 (TG) created
+# 04. Apr 2001 (TG) Sqrt added
+#
+##########################################################################
+
+OUTPUTFILE="UnaryOperators.hh"
+
+# Format: Class name, Function name, operator
+OPLIST="
+Minus,operator-,-
+Fabs,fabs,fabs
+Sqr,sqr,square
+Sqrt,sqrt,sqrt
+"
+
+# generate code:
+(
+echo "\
+#ifndef __UNARYOPERATORS_HH
+#define __UNARYOPERATORS_HH
+//======================================================
+//
+// ATTENTION: This file was automatically generated,
+// do not edit!
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+//======================================================
+
+template <class T, unsigned int D> class SVector;
+template <class T, unsigned int D1, unsigned int D2> class SMatrix;
+"
+
+ for i in $OPLIST; do
+ CNAM=`echo $i | cut -d, -f1`
+ FNAM=`echo $i | cut -d, -f2`
+ OPER=`echo $i | cut -d, -f3`
+
+ echo "
+//==============================================================================
+// $CNAM
+//==============================================================================
+template <class T>
+class $CNAM {
+public:
+ static inline T apply(const T& rhs) {
+ return ${OPER}(rhs);
+ }
+};
+
+//==============================================================================
+// $FNAM (Expr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<UnaryOp<${CNAM}<T>, Expr<A,T,D>, T>, T, D>
+ ${FNAM}(const Expr<A,T,D>& rhs) {
+ typedef UnaryOp<${CNAM}<T>, Expr<A,T,D>, T> ${CNAM}UnaryOp;
+
+ return Expr<${CNAM}UnaryOp,T,D>(${CNAM}UnaryOp(${CNAM}<T>(),rhs));
+}
+
+
+//==============================================================================
+// $FNAM (SVector, unary)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<UnaryOp<${CNAM}<T>, SVector<T,D>, T>, T, D>
+ ${FNAM}(const SVector<T,D>& rhs) {
+ typedef UnaryOp<${CNAM}<T>, SVector<T,D>, T> ${CNAM}UnaryOp;
+
+ return Expr<${CNAM}UnaryOp,T,D>(${CNAM}UnaryOp(${CNAM}<T>(),rhs));
+}
+
+//==============================================================================
+// $FNAM (MatrixExpr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<${CNAM}<T>, Expr<A,T,D,D2>, T>, T, D, D2>
+ ${FNAM}(const Expr<A,T,D,D2>& rhs) {
+ typedef UnaryOp<${CNAM}<T>, Expr<A,T,D,D2>, T> ${CNAM}UnaryOp;
+
+ return Expr<${CNAM}UnaryOp,T,D,D2>(${CNAM}UnaryOp(${CNAM}<T>(),rhs));
+}
+
+
+//==============================================================================
+// $FNAM (SMatrix, unary)
+//==============================================================================
+template <class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<${CNAM}<T>, SMatrix<T,D,D2>, T>, T, D, D2>
+ ${FNAM}(const SMatrix<T,D,D2>& rhs) {
+ typedef UnaryOp<${CNAM}<T>, SMatrix<T,D,D2>, T> ${CNAM}UnaryOp;
+
+ return Expr<${CNAM}UnaryOp,T,D,D2>(${CNAM}UnaryOp(${CNAM}<T>(),rhs));
+}
+"
+ done
+
+echo "#endif"
+) > $OUTPUTFILE
+
+echo "$OUTPUTFILE generated"
diff --git a/smatrix/Module.mk b/smatrix/Module.mk
new file mode 100644
index 0000000000000000000000000000000000000000..279699d339f1418dbf7664030183e259bc57bdef
--- /dev/null
+++ b/smatrix/Module.mk
@@ -0,0 +1,86 @@
+# Module.mk for mathcore module
+# Copyright (c) 2000 Rene Brun and Fons Rademakers
+#
+# Author: Fons Rademakers, 20/6/2005
+
+MODDIR := smatrix
+MODDIRS := $(MODDIR)/src
+MODDIRI := $(MODDIR)/inc
+
+SMATRIXDIR := $(MODDIR)
+SMATRIXDIRS := $(SMATRIXDIR)/src
+SMATRIXDIRI := $(SMATRIXDIR)/inc/Math
+
+##### libSmatrix #####
+SMATRIXL := $(MODDIRI)/LinkDef.h
+#SMATRIXLINC :=
+SMATRIXDS := $(MODDIRS)/G__Smatrix.cxx
+SMATRIXDO := $(SMATRIXDS:.cxx=.o)
+SMATRIXDH := $(SMATRIXDS:.cxx=.h)
+
+SMATRIXDH1 := $(MODDIRI)/Math/SMatrix.h $(MODDIRI)/Math/SVector.h
+
+
+
+SMATRIXH1 := $(filter-out $(MODDIRI)/Math/LinkDef%, $(wildcard $(MODDIRI)/Math/*.h))
+SMATRIXH2 := $(filter-out $(MODDIRI)/Math/LinkDef%, $(wildcard $(MODDIRI)/Math/*.icc))
+SMATRIXH := $(SMATRIXH1) $(SMATRIXH2)
+SMATRIXS := $(filter-out $(MODDIRS)/G__%,$(wildcard $(MODDIRS)/*.cxx))
+SMATRIXO := $(SMATRIXS:.cxx=.o)
+
+SMATRIXDEP := $(SMATRIXO:.o=.d) $(SMATRIXDO:.o=.d)
+
+SMATRIXLIB := $(LPATH)/libSmatrix.$(SOEXT)
+
+# used in the main Makefile
+ALLHDRS += $(patsubst $(MODDIRI)/Math/%.h,include/Math/%.h,$(SMATRIXH1))
+ALLHDRS += $(patsubst $(MODDIRI)/Math/%.icc,include/Math/%.icc,$(SMATRIXH2))
+ALLLIBS += $(SMATRIXLIB)
+
+# include all dependency files
+INCLUDEFILES += $(SMATRIXDEP)
+
+##### local rules #####
+include/Math/%.h: $(SMATRIXDIRI)/%.h
+ cp $< $@
+
+include/Math/%.icc: $(SMATRIXDIRI)/%.icc
+ cp $< $@
+
+$(SMATRIXLIB): $(SMATRIXO) $(SMATRIXDO) $(MAINLIBS)
+ @$(MAKELIB) $(PLATFORM) $(LD) "$(LDFLAGS)" \
+ "$(SOFLAGS)" libSmatrix.$(SOEXT) $@ \
+ "$(SMATRIXO) $(SMATRIXDO)" \
+ "$(SMATRIXLIBEXTRA)"
+
+
+$(SMATRIXDS): $(SMATRIXDH1) $(SMATRIXL) $(SMATRIXLINC) $(ROOTCINTTMP)
+ @echo "Generating dictionary $@..."
+ @echo "for files $(SMATRIXDH1)"
+ $(ROOTCINTTMP) -f $@ -c $(SMATRIXDH1) $(SMATRIXL)
+# python reflex/python/genreflex/genreflex.py $(SMATRIXDIRS)/Dict.h -I$(SMATRIXDIRI) --selection_file=$(SMATRIXDIRS)/Selection.xml -o $(SMATRIXDIRS)/G__Smatrix.cxx
+
+
+$(SMATRIXDO): $(SMATRIXDS)
+ $(CXX) $(NOOPT) $(CXXFLAGS) -I. -I$(SMATRIXDIRI) -o $@ -c $<
+
+all-smatrix: $(SMATRIXLIB)
+
+map-smatrix: $(RLIBMAP)
+ $(RLIBMAP) -r $(ROOTMAP) -l $(SMATRIXLIB) \
+ -d $(SMATRIXLIBDEP) -c $(SMATRIXL) $(SMATRIXLINC)
+
+map:: map-smatrix
+
+clean-smatrix:
+ @rm -f $(SMATRIXO) $(SMATRIXDO)
+
+clean:: clean-smatrix
+
+distclean-smatrix: clean-smatrix
+ @rm -f $(SMATRIXDEP) $(SMATRIXDS) $(SMATRIXDH) $(SMATRIXLIB)
+ @rm -rf include/Math
+
+distclean:: distclean-smatrix
+
+
diff --git a/smatrix/inc/LinkDef.h b/smatrix/inc/LinkDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e50a4863f27de4a8785f5ef19876f00b86f3488
--- /dev/null
+++ b/smatrix/inc/LinkDef.h
@@ -0,0 +1,31 @@
+// @(#)root/smatrix:$Name: $:$Id: LinkDef.h,v 1.1 2005/10/28 15:58:38 moneta Exp $
+// Authors: L. Moneta 2005
+
+
+
+
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+//#pragma link C++ namespace tvmet;
+
+//#pragma link C++ typedef value_type;
+
+#pragma link C++ namespace ROOT::Math;
+
+//generate up to 5x5
+#pragma link C++ class ROOT::Math::SMatrix<double,2,2>+;
+#pragma link C++ class ROOT::Math::SMatrix<double,3,3>+;
+#pragma link C++ class ROOT::Math::SMatrix<double,4,4>+;
+#pragma link C++ class ROOT::Math::SMatrix<double,5,5>+;
+
+#pragma link C++ class ROOT::Math::SVector<double,2>+;
+#pragma link C++ class ROOT::Math::SVector<double,3>+;
+#pragma link C++ class ROOT::Math::SVector<double,4>+;
+#pragma link C++ class ROOT::Math::SVector<double,5>+;
+
+
+#endif
diff --git a/smatrix/inc/Math/BinaryOperators.h b/smatrix/inc/Math/BinaryOperators.h
new file mode 100644
index 0000000000000000000000000000000000000000..b9051e85906cfe836b97870464a9de9a662fadbf
--- /dev/null
+++ b/smatrix/inc/Math/BinaryOperators.h
@@ -0,0 +1,838 @@
+// @(#)root/smatrix:$Name: $:$Id: BinaryOperators.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_BinaryOperators
+#define ROOT_Math_BinaryOperators
+//======================================================
+//
+// ATTENTION: This file was automatically generated,
+// do not edit!
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+//======================================================
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+template <class T, unsigned int D> class SVector;
+template <class T, unsigned int D1, unsigned int D2> class SMatrix;
+
+
+//==============================================================================
+// AddOp
+//==============================================================================
+template <class T>
+class AddOp {
+public:
+ static inline T apply(const T& lhs, const T& rhs) {
+ return lhs + rhs;
+ }
+};
+
+
+//==============================================================================
+// operator+ (SVector, binary)
+//==============================================================================
+template < class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, SVector<T,D>, SVector<T,D>, T>, T, D>
+ operator+(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, SVector<T,D>, SVector<T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SVector, binary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, Expr<A,T,D>, SVector<T,D>, T>, T, D>
+ operator+(const Expr<A,T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<A,T,D>, SVector<T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SVector, binary)
+//==============================================================================
+template < class A, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, SVector<T,D>, Expr<A,T,D>, T>, T, D>
+ operator+(const SVector<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, SVector<T,D>, Expr<A,T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SVector, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, Expr<A,T,D>, Expr<B,T,D>, T>, T, D>
+ operator+(const Expr<A,T,D>& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<A,T,D>, Expr<B,T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, SVector<T,D>, Constant<A>, T>, T, D>
+ operator+(const SVector<T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<AddOp<T>, SVector<T,D>, Constant<A>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator+ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, Constant<A>, SVector<T,D>, T>, T, D>
+ operator+(const A& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, Constant<A>, SVector<T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator+ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, Expr<B,T,D>, Constant<A>, T>, T, D>
+ operator+(const Expr<B,T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<B,T,D>, Constant<A>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator+ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<AddOp<T>, Constant<A>, Expr<B,T,D>, T>, T, D>
+ operator+(const A& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<AddOp<T>, Constant<A>, Expr<B,T,D>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D>(AddOpBinOp(AddOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator+(const SMatrix<T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator+(const Expr<A,T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary)
+//==============================================================================
+template < class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T>, T, D, D2>
+ operator+(const SMatrix<T,D,D2>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator+(const Expr<A,T,D,D2>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, SMatrix<T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator+(const SMatrix<T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<AddOp<T>, SMatrix<T,D,D2>, Constant<A>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator+ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, Constant<A>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator+(const A& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, Constant<A>, SMatrix<T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator+ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, Expr<B,T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator+(const Expr<B,T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<AddOp<T>, Expr<B,T,D,D2>, Constant<A>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator+ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<AddOp<T>, Constant<A>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator+(const A& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<AddOp<T>, Constant<A>, Expr<B,T,D,D2>, T> AddOpBinOp;
+
+ return Expr<AddOpBinOp,T,D,D2>(AddOpBinOp(AddOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// MinOp
+//==============================================================================
+template <class T>
+class MinOp {
+public:
+ static inline T apply(const T& lhs, const T& rhs) {
+ return lhs - rhs;
+ }
+};
+
+
+//==============================================================================
+// operator- (SVector, binary)
+//==============================================================================
+template < class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, SVector<T,D>, SVector<T,D>, T>, T, D>
+ operator-(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, SVector<T,D>, SVector<T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, binary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, Expr<A,T,D>, SVector<T,D>, T>, T, D>
+ operator-(const Expr<A,T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<A,T,D>, SVector<T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, binary)
+//==============================================================================
+template < class A, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, SVector<T,D>, Expr<A,T,D>, T>, T, D>
+ operator-(const SVector<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, SVector<T,D>, Expr<A,T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, Expr<A,T,D>, Expr<B,T,D>, T>, T, D>
+ operator-(const Expr<A,T,D>& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<A,T,D>, Expr<B,T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, SVector<T,D>, Constant<A>, T>, T, D>
+ operator-(const SVector<T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<MinOp<T>, SVector<T,D>, Constant<A>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator- (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, Constant<A>, SVector<T,D>, T>, T, D>
+ operator-(const A& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, Constant<A>, SVector<T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, Expr<B,T,D>, Constant<A>, T>, T, D>
+ operator-(const Expr<B,T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<B,T,D>, Constant<A>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator- (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MinOp<T>, Constant<A>, Expr<B,T,D>, T>, T, D>
+ operator-(const A& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<MinOp<T>, Constant<A>, Expr<B,T,D>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D>(MinOpBinOp(MinOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator-(const SMatrix<T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator-(const Expr<A,T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary)
+//==============================================================================
+template < class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T>, T, D, D2>
+ operator-(const SMatrix<T,D,D2>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator-(const Expr<A,T,D,D2>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, SMatrix<T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator-(const SMatrix<T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<MinOp<T>, SMatrix<T,D,D2>, Constant<A>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator- (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, Constant<A>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator-(const A& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, Constant<A>, SMatrix<T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, Expr<B,T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator-(const Expr<B,T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<MinOp<T>, Expr<B,T,D,D2>, Constant<A>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator- (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MinOp<T>, Constant<A>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator-(const A& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<MinOp<T>, Constant<A>, Expr<B,T,D,D2>, T> MinOpBinOp;
+
+ return Expr<MinOpBinOp,T,D,D2>(MinOpBinOp(MinOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// MulOp
+//==============================================================================
+template <class T>
+class MulOp {
+public:
+ static inline T apply(const T& lhs, const T& rhs) {
+ return lhs * rhs;
+ }
+};
+
+
+//==============================================================================
+// operator* (SVector, binary)
+//==============================================================================
+template < class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, SVector<T,D>, SVector<T,D>, T>, T, D>
+ operator*(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, SVector<T,D>, SVector<T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator* (SVector, binary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, Expr<A,T,D>, SVector<T,D>, T>, T, D>
+ operator*(const Expr<A,T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<A,T,D>, SVector<T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator* (SVector, binary)
+//==============================================================================
+template < class A, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, SVector<T,D>, Expr<A,T,D>, T>, T, D>
+ operator*(const SVector<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, SVector<T,D>, Expr<A,T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator* (SVector, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, Expr<A,T,D>, Expr<B,T,D>, T>, T, D>
+ operator*(const Expr<A,T,D>& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<A,T,D>, Expr<B,T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator* (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, SVector<T,D>, Constant<A>, T>, T, D>
+ operator*(const SVector<T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<MulOp<T>, SVector<T,D>, Constant<A>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator* (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, Constant<A>, SVector<T,D>, T>, T, D>
+ operator*(const A& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, Constant<A>, SVector<T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator* (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, Expr<B,T,D>, Constant<A>, T>, T, D>
+ operator*(const Expr<B,T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<B,T,D>, Constant<A>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator* (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<MulOp<T>, Constant<A>, Expr<B,T,D>, T>, T, D>
+ operator*(const A& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<MulOp<T>, Constant<A>, Expr<B,T,D>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D>(MulOpBinOp(MulOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// times (SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ times(const SMatrix<T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// times (SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ times(const Expr<A,T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// times (SMatrix, binary)
+//==============================================================================
+template < class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T>, T, D, D2>
+ times(const SMatrix<T,D,D2>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// times (SMatrix, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T>, T, D, D2>
+ times(const Expr<A,T,D,D2>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator* (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, SMatrix<T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator*(const SMatrix<T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<MulOp<T>, SMatrix<T,D,D2>, Constant<A>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator* (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, Constant<A>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator*(const A& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, Constant<A>, SMatrix<T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator* (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, Expr<B,T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator*(const Expr<B,T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<MulOp<T>, Expr<B,T,D,D2>, Constant<A>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator* (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<MulOp<T>, Constant<A>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator*(const A& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<MulOp<T>, Constant<A>, Expr<B,T,D,D2>, T> MulOpBinOp;
+
+ return Expr<MulOpBinOp,T,D,D2>(MulOpBinOp(MulOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// DivOp
+//==============================================================================
+template <class T>
+class DivOp {
+public:
+ static inline T apply(const T& lhs, const T& rhs) {
+ return lhs / rhs;
+ }
+};
+
+
+//==============================================================================
+// operator/ (SVector, binary)
+//==============================================================================
+template < class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, SVector<T,D>, SVector<T,D>, T>, T, D>
+ operator/(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, SVector<T,D>, SVector<T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SVector, binary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Expr<A,T,D>, SVector<T,D>, T>, T, D>
+ operator/(const Expr<A,T,D>& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<A,T,D>, SVector<T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SVector, binary)
+//==============================================================================
+template < class A, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, SVector<T,D>, Expr<A,T,D>, T>, T, D>
+ operator/(const SVector<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, SVector<T,D>, Expr<A,T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SVector, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Expr<A,T,D>, Expr<B,T,D>, T>, T, D>
+ operator/(const Expr<A,T,D>& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<A,T,D>, Expr<B,T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, SVector<T,D>, Constant<A>, T>, T, D>
+ operator/(const SVector<T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<DivOp<T>, SVector<T,D>, Constant<A>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator/ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Constant<A>, SVector<T,D>, T>, T, D>
+ operator/(const A& lhs, const SVector<T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, Constant<A>, SVector<T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator/ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Expr<B,T,D>, Constant<A>, T>, T, D>
+ operator/(const Expr<B,T,D>& lhs, const A& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<B,T,D>, Constant<A>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator/ (SVector, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Constant<A>, Expr<B,T,D>, T>, T, D>
+ operator/(const A& lhs, const Expr<B,T,D>& rhs) {
+ typedef BinaryOp<DivOp<T>, Constant<A>, Expr<B,T,D>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator/(const SMatrix<T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, SMatrix<T,D,D2>, SMatrix<T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator/(const Expr<A,T,D,D2>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<A,T,D,D2>, SMatrix<T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary)
+//==============================================================================
+template < class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T>, T, D, D2>
+ operator/(const SMatrix<T,D,D2>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, SMatrix<T,D,D2>, Expr<A,T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator/(const Expr<A,T,D,D2>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<A,T,D,D2>, Expr<B,T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, SMatrix<T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator/(const SMatrix<T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<DivOp<T>, SMatrix<T,D,D2>, Constant<A>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator/ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, Constant<A>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator/(const A& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, Constant<A>, SMatrix<T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+//==============================================================================
+// operator/ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, Expr<B,T,D,D2>, Constant<A>, T>, T, D, D2>
+ operator/(const Expr<B,T,D,D2>& lhs, const A& rhs) {
+ typedef BinaryOp<DivOp<T>, Expr<B,T,D,D2>, Constant<A>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),lhs,Constant<A>(rhs)));
+}
+
+//==============================================================================
+// operator/ (SMatrix, binary, Constant)
+//==============================================================================
+template <class A, class B, class T, unsigned int D, unsigned int D2>
+inline Expr<BinaryOp<DivOp<T>, Constant<A>, Expr<B,T,D,D2>, T>, T, D, D2>
+ operator/(const A& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef BinaryOp<DivOp<T>, Constant<A>, Expr<B,T,D,D2>, T> DivOpBinOp;
+
+ return Expr<DivOpBinOp,T,D,D2>(DivOpBinOp(DivOp<T>(),Constant<A>(lhs),rhs));
+}
+
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+#endif /*ROOT_Math_BinaryOperators */
diff --git a/smatrix/inc/Math/Dfact.h b/smatrix/inc/Math/Dfact.h
new file mode 100644
index 0000000000000000000000000000000000000000..1389519b0d31da05357ec75acbd948faf2b33026
--- /dev/null
+++ b/smatrix/inc/Math/Dfact.h
@@ -0,0 +1,152 @@
+// @(#)root/smatrix:$Name: $:$Id: Dfact.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dfact
+#define ROOT_Math_Dfact
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 02. Apr 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Determinant of a square matrix
+// Code was taken from CERNLIB::kernlib dfact function, translated
+// from FORTRAN to C++ and optimized.
+//
+// changes:
+// 02 Apr 2001 (TG) creation
+//
+// ********************************************************************
+
+#include <cmath>
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+/** Dfact.
+ Function to compute the determinant from a square matrix ($\det(A)$) of
+ dimension $idim$ and order $n$.
+
+ @author T. Glebe
+*/
+template <class Matrix, unsigned int n, unsigned int idim>
+bool Dfact(Matrix& rhs, typename Matrix::value_type& det) {
+
+#ifdef XXX
+ if (idim < n || n <= 0) {
+ return false;
+ }
+#endif
+
+
+ /* Initialized data */
+ // const typename Matrix::value_type* A = rhs.Array();
+ typename Matrix::value_type* a = rhs.Array();
+
+ /* Local variables */
+ static unsigned int nxch, i, j, k, l;
+ static typename Matrix::value_type p, q, tf;
+
+ /* Parameter adjustments */
+ a -= idim + 1;
+
+ /* Function Body */
+
+ // fact.inc
+
+ nxch = 0;
+ det = 1.;
+ for (j = 1; j <= n; ++j) {
+ const unsigned int ji = j * idim;
+ const unsigned int jj = j + ji;
+
+ k = j;
+ p = std::fabs(a[jj]);
+
+ if (j != n) {
+ for (i = j + 1; i <= n; ++i) {
+ q = std::fabs(a[i + ji]);
+ if (q > p) {
+ k = i;
+ p = q;
+ }
+ } // for i
+ if (k != j) {
+ for (l = 1; l <= n; ++l) {
+ const unsigned int li = l*idim;
+ const unsigned int jli = j + li;
+ const unsigned int kli = k + li;
+ tf = a[jli];
+ a[jli] = a[kli];
+ a[kli] = tf;
+ } // for l
+ ++nxch;
+ } // if k != j
+ } // if j!=n
+
+ if (p <= 0.) {
+ det = 0;
+ return false;
+ }
+
+ det *= a[jj];
+#ifdef XXX
+ t = std::fabs(det);
+ if (t < 1e-19 || t > 1e19) {
+ det = 0;
+ return false;
+ }
+#endif
+
+ a[jj] = 1. / a[jj];
+ if (j == n) {
+ continue;
+ }
+
+ const unsigned int jm1 = j - 1;
+ const unsigned int jpi = (j + 1) * idim;
+ const unsigned int jjpi = j + jpi;
+
+ for (k = j + 1; k <= n; ++k) {
+ const unsigned int ki = k * idim;
+ const unsigned int jki = j + ki;
+ const unsigned int kji = k + jpi;
+ if (j != 1) {
+ for (i = 1; i <= jm1; ++i) {
+ const unsigned int ii = i * idim;
+ a[jki] -= a[i + ki] * a[j + ii];
+ a[kji] -= a[i + jpi] * a[k + ii];
+ } // for i
+ }
+ a[jki] *= a[jj];
+ a[kji] -= a[jjpi] * a[k + ji];
+ } // for k
+ } // for j
+
+ if (nxch % 2 != 0) {
+ det = -(det);
+ }
+ return true;
+} // end of Dfact
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Dfact */
diff --git a/smatrix/inc/Math/Dfactir.h b/smatrix/inc/Math/Dfactir.h
new file mode 100644
index 0000000000000000000000000000000000000000..be4b6b354a476ae48b04697488b1a99b51ce62ae
--- /dev/null
+++ b/smatrix/inc/Math/Dfactir.h
@@ -0,0 +1,155 @@
+// @(#)root/smatrix:$Name: $:$Id: Dfactir.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dfactir
+#define ROOT_Math_Dfactir
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 02. Apr 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Determinant of a square matrix, needed for Dfinv()
+// Code was taken from CERNLIB::kernlib dfact function, translated
+// from FORTRAN to C++ and optimized.
+//
+// changes:
+// 02 Apr 2001 (TG) creation
+//
+// ********************************************************************
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+/** Dfactir.
+ Function to compute the determinant from a square matrix ($\det(A)$) of
+ dimension $idim$ and order $n$. A working area $ir$ is returned which is
+ needed by the Dfinv routine.
+
+ @author T. Glebe
+*/
+template <class Matrix, unsigned int n, unsigned int idim>
+bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
+ // int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
+{
+
+#ifdef XXX
+ if (idim < n || n <= 0) {
+ return false;
+ }
+#endif
+
+
+ /* Initialized data */
+ typename Matrix::value_type* a = rhs.Array();
+
+ /* Local variables */
+ static unsigned int nxch, i, j, k, l;
+ static typename Matrix::value_type p, q, tf;
+
+ /* Parameter adjustments */
+ a -= idim + 1;
+ --ir;
+
+ /* Function Body */
+
+ // fact.inc
+ nxch = 0;
+ det = 1.;
+ for (j = 1; j <= n; ++j) {
+ const unsigned int ji = j * idim;
+ const unsigned int jj = j + ji;
+
+ k = j;
+ p = std::fabs(a[jj]);
+
+ if (j != n) {
+ for (i = j + 1; i <= n; ++i) {
+ q = std::fabs(a[i + ji]);
+ if (q > p) {
+ k = i;
+ p = q;
+ }
+ } // for i
+
+ if (k != j) {
+ for (l = 1; l <= n; ++l) {
+ const unsigned int li = l*idim;
+ const unsigned int jli = j + li;
+ const unsigned int kli = k + li;
+ tf = a[jli];
+ a[jli] = a[kli];
+ a[kli] = tf;
+ } // for l
+ ++nxch;
+ ir[nxch] = (j << 12) + k;
+ } // if k != j
+ } // if j!=n
+
+ if (p <= 0.) {
+ det = 0;
+ return false;
+ }
+
+ det *= a[jj];
+#ifdef XXX
+ t = std::fabs(det);
+ if (t < 1e-19 || t > 1e19) {
+ det = 0;
+ return false;
+ }
+#endif
+
+ a[jj] = 1. / a[jj];
+ if (j == n) {
+ continue;
+ }
+
+ const unsigned int jm1 = j - 1;
+ const unsigned int jpi = (j + 1) * idim;
+ const unsigned int jjpi = j + jpi;
+
+ for (k = j + 1; k <= n; ++k) {
+ const unsigned int ki = k * idim;
+ const unsigned int jki = j + ki;
+ const unsigned int kji = k + jpi;
+ if (j != 1) {
+ for (i = 1; i <= jm1; ++i) {
+ const unsigned int ii = i * idim;
+ a[jki] -= a[i + ki] * a[j + ii];
+ a[kji] -= a[i + jpi] * a[k + ii];
+ } // for i
+ }
+ a[jki] *= a[jj];
+ a[kji] -= a[jjpi] * a[k + ji];
+ } // for k
+ } // for j
+
+ if (nxch % 2 != 0) {
+ det = -(det);
+ }
+ ir[n] = nxch;
+ return true;
+} // end of Dfact
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Dfactir */
diff --git a/smatrix/inc/Math/Dfinv.h b/smatrix/inc/Math/Dfinv.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6cf0277f217cc307848fa575df85d207bc90c0e
--- /dev/null
+++ b/smatrix/inc/Math/Dfinv.h
@@ -0,0 +1,146 @@
+// @(#)root/smatrix:$Name: $:$Id: Dfinv.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dfinv
+#define ROOT_Math_Dfinv
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 03. Apr 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Matrix inversion
+// Code was taken from CERNLIB::kernlib dfinv function, translated
+// from FORTRAN to C++ and optimized.
+//
+// changes:
+// 03 Apr 2001 (TG) creation
+//
+// ********************************************************************
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+
+/** Dfinv.
+ Function to compute the inverse of a square matrix ($A^{-1}$) of
+ dimension $idim$ and order $n$. The routine Dfactir must be called
+ before Dfinv!
+
+ @author T. Glebe
+*/
+template <class Matrix, unsigned int n, unsigned int idim>
+bool Dfinv(Matrix& rhs, unsigned int* ir) {
+#ifdef XXX
+ if (idim < n || n <= 0 || n==1) {
+ return false;
+ }
+#endif
+
+ typename Matrix::value_type* a = rhs.Array();
+
+ /* Local variables */
+ static unsigned int nxch, i, j, k, m, ij;
+ static unsigned int im2, nm1, nmi;
+ static typename Matrix::value_type s31, s34, ti;
+
+ /* Parameter adjustments */
+ a -= idim + 1;
+ --ir;
+
+ /* Function Body */
+
+ /* finv.inc */
+
+ a[idim + 2] = -a[(idim << 1) + 2] * (a[idim + 1] * a[idim + 2]);
+ a[(idim << 1) + 1] = -a[(idim << 1) + 1];
+
+ if (n != 2) {
+ for (i = 3; i <= n; ++i) {
+ const unsigned int ii = i * idim;
+ const unsigned int iii = i + ii;
+ const unsigned int imi = ii - idim;
+ const unsigned int iimi = i + imi;
+ im2 = i - 2;
+ for (j = 1; j <= im2; ++j) {
+ const unsigned int ji = j * idim;
+ const unsigned int jii = j + ii;
+ s31 = 0.;
+ for (k = j; k <= im2; ++k) {
+ s31 += a[k + ji] * a[i + k * idim];
+ a[jii] += a[j + (k + 1) * idim] * a[k + 1 + ii];
+ } // for k
+ a[i + ji] = -a[iii] * (a[i - 1 + ji] * a[iimi] + s31);
+ a[jii] *= -1;
+ } // for j
+ a[iimi] = -a[iii] * (a[i - 1 + imi] * a[iimi]);
+ a[i - 1 + ii] *= -1;
+ } // for i
+ } // if n!=2
+
+ nm1 = n - 1;
+ for (i = 1; i <= nm1; ++i) {
+ const unsigned int ii = i * idim;
+ nmi = n - i;
+ for (j = 1; j <= i; ++j) {
+ const unsigned int ji = j * idim;
+ const unsigned int iji = i + ji;
+ for (k = 1; k <= nmi; ++k) {
+ a[iji] += a[i + k + ji] * a[i + (i + k) * idim];
+ } // for k
+ } // for j
+
+ for (j = 1; j <= nmi; ++j) {
+ const unsigned int ji = j * idim;
+ s34 = 0.;
+ for (k = j; k <= nmi; ++k) {
+ s34 += a[i + k + ii + ji] * a[i + (i + k) * idim];
+ } // for k
+ a[i + ii + ji] = s34;
+ } // for j
+ } // for i
+
+ nxch = ir[n];
+ if (nxch == 0) {
+ return true;
+ }
+
+ for (m = 1; m <= nxch; ++m) {
+ k = nxch - m + 1;
+ ij = ir[k];
+ i = ij / 4096;
+ j = ij % 4096;
+ const unsigned int ii = i * idim;
+ const unsigned int ji = j * idim;
+ for (k = 1; k <= n; ++k) {
+ ti = a[k + ii];
+ a[k + ii] = a[k + ji];
+ a[k + ji] = ti;
+ } // for k
+ } // for m
+
+ return true;
+} // Dfinv
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Dfinv */
diff --git a/smatrix/inc/Math/Dinv.h b/smatrix/inc/Math/Dinv.h
new file mode 100644
index 0000000000000000000000000000000000000000..8db4a0d97878dba3b2c3743fe9b46ef881305fb4
--- /dev/null
+++ b/smatrix/inc/Math/Dinv.h
@@ -0,0 +1,246 @@
+// @(#)root/smatrix:$Name: $:$Id: Dinv.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dinv
+#define ROOT_Math_Dinv
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 03. Apr 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: square Matrix inversion
+// Code was taken from CERNLIB::kernlib dfinv function, translated
+// from FORTRAN to C++ and optimized.
+// n: Order of the square matrix
+// idim: First dimension of array A
+//
+// changes:
+// 03 Apr 2001 (TG) creation
+//
+// ********************************************************************
+#include "Math/Dfactir.h"
+#include "Math/Dfinv.h"
+
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+/** Inverter.
+ Class to specialize calls to Dinv. Dinv computes the inverse of a square
+ matrix if dimension $idim$ and order $n$. The content of the matrix will be
+ replaced by its inverse. In case the inversion fails, the matrix content is
+ destroyed. Invert specializes Dinv by the matrix order. E.g. if the order
+ of the matrix is two, the routine Invert<2> is called which implements
+ Cramers rule.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Inverter class
+//==============================================================================
+template <unsigned int idim, unsigned int n = idim>
+class Inverter {
+public:
+ ///
+ template <class Matrix>
+ static bool Dinv(Matrix& rhs) {
+
+#ifdef XXX
+ if (n < 1 || n > idim) {
+ return false;
+ }
+#endif
+
+ /* Initialized data */
+ static unsigned int work[n];
+ for(unsigned int i=0; i<n; ++i) work[i] = 0;
+
+ static typename Matrix::value_type det = 0;
+
+ /* Function Body */
+
+ /* N.GT.3 CASES. FACTORIZE MATRIX AND INVERT. */
+ if (Dfactir<Matrix,n,idim>(rhs,det,work) == false) {
+ std::cerr << "Dfactir failed!!" << std::endl;
+ return false;
+ }
+ return Dfinv<Matrix,n,idim>(rhs,work);
+ } // Dinv
+}; // class Inverter
+
+
+/** Inverter<0>.
+ In case of zero order, do nothing.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Invert<0>
+//==============================================================================
+template <>
+class Inverter<0> {
+public:
+ ///
+ template <class Matrix>
+ inline static bool Dinv(Matrix& rhs) { return true; }
+};
+
+
+/** Invert<1>.
+ $1\times1$ (sub-)matrix. $a_{11} \to 1/a_{11}$
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Invert<1>
+//==============================================================================
+template <>
+class Inverter<1> {
+public:
+ ///
+ template <class Matrix>
+ static bool Dinv(Matrix& rhs) {
+ typename Matrix::value_type* a = rhs.Array();
+
+ if (a[0] == 0.) {
+ return false;
+ }
+ a[0] = 1. / a[0];
+ return true;
+ }
+};
+
+
+/** Invert<2>.
+ $2\times2$ (sub-)matrix. Use Cramers rule.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Invert<2>: Cramers rule
+//==============================================================================
+template <>
+class Inverter<2> {
+public:
+ ///
+ template <class Matrix>
+ static bool Dinv(Matrix& rhs) {
+
+ typename Matrix::value_type* a = rhs.Array();
+ typename Matrix::value_type det = a[0] * a[3] - a[2] * a[1];
+
+ if (det == 0.) { return false; }
+
+ typename Matrix::value_type s = 1. / det;
+ typename Matrix::value_type c11 = s * a[3];
+
+ a[2] = -s * a[2];
+ a[1] = -s * a[1];
+ a[3] = s * a[0];
+ a[0] = c11;
+ return true;
+ }
+};
+
+
+/** Invert<3>.
+ $3\times3$ (sub-)matrix. Use pivotisation.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Inverter<3>
+//==============================================================================
+template <>
+class Inverter<3> {
+public:
+ ///
+ template <class Matrix>
+ static bool Dinv(Matrix& rhs) {
+
+ typename Matrix::value_type* a = rhs.Array();
+
+ static typename Matrix::value_type t1, t2, t3, temp, s;
+ static typename Matrix::value_type c11, c12, c13, c21, c22, c23, c31, c32, c33, det;
+
+
+ /* COMPUTE COFACTORS. */
+ c11 = a[4] * a[8] - a[7] * a[5];
+ c12 = a[7] * a[2] - a[1] * a[8];
+ c13 = a[1] * a[5] - a[4] * a[2];
+ c21 = a[5] * a[6] - a[8] * a[3];
+ c22 = a[8] * a[0] - a[2] * a[6];
+ c23 = a[2] * a[3] - a[5] * a[0];
+ c31 = a[3] * a[7] - a[6] * a[4];
+ c32 = a[6] * a[1] - a[0] * a[7];
+ c33 = a[0] * a[4] - a[3] * a[1];
+
+ t1 = std::fabs(a[0]);
+ t2 = std::fabs(a[1]);
+ t3 = std::fabs(a[2]);
+
+ /* (SET TEMP=PIVOT AND DET=PIVOT*DET.) */
+ if(t1 < t2) {
+ if (t3 < t2) {
+ /* (PIVOT IS A21) */
+ temp = a[1];
+ det = c13 * c32 - c12 * c33;
+ } else {
+ /* (PIVOT IS A31) */
+ temp = a[2];
+ det = c23 * c12 - c22 * c13;
+ }
+ } else {
+ if(t3 < t1) {
+ /* (PIVOT IS A11) */
+ temp = a[0];
+ det = c22 * c33 - c23 * c32;
+ } else {
+ /* (PIVOT IS A31) */
+ temp = a[2];
+ det = c23 * c12 - c22 * c13;
+ }
+ }
+
+ if (det == 0.) {
+ return false;
+ }
+
+ /* SET ELEMENTS OF INVERSE IN A. */
+ s = temp / det;
+ a[0] = s * c11;
+ a[3] = s * c21;
+ a[6] = s * c31;
+ a[1] = s * c12;
+ a[4] = s * c22;
+ a[7] = s * c32;
+ a[2] = s * c13;
+ a[5] = s * c23;
+ a[8] = s * c33;
+ return true;
+ }
+};
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Dinv */
diff --git a/smatrix/inc/Math/Dsfact.h b/smatrix/inc/Math/Dsfact.h
new file mode 100644
index 0000000000000000000000000000000000000000..abd6280611591d1a16db96b5b79d5914da92daaa
--- /dev/null
+++ b/smatrix/inc/Math/Dsfact.h
@@ -0,0 +1,112 @@
+// @(#)root/smatrix:$Name: $:$Id: Dsfact.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dsfact
+#define ROOT_Math_Dsfact
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 22. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Determinant of a symmetric, positive definite matrix.
+// Code was taken from CERNLIB::kernlib dsfact function, translated
+// from FORTRAN to C++ and optimized.
+//
+// changes:
+// 22 Mar 2001 (TG) creation
+// 18 Apr 2001 (TG) removed internal copying of array.
+//
+// ********************************************************************
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+
+/** Dsfact.
+ Compute determinant of a symmetric, positive definite matrix of dimension
+ $idim$ and order $n$.
+
+ @author T. Glebe
+*/
+template <class Matrix, unsigned int n, unsigned int idim>
+bool Dsfact(Matrix& rhs, typename Matrix::value_type& det) {
+
+#ifdef XXX
+ /* Function Body */
+ if (idim < n || n <= 0) {
+ return false;
+ }
+#endif
+
+ typename Matrix::value_type* a = rhs.Array();
+
+#ifdef XXX
+ const typename Matrix::value_type* A = rhs.Array();
+ typename Matrix::value_type array[Matrix::size];
+ typename Matrix::value_type* a = array;
+
+ // copy contents of matrix to working place
+ for(unsigned int i=0; i<Matrix::size; ++i) {
+ array[i] = A[i];
+ }
+#endif
+
+ /* Local variables */
+ static unsigned int i, j, l;
+
+ /* Parameter adjustments */
+ a -= idim + 1;
+
+ /* sfactd.inc */
+ det = 1.;
+ for (j = 1; j <= n; ++j) {
+ const unsigned int ji = j * idim;
+ const unsigned int jj = j + ji;
+
+ if (a[jj] <= 0.) {
+ det = 0.;
+ return false;
+ }
+
+ const unsigned int jp1 = j + 1;
+ const unsigned int jpi = jp1 * idim;
+
+ det *= a[jj];
+ a[jj] = 1. / a[jj];
+
+ for (l = jp1; l <= n; ++l) {
+ a[j + l * idim] = a[jj] * a[l + ji];
+
+ const unsigned int lj = l + jpi;
+
+ for (i = 1; i <= j; ++i) {
+ a[lj] -= a[l + i * idim] * a[i + jpi];
+ } // for i
+ } // for l
+ } // for j
+
+ return true;
+} // end of Dsfact
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+#endif /* ROOT_Math_Dsfact */
+
diff --git a/smatrix/inc/Math/Dsinv.h b/smatrix/inc/Math/Dsinv.h
new file mode 100644
index 0000000000000000000000000000000000000000..e4b4a13e32ed31fd00893464ded5f926ecb8c066
--- /dev/null
+++ b/smatrix/inc/Math/Dsinv.h
@@ -0,0 +1,138 @@
+// @(#)root/smatrix:$Name: $:$Id: Dsinv.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Dsinv
+#define ROOT_Math_Dsinv
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 22. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Inversion of a symmetric, positive definite matrix.
+// Code was taken from CERNLIB::kernlib dsinv function, translated
+// from FORTRAN to C++ and optimized.
+//
+// changes:
+// 22 Mar 2001 (TG) creation
+//
+// ********************************************************************
+
+
+namespace ROOT {
+
+ namespace Math {
+
+/** Dsinv.
+ Compute inverse of a symmetric, positive definite matrix of dimension
+ $idim$ and order $n$.
+
+ @author T. Glebe
+*/
+template <class T, int n, int idim>
+bool Dsinv(T* a) {
+
+ /* Local variables */
+ static int i, j, k, l;
+ static T s31, s32;
+ static int jm1, jp1;
+
+ /* Parameter adjustments */
+ a -= idim + 1;
+
+ /* Function Body */
+ if (idim < n || n <= 1) {
+ return false;
+ }
+
+ /* sfact.inc */
+ for (j = 1; j <= n; ++j) {
+ const int ja = j * idim;
+ const int jj = j + ja;
+ const int ja1 = ja + idim;
+
+ if (a[jj] <= 0.) { return false; }
+ a[jj] = 1. / a[jj];
+ if (j == n) { break; }
+
+ for (l = j + 1; l <= n; ++l) {
+ a[j + l * idim] = a[jj] * a[l + ja];
+ const int lj = l + ja1;
+ for (i = 1; i <= j; ++i) {
+ a[lj] -= a[l + i * idim] * a[i + ja1];
+ }
+ }
+ }
+
+ /* sfinv.inc */
+ // idim << 1 is equal to idim * 2
+ // compiler will compute the arguments!
+ a[(idim << 1) + 1] = -a[(idim << 1) + 1];
+ a[idim + 2] = a[(idim << 1) + 1] * a[(idim << 1) + 2];
+
+ if(n > 2) {
+
+ for (j = 3; j <= n; ++j) {
+ const int jm2 = j - 2;
+ const int ja = j * idim;
+ const int jj = j + ja;
+ const int j1 = j - 1 + ja;
+
+ for (k = 1; k <= jm2; ++k) {
+ s31 = a[k + ja];
+
+ for (i = k; i <= jm2; ++i) {
+ s31 += a[k + (i + 1) * idim] * a[i + 1 + ja];
+ } // for i
+ a[k + ja] = -s31;
+ a[j + k * idim] = -s31 * a[jj];
+ } // for k
+ a[j1] *= -1;
+ // a[j1] = -a[j1];
+ a[jj - idim] = a[j1] * a[jj];
+ } // for j
+ } // if (n>2)
+
+ j = 1;
+ do {
+ const int jad = j * idim;
+ const int jj = j + jad;
+
+ jp1 = j + 1;
+ for (i = jp1; i <= n; ++i) {
+ a[jj] += a[j + i * idim] * a[i + jad];
+ } // for i
+
+ jm1 = j;
+ j = jp1;
+ const int ja = j * idim;
+
+ for (k = 1; k <= jm1; ++k) {
+ s32 = 0.;
+ for (i = j; i <= n; ++i) {
+ s32 += a[k + i * idim] * a[i + ja];
+ } // for i
+ a[k + ja] = a[j + k * idim] = s32;
+ } // for k
+ } while(j < n);
+
+ return true;
+} // end of Dsinv
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+#endif /* ROOT_Math_Dsinv */
diff --git a/smatrix/inc/Math/Expression.h b/smatrix/inc/Math/Expression.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc29e7c5d1a3b7f1d613865b848c2549ee203bef
--- /dev/null
+++ b/smatrix/inc/Math/Expression.h
@@ -0,0 +1,240 @@
+// @(#)root/smatrix:$Name: $:$Id: Expression.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Expression
+#define ROOT_Math_Expression
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 19. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Expression Template Elements for SVector
+//
+// changes:
+// 19 Mar 2001 (TG) creation
+// 20 Mar 2001 (TG) added rows(), cols() to Expr
+// 21 Mar 2001 (TG) added Expr::value_type
+// 11 Apr 2001 (TG) rows(), cols() replaced by rows, cols
+// 10 Okt 2001 (TG) added print() and operator<<() for Expr class
+//
+// ********************************************************************
+
+/** Expr.
+ An Expression wrapper class.
+
+ @memo Expr.
+ @author T. Glebe
+*/
+//==============================================================================
+// Expr: class representing SVector expressions
+//=============================================================================
+
+// originally references where used as data members of operation wrappers
+// problem appear on Windows because original ref goes out of scope use then by value by default (like in tvmet)
+//#define SMATRIX_USE_REFERENCES
+
+
+#include <iomanip>
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+template <class ExprType, class T, unsigned int D, unsigned int D2 = 0>
+class Expr {
+public:
+ typedef T value_type;
+
+ ///
+ Expr(const ExprType& rhs) :
+ rhs_(rhs) {}
+
+ ///
+ ~Expr() {}
+
+ ///
+ inline T apply(unsigned int i) const {
+ return rhs_.apply(i);
+ }
+
+#ifdef OLD_IMPL
+ ///
+ static const unsigned int rows = D;
+ ///
+ static const unsigned int cols = D2;
+#else
+ // use enumerations
+ enum {
+ ///
+ rows = D,
+ ///
+ cols = D2,
+ };
+#endif
+
+ /// used by operator<<()
+ std::ostream& print(std::ostream& os) const {
+ os.setf(std::ios::right,std::ios::adjustfield);
+
+ if(D2 == 0) {
+ unsigned int i=0;
+ for(; i<D-1; ++i) {
+ os << apply(i) << ", ";
+ }
+ os << apply(i);
+ } else {
+ os << "[ ";
+ for (unsigned int i=0; i < D; ++i) {
+ for (unsigned int j=0; j < D2; ++j) {
+ os << std::setw(12) << apply(i*D2+j);
+ if ((!((j+1)%12)) && (j < D2-1))
+ os << std::endl << " ...";
+ }
+ if (i != D - 1)
+ os << std::endl << " ";
+ }
+ os << " ]";
+ } // if D2==0
+
+ return os;
+ }
+
+private:
+ ExprType rhs_; // cannot be a reference!
+};
+
+//==============================================================================
+// operator<<
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline std::ostream& operator<<(std::ostream& os, const Expr<A,T,D1,D2>& rhs) {
+ return rhs.print(os);
+}
+
+/** BinaryOp.
+ A class representing binary operators in the parse tree.
+
+ @memo BinaryOp
+ @author T. Glebe
+*/
+//==============================================================================
+// BinaryOp
+//==============================================================================
+template <class Operator, class LHS, class RHS, class T>
+class BinaryOp {
+public:
+ ///
+ BinaryOp( Operator /* op */, const LHS& lhs, const RHS& rhs) :
+ lhs_(lhs), rhs_(rhs) {}
+
+ ///
+ ~BinaryOp() {}
+
+ ///
+ inline T apply(unsigned int i) const {
+ return Operator::apply(lhs_.apply(i), rhs_.apply(i));
+ }
+
+protected:
+
+#ifdef SMATRIX_USE_REFERENCES
+ const LHS& lhs_;
+ const RHS& rhs_;
+#else
+ const LHS lhs_;
+ const RHS rhs_;
+#endif
+
+};
+
+
+/** UnaryOp.
+ A class representing unary operators in the parse tree.
+
+ @memo UnaryOp
+ @author T. Glebe
+*/
+//==============================================================================
+// UnaryOp
+//==============================================================================
+template <class Operator, class RHS, class T>
+class UnaryOp {
+public:
+ ///
+ UnaryOp( Operator /* op */ , const RHS& rhs) :
+ rhs_(rhs) {}
+
+ ///
+ ~UnaryOp() {}
+
+ ///
+ inline T apply(unsigned int i) const {
+ return Operator::apply(rhs_.apply(i));
+ }
+
+protected:
+
+#ifdef SMATRIX_USE_REFERENCES
+ const RHS& rhs_;
+#else
+ const RHS rhs_;
+#endif
+
+};
+
+
+/** Constant.
+ A class representing constant expressions (literals) in the parse tree.
+
+ @memo Constant
+ @author T. Glebe
+*/
+//==============================================================================
+// Constant
+//==============================================================================
+template <class T>
+class Constant {
+public:
+ ///
+ Constant( const T& rhs ) :
+ rhs_(rhs) {}
+
+ ///
+ ~Constant() {}
+
+ ///
+ inline T apply(unsigned int /*i */ ) const { return rhs_; }
+
+protected:
+
+#ifdef SMATRIX_USE_REFERENCES
+ const T& rhs_;
+#else
+ const T rhs_;
+#endif
+
+};
+
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Expression */
diff --git a/smatrix/inc/Math/Functions.h b/smatrix/inc/Math/Functions.h
new file mode 100644
index 0000000000000000000000000000000000000000..ace5fb9de74d440c5eabe2443d3da2fa9bae1fe9
--- /dev/null
+++ b/smatrix/inc/Math/Functions.h
@@ -0,0 +1,392 @@
+// @(#)root/smatrix:$Name: $:$Id: Functions.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_Functions
+#define ROOT_Math_Functions
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 16. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Functions which are not applied like a unary operator
+//
+// changes:
+// 16 Mar 2001 (TG) creation
+// 03 Apr 2001 (TG) minimum added, doc++ comments added
+// 07 Apr 2001 (TG) Lmag2, Lmag added
+// 19 Apr 2001 (TG) added #include <cmath>
+// 24 Apr 2001 (TG) added sign()
+// 26 Jul 2001 (TG) round() added
+// 27 Sep 2001 (TG) added Expr declaration
+//
+// ********************************************************************
+#include <cmath>
+#include "Math/Expression.h"
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+template <class T, unsigned int D> class SVector;
+
+
+/** square.
+ Template to compute $x\cdot x$
+
+ @author T. Glebe
+*/
+//==============================================================================
+// square: x*x
+//==============================================================================
+template <class T>
+inline const T Square(const T& x) { return x*x; }
+
+/** maximum.
+ Template to compute $\max(i,j)$
+
+ @author T. Glebe
+*/
+//==============================================================================
+// maximum
+//==============================================================================
+template <class T>
+inline const T Maximum(const T& lhs, const T& rhs) {
+ return (lhs > rhs) ? lhs : rhs;
+}
+
+/** minimum.
+ Template to compute $\min(i,j)$
+
+ @author T. Glebe
+*/
+//==============================================================================
+// minimum
+//==============================================================================
+template <class T>
+inline const T Minimum(const T& lhs, const T& rhs) {
+ return (lhs < rhs) ? lhs : rhs;
+}
+
+/** round.
+ Template to compute nearest integer value.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// round
+//==============================================================================
+template <class T>
+inline int Round(const T& x) {
+ return (x-static_cast<int>(x) < 0.5) ? static_cast<int>(x) : static_cast<int>(x+1);
+}
+
+
+/** sign.
+ Template to compute the sign of a number $\textrm{sgn}(i)$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// sign
+//==============================================================================
+template <class T>
+inline const int Sign(const T& x) { return (x==0)? 0 : (x<0)? -1 : 1; }
+
+//==============================================================================
+// meta_dot
+//==============================================================================
+template <unsigned int I>
+struct meta_dot {
+ template <class A, class B, class T>
+ static inline T f(const A& lhs, const B& rhs, const T& x) {
+ return lhs.apply(I) * rhs.apply(I) + meta_dot<I-1>::f(lhs,rhs,x);
+ }
+};
+
+
+//==============================================================================
+// meta_dot<0>
+//==============================================================================
+template <>
+struct meta_dot<0> {
+ template <class A, class B, class T>
+ static inline T f(const A& lhs, const B& rhs, const T& /*x */) {
+ return lhs.apply(0) * rhs.apply(0);
+ }
+};
+
+
+/** dot.
+ Template to compute $\vec{a}\cdot\vec{b} = \sum_i a_i\cdot b_i$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// dot
+//==============================================================================
+template <class T, unsigned int D>
+inline T Dot(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
+ return meta_dot<D-1>::f(lhs,rhs, T());
+}
+
+//==============================================================================
+// dot
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Dot(const SVector<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ return meta_dot<D-1>::f(lhs,rhs, T());
+}
+
+//==============================================================================
+// dot
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Dot(const Expr<A,T,D>& lhs, const SVector<T,D>& rhs) {
+ return meta_dot<D-1>::f(lhs,rhs, T());
+}
+
+
+//==============================================================================
+// dot
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline T Dot(const Expr<A,T,D>& lhs, const Expr<B,T,D>& rhs) {
+ return meta_dot<D-1>::f(rhs,lhs, T());
+}
+
+
+//==============================================================================
+// meta_mag
+//==============================================================================
+template <unsigned int I>
+struct meta_mag {
+ template <class A, class T>
+ static inline T f(const A& rhs, const T& x) {
+ return Square(rhs.apply(I)) + meta_mag<I-1>::f(rhs, x);
+ }
+};
+
+
+//==============================================================================
+// meta_mag<0>
+//==============================================================================
+template <>
+struct meta_mag<0> {
+ template <class A, class T>
+ static inline T f(const A& rhs, const T& ) {
+ return Square(rhs.apply(0));
+ }
+};
+
+
+/** mag2.
+ Template to compute $|\vec{v}|^2 = \sum_iv_i^2$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// mag2
+//==============================================================================
+template <class T, unsigned int D>
+inline T Mag2(const SVector<T,D>& rhs) {
+ return meta_mag<D-1>::f(rhs, T());
+}
+
+//==============================================================================
+// mag2
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Mag2(const Expr<A,T,D>& rhs) {
+ return meta_mag<D-1>::f(rhs, T());
+}
+
+/** mag.
+ Length of a vector: $|\vec{v}| = \sqrt{\sum_iv_i^2}$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// mag
+//==============================================================================
+template <class T, unsigned int D>
+inline T Mag(const SVector<T,D>& rhs) {
+ return std::sqrt(Mag2(rhs));
+}
+
+//==============================================================================
+// mag
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Mag(const Expr<A,T,D>& rhs) {
+ return std::sqrt(Mag2(rhs));
+}
+
+
+/** Lmag2.
+ Template to compute $|\vec{v}|^2 = v_0^2 - v_1^2 - v_2^2 -v_3^2$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Lmag2
+//==============================================================================
+template <class T>
+inline T Lmag2(const SVector<T,4>& rhs) {
+ return Square(rhs[0]) - Square(rhs[1]) - Square(rhs[2]) - Square(rhs[3]);
+}
+
+//==============================================================================
+// Lmag2
+//==============================================================================
+template <class A, class T>
+inline T Lmag2(const Expr<A,T,4>& rhs) {
+ return Square(rhs.apply(0))
+ - Square(rhs.apply(1)) - Square(rhs.apply(2)) - Square(rhs.apply(3));
+}
+
+/** Lmag.
+ Length of a vector Lorentz-Vector: $|\vec{v}| = \sqrt{v_0^2 -
+ v_1^2 - v_2^2 -v_3^2}$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// Lmag
+//==============================================================================
+template <class T>
+inline T Lmag(const SVector<T,4>& rhs) {
+ return std::sqrt(Lmag2(rhs));
+}
+
+//==============================================================================
+// Lmag
+//==============================================================================
+template <class A, class T>
+inline T Lmag(const Expr<A,T,4>& rhs) {
+ return std::sqrt(Lmag2(rhs));
+}
+
+
+/** cross.
+ Cross product of two 3-dim vectors: $\vec{c} = \vec{a}\times\vec{b}$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// cross product
+//==============================================================================
+template <class T>
+inline SVector<T,3> Cross(const SVector<T,3>& lhs, const SVector<T,3>& rhs) {
+ return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
+ lhs.apply(2)*rhs.apply(1),
+ lhs.apply(2)*rhs.apply(0) -
+ lhs.apply(0)*rhs.apply(2),
+ lhs.apply(0)*rhs.apply(1) -
+ lhs.apply(1)*rhs.apply(0));
+}
+
+//==============================================================================
+// cross product
+//==============================================================================
+template <class A, class T>
+inline SVector<T,3> Cross(const Expr<A,T,3>& lhs, const SVector<T,3>& rhs) {
+ return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
+ lhs.apply(2)*rhs.apply(1),
+ lhs.apply(2)*rhs.apply(0) -
+ lhs.apply(0)*rhs.apply(2),
+ lhs.apply(0)*rhs.apply(1) -
+ lhs.apply(1)*rhs.apply(0));
+}
+
+//==============================================================================
+// cross product
+//==============================================================================
+template <class T, class A>
+inline SVector<T,3> Cross(const SVector<T,3>& lhs, const Expr<A,T,3>& rhs) {
+ return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
+ lhs.apply(2)*rhs.apply(1),
+ lhs.apply(2)*rhs.apply(0) -
+ lhs.apply(0)*rhs.apply(2),
+ lhs.apply(0)*rhs.apply(1) -
+ lhs.apply(1)*rhs.apply(0));
+}
+
+//==============================================================================
+// cross product
+//==============================================================================
+template <class A, class B, class T>
+inline SVector<T,3> Cross(const Expr<A,T,3>& lhs, const Expr<B,T,3>& rhs) {
+ return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
+ lhs.apply(2)*rhs.apply(1),
+ lhs.apply(2)*rhs.apply(0) -
+ lhs.apply(0)*rhs.apply(2),
+ lhs.apply(0)*rhs.apply(1) -
+ lhs.apply(1)*rhs.apply(0));
+}
+
+
+/** Unit.
+ Return a vector of unit lenght: $\vec{e}_v = \vec{v}/|\vec{v}|$.
+
+ @author T. Glebe
+*/
+//==============================================================================
+// unit: returns a unit vector
+//==============================================================================
+template <class T, unsigned int D>
+inline SVector<T,D> Unit(const SVector<T,D>& rhs) {
+ return SVector<T,D>(rhs).Unit();
+}
+
+//==============================================================================
+// unit: returns a unit vector
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline SVector<T,D> Unit(const Expr<A,T,D>& rhs) {
+ return SVector<T,D>(rhs).Unit();
+}
+
+#ifdef XXX
+//==============================================================================
+// unit: returns a unit vector (worse performance)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T>, T, D>
+ unit(const SVector<T,D>& lhs) {
+ typedef BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T> DivOpBinOp;
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
+}
+
+//==============================================================================
+// unit: returns a unit vector (worse performance)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<BinaryOp<DivOp<T>, Expr<A,T,D>, Constant<T>, T>, T, D>
+ unit(const Expr<A,T,D>& lhs) {
+ typedef BinaryOp<DivOp<T>, Expr<A,T,D>, Constant<T>, T> DivOpBinOp;
+ return Expr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
+}
+#endif
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_Functions */
diff --git a/smatrix/inc/Math/MConfig.h b/smatrix/inc/Math/MConfig.h
new file mode 100644
index 0000000000000000000000000000000000000000..1b2ef6a3365c1fa06f9e5a67a3b52cb49ca5b991
--- /dev/null
+++ b/smatrix/inc/Math/MConfig.h
@@ -0,0 +1,21 @@
+// @(#)root/smatrix:$Name: $:$Id: MConfig.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_MConfig_
+#define ROOT_Math_MConfig
+
+// for alpha streams
+#if defined(__alpha) && !defined(linux)
+# include <standards.h>
+# ifndef __USE_STD_IOSTREAM
+# define __USE_STD_IOSTREAM
+# endif
+#endif
+
+
+#if defined(__sun) && !defined(linux)
+#include <stdlib.h>
+#endif
+
+
+#endif
diff --git a/smatrix/inc/Math/MatrixFunctions.h b/smatrix/inc/Math/MatrixFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad2c93b9bcb6b6596a7f5c6abc0b23f5055e575b
--- /dev/null
+++ b/smatrix/inc/Math/MatrixFunctions.h
@@ -0,0 +1,528 @@
+// @(#)root/smatrix:$Name: $:$Id: MatrixFunctions.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_MatrixFunctions
+#define ROOT_Math_MatrixFunctions
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 20. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: Functions/Operators special to Matrix
+//
+// changes:
+// 20 Mar 2001 (TG) creation
+// 20 Mar 2001 (TG) Added Matrix * Vector multiplication
+// 21 Mar 2001 (TG) Added transpose, product
+// 11 Apr 2001 (TG) transpose() speed improvment by removing rows(), cols()
+// through static members of Matrix and Expr class
+//
+// ********************************************************************
+
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+#ifdef XXX
+//==============================================================================
+// SMatrix * SVector
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SVector<T,D1> operator*(const SMatrix<T,D1,D2>& rhs, const SVector<T,D2>& lhs)
+{
+ SVector<T,D1> tmp;
+ for(unsigned int i=0; i<D1; ++i) {
+ const unsigned int rpos = i*D2;
+ for(unsigned int j=0; j<D2; ++j) {
+ tmp[i] += rhs.apply(rpos+j) * lhs.apply(j);
+ }
+ }
+ return tmp;
+}
+#endif
+
+
+//==============================================================================
+// meta_row_dot
+//==============================================================================
+template <unsigned int I>
+struct meta_row_dot {
+ template <class A, class B>
+ static inline typename A::value_type f(const A& lhs, const B& rhs,
+ const unsigned int offset) {
+ return lhs.apply(offset+I) * rhs.apply(I) + meta_row_dot<I-1>::f(lhs,rhs,offset);
+ }
+};
+
+
+//==============================================================================
+// meta_row_dot<0>
+//==============================================================================
+template <>
+struct meta_row_dot<0> {
+ template <class A, class B>
+ static inline typename A::value_type f(const A& lhs, const B& rhs,
+ const unsigned int offset) {
+ return lhs.apply(offset) * rhs.apply(0);
+ }
+};
+
+//==============================================================================
+// VectorMatrixRowOp
+//==============================================================================
+template <class Matrix, class Vector, unsigned int D2>
+class VectorMatrixRowOp {
+public:
+ ///
+ VectorMatrixRowOp(const Matrix& lhs, const Vector& rhs) :
+ lhs_(lhs), rhs_(rhs) {}
+
+ ///
+ ~VectorMatrixRowOp() {}
+
+ /// calc $\sum_{j} a_{ij} * v_j$
+ inline typename Matrix::value_type apply(unsigned int i) const {
+ return meta_row_dot<D2-1>::f(lhs_, rhs_, i*D2);
+ }
+
+protected:
+ const Matrix& lhs_;
+ const Vector& rhs_;
+};
+
+
+//==============================================================================
+// meta_col_dot
+//==============================================================================
+template <unsigned int I>
+struct meta_col_dot {
+ template <class Matrix, class Vector>
+ static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
+ const unsigned int offset) {
+ return lhs.apply(Matrix::cols*I+offset) * rhs.apply(I) +
+ meta_col_dot<I-1>::f(lhs,rhs,offset);
+ }
+};
+
+
+//==============================================================================
+// meta_col_dot<0>
+//==============================================================================
+template <>
+struct meta_col_dot<0> {
+ template <class Matrix, class Vector>
+ static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
+ const unsigned int offset) {
+ return lhs.apply(offset) * rhs.apply(0);
+ }
+};
+
+//==============================================================================
+// VectorMatrixColOp
+//==============================================================================
+template <class Vector, class Matrix, unsigned int D1>
+class VectorMatrixColOp {
+public:
+ ///
+ VectorMatrixColOp(const Vector& lhs, const Matrix& rhs) :
+ lhs_(lhs), rhs_(rhs) {}
+
+ ///
+ ~VectorMatrixColOp() {}
+
+ /// calc $\sum_{j} a_{ij} * v_j$
+ inline typename Matrix::value_type apply(unsigned int i) const {
+ return meta_col_dot<D1-1>::f(rhs_, lhs_, i);
+ }
+
+protected:
+ const Vector& lhs_;
+ const Matrix& rhs_;
+};
+
+//==============================================================================
+// operator*: SMatrix * SVector
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixRowOp<SMatrix<T,D1,D2>,SVector<T,D2>, D2>, T, D1>
+ operator*(const SMatrix<T,D1,D2>& lhs, const SVector<T,D2>& rhs) {
+
+ typedef VectorMatrixRowOp<SMatrix<T,D1,D2>,SVector<T,D2>, D2> VMOp;
+ return Expr<VMOp, T, D1>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: SMatrix * Expr<A,T,D2>
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixRowOp<SMatrix<T,D1,D2>, Expr<A,T,D2>, D2>, T, D1>
+ operator*(const SMatrix<T,D1,D2>& lhs, const Expr<A,T,D2>& rhs) {
+
+ typedef VectorMatrixRowOp<SMatrix<T,D1,D2>,Expr<A,T,D2>, D2> VMOp;
+ return Expr<VMOp, T, D1>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: Expr<A,T,D1,D2> * SVector
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixRowOp<Expr<A,T,D1,D2>, SVector<T,D2>, D2>, T, D1>
+ operator*(const Expr<A,T,D1,D2>& lhs, const SVector<T,D2>& rhs) {
+
+ typedef VectorMatrixRowOp<Expr<A,T,D1,D2>,SVector<T,D2>, D2> VMOp;
+ return Expr<VMOp, T, D1>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: Expr<A,T,D1,D2> * Expr<B,T,D2>
+//==============================================================================
+template <class A, class B, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixRowOp<Expr<A,T,D1,D2>, Expr<B,T,D2>, D2>, T, D1>
+ operator*(const Expr<A,T,D1,D2>& lhs, const Expr<B,T,D2>& rhs) {
+
+ typedef VectorMatrixRowOp<Expr<A,T,D1,D2>,Expr<B,T,D2>, D2> VMOp;
+ return Expr<VMOp, T, D1>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: SVector * SMatrix
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2>, D1>, T, D2>
+ operator*(const SVector<T,D1>& lhs, const SMatrix<T,D1,D2>& rhs) {
+
+ typedef VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2>, D1> VMOp;
+ return Expr<VMOp, T, D2>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: SVector * Expr<A,T,D1,D2>
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2>, D1>, T, D2>
+ operator*(const SVector<T,D1>& lhs, const Expr<A,T,D1,D2>& rhs) {
+
+ typedef VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2>, D1> VMOp;
+ return Expr<VMOp, T, D2>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: Expr<A,T,D1> * SMatrix
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixColOp<Expr<A,T,D1>, SMatrix<T,D1,D2>, D1>, T, D2>
+ operator*(const Expr<A,T,D1>& lhs, const SMatrix<T,D1,D2>& rhs) {
+
+ typedef VectorMatrixColOp<Expr<A,T,D1>, SMatrix<T,D1,D2>, D1> VMOp;
+ return Expr<VMOp, T, D2>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator*: Expr<A,T,D1> * Expr<B,T,D1,D2>
+//==============================================================================
+template <class A, class B, class T, unsigned int D1, unsigned int D2>
+inline Expr<VectorMatrixColOp<Expr<A,T,D1>, Expr<B,T,D1,D2>, D1>, T, D2>
+ operator*(const Expr<A,T,D1>& lhs, const Expr<B,T,D1,D2>& rhs) {
+
+ typedef VectorMatrixColOp<Expr<A,T,D1>, Expr<B,T,D1,D2>, D1> VMOp;
+ return Expr<VMOp, T, D2>(VMOp(lhs,rhs));
+}
+
+//==============================================================================
+// meta_matrix_dot
+//==============================================================================
+template <unsigned int I>
+struct meta_matrix_dot {
+ template <class MatrixA, class MatrixB>
+ static inline typename MatrixA::value_type f(const MatrixA& lhs, const MatrixB& rhs,
+ const unsigned int offset) {
+ return lhs.apply(offset/MatrixB::cols*MatrixA::cols + I) *
+ rhs.apply(MatrixB::cols*I + offset%MatrixB::cols) +
+ meta_matrix_dot<I-1>::f(lhs,rhs,offset);
+ }
+};
+
+
+//==============================================================================
+// meta_matrix_dot<0>
+//==============================================================================
+template <>
+struct meta_matrix_dot<0> {
+ template <class MatrixA, class MatrixB>
+ static inline typename MatrixA::value_type f(const MatrixA& lhs, const MatrixB& rhs,
+ const unsigned int offset) {
+ return lhs.apply(offset/MatrixB::cols*MatrixA::cols) *
+ rhs.apply(offset%MatrixB::cols);
+ }
+};
+
+//==============================================================================
+// MatrixMulOp
+//==============================================================================
+template <class MatrixA, class MatrixB, class T, unsigned int D>
+class MatrixMulOp {
+public:
+ ///
+ MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
+ lhs_(lhs), rhs_(rhs) {}
+
+ ///
+ ~MatrixMulOp() {}
+
+ /// calc $\sum_{j} a_{ik} * b_{kj}$
+ inline T apply(unsigned int i) const {
+ return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
+ }
+
+protected:
+ const MatrixA& lhs_;
+ const MatrixB& rhs_;
+};
+
+
+//==============================================================================
+// operator* (SMatrix * SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<SMatrix<T,D1,D>, SMatrix<T,D,D2>,T,D>, T, D1, D2>
+ operator*(const SMatrix<T,D1,D>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef MatrixMulOp<SMatrix<T,D1,D>, SMatrix<T,D,D2>, T,D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (SMatrix * Expr, binary)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<SMatrix<T,D1,D>, Expr<A,T,D,D2>,T,D>, T, D1, D2>
+ operator*(const SMatrix<T,D1,D>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef MatrixMulOp<SMatrix<T,D1,D>, Expr<A,T,D,D2>,T,D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (Expr * SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<Expr<A,T,D1,D>, SMatrix<T,D,D2>,T,D>, T, D1, D2>
+ operator*(const Expr<A,T,D1,D>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef MatrixMulOp<Expr<A,T,D1,D>, SMatrix<T,D,D2>,T,D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (Expr * Expr, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<Expr<A,T,D1,D>, Expr<B,T,D,D2>,T,D>, T, D1, D2>
+ operator*(const Expr<A,T,D1,D>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef MatrixMulOp<Expr<A,T,D1,D>, Expr<B,T,D,D2>, T,D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+
+#ifdef XXX
+//==============================================================================
+// MatrixMulOp
+//==============================================================================
+template <class MatrixA, class MatrixB, unsigned int D>
+class MatrixMulOp {
+public:
+ ///
+ MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
+ lhs_(lhs), rhs_(rhs) {}
+
+ ///
+ ~MatrixMulOp() {}
+
+ /// calc $\sum_{j} a_{ik} * b_{kj}$
+ inline typename MatrixA::value_type apply(unsigned int i) const {
+ return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
+ }
+
+protected:
+ const MatrixA& lhs_;
+ const MatrixB& rhs_;
+};
+
+
+//==============================================================================
+// operator* (SMatrix * SMatrix, binary)
+//==============================================================================
+template < class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<SMatrix<T,D1,D>, SMatrix<T,D,D2>, D>, T, D1, D2>
+ operator*(const SMatrix<T,D1,D>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef MatrixMulOp<SMatrix<T,D1,D>, SMatrix<T,D,D2>, D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (SMatrix * Expr, binary)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<SMatrix<T,D1,D>, Expr<A,T,D,D2>, D>, T, D1, D2>
+ operator*(const SMatrix<T,D1,D>& lhs, const Expr<A,T,D,D2>& rhs) {
+ typedef MatrixMulOp<SMatrix<T,D1,D>, Expr<A,T,D,D2>, D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (Expr * SMatrix, binary)
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<Expr<A,T,D1,D>, SMatrix<T,D,D2>, D>, T, D1, D2>
+ operator*(const Expr<A,T,D1,D>& lhs, const SMatrix<T,D,D2>& rhs) {
+ typedef MatrixMulOp<Expr<A,T,D1,D>, SMatrix<T,D,D2>, D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+
+//==============================================================================
+// operator* (Expr * Expr, binary)
+//==============================================================================
+template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2>
+inline Expr<MatrixMulOp<Expr<A,T,D1,D>, Expr<B,T,D,D2>, D>, T, D1, D2>
+ operator*(const Expr<A,T,D1,D>& lhs, const Expr<B,T,D,D2>& rhs) {
+ typedef MatrixMulOp<Expr<A,T,D1,D>, Expr<B,T,D,D2>, D> MatMulOp;
+
+ return Expr<MatMulOp,T,D1,D2>(MatMulOp(lhs,rhs));
+}
+#endif
+
+//==============================================================================
+// TransposeOp
+//==============================================================================
+template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
+class TransposeOp {
+public:
+ ///
+ TransposeOp( const Matrix& rhs) :
+ rhs_(rhs) {}
+
+ ///
+ ~TransposeOp() {}
+
+ ///
+ inline T apply(unsigned int i) const {
+ return rhs_.apply( (i%D1)*D2 + i/D1);
+ }
+
+protected:
+ const Matrix& rhs_;
+};
+
+
+//==============================================================================
+// transpose
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline Expr<TransposeOp<SMatrix<T,D1,D2>,T,D1,D2>, T, D2, D1>
+ transpose(const SMatrix<T,D1,D2>& rhs) {
+ typedef TransposeOp<SMatrix<T,D1,D2>,T,D1,D2> MatTrOp;
+
+ return Expr<MatTrOp, T, D2, D1>(MatTrOp(rhs));
+}
+
+//==============================================================================
+// transpose
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2>
+inline Expr<TransposeOp<Expr<A,T,D1,D2>,T,D1,D2>, T, D2, D1>
+ transpose(const Expr<A,T,D1,D2>& rhs) {
+ typedef TransposeOp<Expr<A,T,D1,D2>,T,D1,D2> MatTrOp;
+
+ return Expr<MatTrOp, T, D2, D1>(MatTrOp(rhs));
+}
+
+//==============================================================================
+// product: SMatrix/SVector calculate v^T * A * v
+//==============================================================================
+template <class T, unsigned int D>
+inline T Product(const SMatrix<T,D>& lhs, const SVector<T,D>& rhs) {
+ return Dot(rhs, lhs * rhs);
+}
+
+//==============================================================================
+// product: SVector/SMatrix calculate v^T * A * v
+//==============================================================================
+template <class T, unsigned int D>
+inline T Product(const SVector<T,D>& lhs, const SMatrix<T,D>& rhs) {
+ return Dot(lhs, rhs * lhs);
+}
+
+//==============================================================================
+// product: SMatrix/Expr calculate v^T * A * v
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Product(const SMatrix<T,D>& lhs, const Expr<A,T,D>& rhs) {
+ return Dot(rhs, lhs * rhs);
+}
+
+//==============================================================================
+// product: Expr/SMatrix calculate v^T * A * v
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Product(const Expr<A,T,D>& lhs, const SMatrix<T,D>& rhs) {
+ return Dot(lhs, rhs * lhs);
+}
+
+//==============================================================================
+// product: SVector/Expr calculate v^T * A * v
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Product(const SVector<T,D>& lhs, const Expr<A,T,D,D>& rhs) {
+ return Dot(lhs, rhs * lhs);
+}
+
+//==============================================================================
+// product: Expr/SVector calculate v^T * A * v
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline T Product(const Expr<A,T,D,D>& lhs, const SVector<T,D>& rhs) {
+ return Dot(rhs, lhs * rhs);
+}
+
+//==============================================================================
+// product: Expr/Expr calculate v^T * A * v
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline T Product(const Expr<A,T,D,D>& lhs, const Expr<B,T,D>& rhs) {
+ return Dot(rhs, lhs * rhs);
+}
+
+//==============================================================================
+// product: Expr/Expr calculate v^T * A * v
+//==============================================================================
+template <class A, class B, class T, unsigned int D>
+inline T Product(const Expr<A,T,D>& lhs, const Expr<B,T,D,D>& rhs) {
+ return Dot(lhs, rhs * lhs);
+}
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+#endif /* ROOT_Math_MatrixFunctions */
diff --git a/smatrix/inc/Math/SMatrix.h b/smatrix/inc/Math/SMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..f7adaf532b451cabd35acbeed9edec57a571bafc
--- /dev/null
+++ b/smatrix/inc/Math/SMatrix.h
@@ -0,0 +1,281 @@
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_SMatrix
+#define ROOT_Math_SMatrix
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 20. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: A fixed size two dimensional Matrix class
+//
+// changes:
+// 20 Mar 2001 (TG) creation
+// 21 Mar 2001 (TG) added operators +=, -=, *=, /=
+// 26 Mar 2001 (TG) place_in_row(), place_in_col() added
+// 02 Apr 2001 (TG) non-const Array() added
+// 03 Apr 2001 (TG) invert() added
+// 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added
+// 09 Apr 2001 (TG) CTOR from array added
+// 11 Apr 2001 (TG) rows(), cols(), size() replaced by rows, cols, size
+// 25 Mai 2001 (TG) row(), col() added
+// 04 Sep 2001 (TG) moved inlined functions to .icc file
+// 11 Jan 2002 (TG) added operator==(), operator!=()
+// 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<()
+//
+// ********************************************************************
+// for platform specific configurations
+#include "Math/MConfig.h"
+
+#include <iosfwd>
+
+
+// expression engine
+
+#include "Math/Expression.h"
+
+
+namespace ROOT {
+
+ namespace Math {
+
+ template <class T, unsigned int D> class SVector;
+
+
+
+
+/** SMatrix.
+ A generic fixed size n x m Matrix class.q
+
+ @memo SMatrix
+ @author T. Glebe
+*/
+//==============================================================================
+// SMatrix: column-wise storage
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2 = D1>
+class SMatrix {
+public:
+ /** @name --- Typedefs --- */
+ ///
+ typedef T value_type;
+
+ /** @name --- Constructors --- */
+ ///
+ SMatrix();
+ ///
+ SMatrix(const SMatrix<T,D1,D2>& rhs);
+ ///
+ template <class A>
+ SMatrix(const Expr<A,T,D1,D2>& rhs);
+
+ // skip this methods (they are too ambigous)
+#ifdef OLD_IMPL
+ /// 2nd arg: set only diagonal?
+ SMatrix(const T& rhs, const bool diagonal=false);
+ /// constructor via dyadic product
+ SMatrix(const SVector<T,D1>& rhs);
+ /// constructor via dyadic product
+ template <class A>
+ SMatrix(const Expr<A,T,D1>& rhs);
+
+ /** constructor via array, triag=true: array contains only upper/lower
+ triangular part of a symmetric matrix, len: length of array */
+ template <class T1>
+ //SMatrix(const T1* a, const bool triang=false, const unsigned int len=D1*D2);
+ // to avoid clash with TRootIOconst
+ SMatrix(const T1* a, const bool triang, const unsigned int len=D1*D2);
+
+
+ /// assign from a scalar value
+ SMatrix<T,D1,D2>& operator=(const T& rhs);
+#endif
+
+ ///
+ template <class A>
+ SMatrix<T,D1,D2>& operator=(const Expr<A,T,D1,D2>& rhs);
+
+
+#ifdef OLD_IMPL
+ /// return no. of matrix rows
+ static const unsigned int rows = D1;
+ /// return no. of matrix columns
+ static const unsigned int cols = D2;
+ /// return no of elements: rows*columns
+ static const unsigned int size = D1*D2;
+#else
+ enum {
+ /// return no. of matrix rows
+ rows = D1,
+ /// return no. of matrix columns
+ cols = D2,
+ /// return no of elements: rows*columns
+ size = D1*D2
+ };
+#endif
+ /** @name --- Access functions --- */
+ /// access the parse tree
+ T apply(unsigned int i) const;
+ /// return read-only pointer to internal array
+ const T* Array() const;
+ /// return pointer to internal array
+ T* Array();
+
+ /** @name --- Operators --- */
+ /// element wise comparison
+ bool operator==(const T& rhs) const;
+ /// element wise comparison
+ bool operator!=(const T& rhs) const;
+ /// element wise comparison
+ bool operator==(const SMatrix<T,D1,D2>& rhs) const;
+ /// element wise comparison
+ bool operator!=(const SMatrix<T,D1,D2>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator==(const Expr<A,T,D1,D2>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator!=(const Expr<A,T,D1,D2>& rhs) const;
+
+ /// element wise comparison
+ bool operator>(const T& rhs) const;
+ /// element wise comparison
+ bool operator<(const T& rhs) const;
+ /// element wise comparison
+ bool operator>(const SMatrix<T,D1,D2>& rhs) const;
+ /// element wise comparison
+ bool operator<(const SMatrix<T,D1,D2>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator>(const Expr<A,T,D1,D2>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator<(const Expr<A,T,D1,D2>& rhs) const;
+
+ /// read-only access
+ const T& operator()(unsigned int i, unsigned int j) const;
+ /// read/write access
+ T& operator()(unsigned int i, unsigned int j);
+
+ ///
+ SMatrix<T,D1,D2>& operator+=(const SMatrix<T,D1,D2>& rhs);
+
+ SMatrix<T,D1,D2>& operator-=(const SMatrix<T,D1,D2>& rhs);
+
+ SMatrix<T,D1,D2>& operator*=(const SMatrix<T,D1,D2>& rhs);
+
+ SMatrix<T,D1,D2>& operator/=(const SMatrix<T,D1,D2>& rhs);
+
+
+#ifndef __CINT__
+ ///
+ template <class A>
+ SMatrix<T,D1,D2>& operator+=(const Expr<A,T,D1,D2>& rhs);
+ ///
+ ///
+ template <class A>
+ SMatrix<T,D1,D2>& operator-=(const Expr<A,T,D1,D2>& rhs);
+ ///
+ ///
+ template <class A>
+ SMatrix<T,D1,D2>& operator*=(const Expr<A,T,D1,D2>& rhs);
+ ///
+ ///
+ template <class A>
+ SMatrix<T,D1,D2>& operator/=(const Expr<A,T,D1,D2>& rhs);
+
+#endif
+
+ /** @name --- Expert functions --- */
+ /// invert symmetric, pos. def. Matrix via Dsinv
+ bool Sinvert();
+ /** determinant of symmetrc, pos. def. Matrix via Dsfact. \textbf{Note:} this
+ will destroy the contents of the Matrix!*/
+ bool Sdet(T& det);
+ /// invert square Matrix via Dinv
+ bool Invert();
+ /** determinant of square Matrix via Dfact. \textbf{Note:} this will destroy
+ the contents of the Matrix! */
+ bool Det(T& det);
+ /// place a vector in a Matrix row
+ template <unsigned int D>
+ SMatrix<T,D1,D2>& Place_in_row(const SVector<T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// place a vector expression in a Matrix row
+ template <class A, unsigned int D>
+ SMatrix<T,D1,D2>& Place_in_row(const Expr<A,T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// place a vector in a Matrix column
+ template <unsigned int D>
+ SMatrix<T,D1,D2>& Place_in_col(const SVector<T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// place a vector expression in a Matrix column
+ template <class A, unsigned int D>
+ SMatrix<T,D1,D2>& Place_in_col(const Expr<A,T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// place a matrix in this matrix
+ template <unsigned int D3, unsigned int D4>
+ SMatrix<T,D1,D2>& Place_at(const SMatrix<T,D3,D4>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// place a matrix expression in this matrix
+ template <class A, unsigned int D3, unsigned int D4>
+ SMatrix<T,D1,D2>& Place_at(const Expr<A,T,D3,D4>& rhs,
+ const unsigned int row,
+ const unsigned int col);
+ /// return a Matrix row as a vector
+ SVector<T,D2> Row(const unsigned int therow) const;
+ /// return a Matrix column as a vector
+ SVector<T,D1> Col(const unsigned int thecol) const;
+ /// used by operator<<()
+ std::ostream& Print(std::ostream& os) const;
+
+private:
+ T fArray[D1*D2];
+}; // end of class SMatrix
+
+
+
+//==============================================================================
+// operator<<
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline std::ostream& operator<<(std::ostream& os, const ROOT::Math::SMatrix<T,D1,D2>& rhs) {
+ return rhs.Print(os);
+}
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+
+
+
+#ifndef __CINT__
+
+#include "Math/SMatrix.icc"
+// include Matrix-Vector multiplication
+#include "Math/MatrixFunctions.h"
+
+#endif //__CINT__
+
+#endif /* ROOT_Math_SMatrix */
diff --git a/smatrix/inc/Math/SMatrix.icc b/smatrix/inc/Math/SMatrix.icc
new file mode 100644
index 0000000000000000000000000000000000000000..0c306ba0f8df1662db0ddb8c01099bdd16be4edf
--- /dev/null
+++ b/smatrix/inc/Math/SMatrix.icc
@@ -0,0 +1,597 @@
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.iccv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_SMatrix_icc
+#define ROOT_Math_SMatrix_icc
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 21. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: SMatrix implementation file
+//
+// changes:
+// 21 Mar 2001 (TG) creation
+// 26 Mar 2001 (TG) place_in_row(), place_in_col() added
+// 03 Apr 2001 (TG) invert() added
+// 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added
+// 09 Apr 2001 (TG) CTOR from array added
+// 25 Mai 2001 (TG) row(), col() added
+// 11 Jul 2001 (TG) added #include Functions.hh
+// 11 Jan 2002 (TG) added operator==(), operator!=()
+// 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<()
+//
+// ********************************************************************
+#include <iostream>
+#include <iomanip>
+#include <assert.h>
+#include "Math/Dsinv.h"
+#include "Math/Dsfact.h"
+#include "Math/Dfact.h"
+#include "Math/Dinv.h"
+#include "Math/Functions.h"
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+//==============================================================================
+// Constructors
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>::SMatrix() {
+ // operator=(0); // if operator=(T ) is defined
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] = 0;
+ }
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>::SMatrix(const SMatrix<T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] = rhs.fArray[i];
+ }
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>::SMatrix(const Expr<A,T,D1,D2>& rhs) {
+ operator=(rhs);
+}
+
+#ifdef OLD_IMPL
+
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>::SMatrix(const T& rhs, const bool diagonal) {
+ if(diagonal) {
+ // symmetric matrix!
+ assert(D1==D2);
+ // 50% faster than with nested loops
+ for(unsigned int i=0; i<D1*D2; ++i)
+ fArray[i] = 0;
+ for(unsigned int i=0; i<D1; ++i)
+ fArray[i*D1+i] = rhs;
+ } else {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] = rhs;
+ }
+ }
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>::SMatrix(const SVector<T,D1>& rhs) {
+ for(unsigned int i=0; i<D1; ++i) {
+ const unsigned int row = i*D1;
+ for(unsigned int j=0; j<D1; ++j) {
+ fArray[row+j] = rhs[i]*rhs[j];
+ }
+ }
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>::SMatrix(const Expr<A,T,D1>& rhs) {
+ for(unsigned int i=0; i<D1; ++i) {
+ const unsigned int row = i*D1;
+ for(unsigned int j=0; j<D1; ++j) {
+ fArray[row+j] = rhs.apply(i)*rhs.apply(j);
+ }
+ }
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class T1>
+SMatrix<T,D1,D2>::SMatrix(const T1* a, const bool triang, const unsigned int len) {
+ // fill from array with lower triangular matrix?
+ if(triang==false) {
+ for(unsigned int i=0; i<len; ++i) {
+ fArray[i] = a[i];
+ }
+ if(len<D1*D2)
+ for(unsigned int i=len; i<D1*D2; ++i)
+ fArray[i] = 0.;
+ } else {
+ // symmetric matrix required!
+ assert(D1==D2);
+
+ unsigned int k=0;
+ if(len == (Square(D1)+D1)/2) {
+ for(unsigned int i=0; i<D1; ++i) {
+ const unsigned int row = i*D1;
+ for(unsigned int j=0; j<=i && k<len; ++j) {
+ fArray[row+j] = fArray[j*D1+i] = a[k++];
+ }
+ }
+ } else {
+ unsigned int k=0;
+ for(unsigned int i=0; i<D1; ++i) {
+ const unsigned int row = i*D1;
+ for(unsigned int j=0; j<=i; ++j) {
+ if(k<len)
+ fArray[row+j] = fArray[j*D1+i] = a[k++];
+ else
+ fArray[row+j] = fArray[j*D1+i] = 0.;
+ } // for j
+ } // for i
+ } // if len ==
+ } // if triang
+}
+
+
+//==============================================================================
+// operator=
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator=(const T& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] = rhs;
+ }
+ return *this;
+}
+
+#endif
+
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator=(const Expr<A,T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] = rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// operator+=
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator+=(const SMatrix<T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] += rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator+=(const Expr<A,T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] += rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// operator-=
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator-=(const SMatrix<T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] -= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator-=(const Expr<A,T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] -= rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// operator*=
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator*=(const SMatrix<T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] *= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator*=(const Expr<A,T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] *= rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// operator/=
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator/=(const SMatrix<T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] /= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator/=(const Expr<A,T,D1,D2>& rhs) {
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ fArray[i] /= rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// operator== (element wise comparison)
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator==(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] == rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator==(const SMatrix<T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] == rhs.fArray[i]);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+bool SMatrix<T,D1,D2>::operator==(const Expr<A,T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] == rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// operator!= (element wise comparison)
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::operator!=(const T& rhs) const {
+ return !operator==(rhs);
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::operator!=(const SMatrix<T,D1,D2>& rhs) const {
+ return !operator==(rhs);
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+inline bool SMatrix<T,D1,D2>::operator!=(const Expr<A,T,D1,D2>& rhs) const {
+ return !operator==(rhs);
+}
+
+
+//==============================================================================
+// operator> (element wise comparison)
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator>(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] > rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator>(const SMatrix<T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] > rhs.fArray[i]);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+bool SMatrix<T,D1,D2>::operator>(const Expr<A,T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] > rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// operator< (element wise comparison)
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator<(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] < rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+bool SMatrix<T,D1,D2>::operator<(const SMatrix<T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] < rhs.fArray[i]);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <class A>
+bool SMatrix<T,D1,D2>::operator<(const Expr<A,T,D1,D2>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D1*D2; ++i) {
+ rc = rc && (fArray[i] < rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// sinvert
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::Sinvert() {
+ return Dsinv<T,D1,D1>(fArray);
+}
+
+//==============================================================================
+// sdet
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::Sdet(T& det) {
+ return Dsfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
+}
+
+//==============================================================================
+// invert
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::Invert() {
+ return Inverter<D2>::Dinv(*this);
+}
+
+//==============================================================================
+// det
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline bool SMatrix<T,D1,D2>::Det(T& det) {
+ return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
+}
+
+//==============================================================================
+// place_in_row
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int D>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_row(const SVector<T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+
+ assert(col+D <= D2);
+
+ for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
+ fArray[i] = rhs.apply(j);
+ }
+ return *this;
+}
+
+//==============================================================================
+// place_in_row
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <class A, unsigned int D>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_row(const Expr<A,T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+
+ assert(col+D <= D2);
+
+ for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
+ fArray[i] = rhs.apply(j);
+ }
+ return *this;
+}
+
+//==============================================================================
+// place_in_col
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int D>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_col(const SVector<T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+
+ assert(row+D <= D1);
+
+ for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
+ fArray[i] = rhs.apply(j);
+ }
+ return *this;
+}
+
+//==============================================================================
+// place_in_col
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <class A, unsigned int D>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_col(const Expr<A,T,D>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+
+ assert(row+D <= D1);
+
+ for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
+ fArray[i] = rhs.apply(j);
+ }
+ return *this;
+}
+
+//==============================================================================
+// place_at
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int D3, unsigned int D4>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_at(const SMatrix<T,D3,D4>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+ assert(row+D3 <= D1 && col+D4 <= D2);
+
+ const unsigned int offset = row*D2+col;
+
+ for(unsigned int i=0; i<D3*D4; ++i) {
+ fArray[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
+ }
+
+ return *this;
+}
+
+//==============================================================================
+// place_at
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <class A, unsigned int D3, unsigned int D4>
+SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_at(const Expr<A,T,D3,D4>& rhs,
+ const unsigned int row,
+ const unsigned int col) {
+ assert(row+D3 <= D1 && col+D4 <= D2);
+
+ const unsigned int offset = row*D2+col;
+
+ for(unsigned int i=0; i<D3*D4; ++i) {
+ fArray[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
+ }
+
+ return *this;
+}
+
+//==============================================================================
+// row
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SVector<T,D2> SMatrix<T,D1,D2>::Row(const unsigned int therow) const {
+
+ const unsigned int offset = therow*D2;
+
+ static SVector<T,D2> tmp;
+ for(unsigned int i=0; i<D2; ++i) {
+ tmp[i] = fArray[offset+i];
+ }
+ return tmp;
+}
+
+//==============================================================================
+// col
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+SVector<T,D1> SMatrix<T,D1,D2>::Col(const unsigned int thecol) const {
+
+ static SVector<T,D1> tmp;
+ for(unsigned int i=0; i<D1; ++i) {
+ tmp[i] = fArray[thecol+i*D2];
+ }
+ return tmp;
+}
+
+//==============================================================================
+// print
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+std::ostream& SMatrix<T,D1,D2>::Print(std::ostream& os) const {
+ os.setf(std::ios::right,std::ios::adjustfield);
+ // os.setf(ios::fixed);
+
+ os << "[ ";
+ for (unsigned int i=0; i < D1; ++i) {
+ for (unsigned int j=0; j < D2; ++j) {
+ os << std::setw(12) << fArray[i*D2+j];
+ if ((!((j+1)%12)) && (j < D2-1))
+ os << std::endl << " ...";
+ }
+ if (i != D1 - 1)
+ os << std::endl << " ";
+ }
+ os << " ]";
+
+ return os;
+}
+
+//==============================================================================
+// Access functions
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline T SMatrix<T,D1,D2>::apply(unsigned int i) const { return fArray[i]; }
+
+template <class T, unsigned int D1, unsigned int D2>
+inline const T* SMatrix<T,D1,D2>::Array() const { return fArray; }
+
+template <class T, unsigned int D1, unsigned int D2>
+inline T* SMatrix<T,D1,D2>::Array() { return fArray; }
+
+//==============================================================================
+// Operators
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+inline const T& SMatrix<T,D1,D2>::operator()(unsigned int i, unsigned int j) const {
+ return fArray[i*D2+j];
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+inline T& SMatrix<T,D1,D2>::operator()(unsigned int i, unsigned int j) {
+ return fArray[i*D2+j];
+}
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif
diff --git a/smatrix/inc/Math/SVector.h b/smatrix/inc/Math/SVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..57ea4909d3485575f35aee110d7d7dbd71dd0e06
--- /dev/null
+++ b/smatrix/inc/Math/SVector.h
@@ -0,0 +1,242 @@
+// @(#)root/smatrix:$Name: $:$Id: SVector.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_SVector
+#define ROOT_Math_SVector
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 16. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: A fixed size Vector class
+//
+// changes:
+// 16 Mar 2001 (TG) creation
+// 21 Mar 2001 (TG) SVector::value_type added
+// 21 Mar 2001 (TG) added operators +=, -=, *=, /=
+// 26 Mar 2001 (TG) added place_at()
+// 03 Apr 2001 (TG) Array() added
+// 06 Apr 2001 (TG) CTORS added
+// 07 Apr 2001 (TG) CTORS added
+// 22 Aug 2001 (TG) CTOR(T*,len) added
+// 04 Sep 2001 (TG) moved inlined functions to .icc file
+// 14 Jan 2002 (TG) added operator==(), operator!=(), operator>(), operator<()
+//
+// ********************************************************************
+
+#include "Math/MConfig.h"
+
+#include <iosfwd>
+
+
+// expression engine
+#include "Math/Expression.h"
+
+
+
+namespace ROOT {
+
+ namespace Math {
+
+ template <class A, class T, unsigned int D, unsigned int D2> class Expr;
+
+
+/** SVector.
+ A generic fixed size Vector class.
+
+ @memo SVector
+ @author T. Glebe
+*/
+//==============================================================================
+// SVector
+//==============================================================================
+template <class T, unsigned int D>
+class SVector {
+public:
+ /** @name --- Typedefs --- */
+ ///
+ typedef T value_type;
+
+ /** @name --- Constructors --- */
+ ///
+ SVector();
+ ///
+ template <class A>
+ SVector(const Expr<A,T,D>& rhs);
+ ///
+ SVector(const SVector<T,D>& rhs);
+ /// $D1\le D$ required!
+ template <unsigned int D1>
+ SVector(const SVector<T,D1>& rhs);
+ /// $D1\le D-1$ required!
+ template <unsigned int D1>
+ SVector(const T& a1, const SVector<T,D1>& rhs);
+ /// fill from array, len must be equal to D!
+ SVector(const T* a, unsigned int len);
+ ///
+ SVector(const T& rhs);
+ ///
+ SVector(const T& a1, const T& a2);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8,
+ const T& a9);
+ ///
+ SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8,
+ const T& a9, const T& a10);
+
+ ///
+ SVector<T,D>& operator=(const T& rhs);
+ ///
+ template <class A>
+ SVector<T,D>& operator=(const Expr<A,T,D>& rhs);
+
+ /** @name --- Access functions --- */
+ /// return dimension $D$
+ inline static unsigned int Dim() { return D; }
+ /// access the parse tree
+ T apply(unsigned int i) const;
+ /// return read-only pointer to internal array
+ const T* Array() const;
+ /// return pointer to internal array
+ T* Array();
+
+ /** @name --- Operators --- */
+ /// element wise comparison
+ bool operator==(const T& rhs) const;
+ /// element wise comparison
+ bool operator!=(const T& rhs) const;
+ /// element wise comparison
+ bool operator==(const SVector<T,D>& rhs) const;
+ /// element wise comparison
+ bool operator!=(const SVector<T,D>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator==(const Expr<A,T,D>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator!=(const Expr<A,T,D>& rhs) const;
+
+ /// element wise comparison
+ bool operator>(const T& rhs) const;
+ /// element wise comparison
+ bool operator<(const T& rhs) const;
+ /// element wise comparison
+ bool operator>(const SVector<T,D>& rhs) const;
+ /// element wise comparison
+ bool operator<(const SVector<T,D>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator>(const Expr<A,T,D>& rhs) const;
+ /// element wise comparison
+ template <class A>
+ bool operator<(const Expr<A,T,D>& rhs) const;
+
+ /// read-only access
+ const T& operator[](unsigned int i) const;
+ /// read-only access
+ const T& operator()(unsigned int i) const;
+ /// read/write access
+ T& operator[](unsigned int i);
+ /// read/write access
+ T& operator()(unsigned int i);
+
+ ///
+ SVector<T,D>& operator+=(const SVector<T,D>& rhs);
+ ///
+ SVector<T,D>& operator-=(const SVector<T,D>& rhs);
+ ///
+ SVector<T,D>& operator*=(const SVector<T,D>& rhs);
+ ///
+ SVector<T,D>& operator/=(const SVector<T,D>& rhs);
+
+
+#ifndef __CINT__
+ ///
+ template <class A>
+ SVector<T,D>& operator+=(const Expr<A,T,D>& rhs);
+ ///
+ template <class A>
+ SVector<T,D>& operator-=(const Expr<A,T,D>& rhs);
+ ///
+ template <class A>
+ SVector<T,D>& operator*=(const Expr<A,T,D>& rhs);
+ ///
+ template <class A>
+ SVector<T,D>& operator/=(const Expr<A,T,D>& rhs);
+
+#endif
+
+ /** @name --- Expert functions --- */
+ /// transform vector into a vector of lenght 1
+ SVector<T,D>& Unit();
+ /// place a sub-vector starting at <row>
+ template <unsigned int D2>
+ SVector<T,D>& Place_at(const SVector<T,D2>& rhs, const unsigned int row);
+ /// place a sub-vector starting at <row>
+ template <class A, unsigned int D2>
+ SVector<T,D>& Place_at(const Expr<A,T,D2>& rhs, const unsigned int row);
+ /// used by operator<<()
+ std::ostream& Print(std::ostream& os) const;
+
+private:
+ T fArray[D];
+}; // end of class SVector
+
+
+//==============================================================================
+// operator<<
+//==============================================================================
+template <class T, unsigned int D>
+std::ostream& operator<<(std::ostream& os, const ROOT::Math::SVector<T,D>& rhs);
+
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#ifndef __CINT__
+
+// include implementation file
+#include "Math/SVector.icc"
+
+// include operators and functions
+#include "Math/UnaryOperators.h"
+#include "Math/BinaryOperators.h"
+#include "Math/Functions.h"
+
+#endif // __CINT__
+
+
+#endif /* ROOT_Math_SVector */
diff --git a/smatrix/inc/Math/SVector.icc b/smatrix/inc/Math/SVector.icc
new file mode 100644
index 0000000000000000000000000000000000000000..ed07b07af7eec70c5183af2fbb62cdc3c121526c
--- /dev/null
+++ b/smatrix/inc/Math/SVector.icc
@@ -0,0 +1,465 @@
+// @(#)root/smatrix:$Name: $:$Id: SVector.iccv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_SVector_icc
+#define ROOT_Math_SVector_icc
+// ********************************************************************
+//
+// source:
+//
+// type: source code
+//
+// created: 21. Mar 2001
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+// Description: A fixed size Vector class
+//
+// changes:
+// 21 Mar 2001 (TG) creation
+// 26 Mar 2001 (TG) added place_at()
+// 06 Apr 2001 (TG) CTORS added
+// 07 Apr 2001 (TG) CTORS added
+// 22 Aug 2001 (TG) CTOR(T*,len) added
+// 14 Jan 2002 (TG) added operator==(), operator!=(), operator>(), operator<()
+//
+// ********************************************************************
+#include <iostream>
+#include <assert.h>
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+//==============================================================================
+// Constructors
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>::SVector() {
+ operator=(0);
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>::SVector(const Expr<A,T,D>& rhs) {
+ operator=(rhs);
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const SVector<T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i)
+ fArray[i] = rhs.fArray[i];
+}
+
+template <class T, unsigned int D>
+template <unsigned int D1>
+SVector<T,D>::SVector(const SVector<T,D1>& rhs) {
+ for(unsigned int i=0; i<D1; ++i)
+ fArray[i] = rhs.fArray[i];
+}
+
+template <class T, unsigned int D>
+template <unsigned int D1>
+SVector<T,D>::SVector(const T& a1, const SVector<T,D1>& rhs) {
+ fArray[0] = a1;
+ for(unsigned int i=0; i<D1; ++i)
+ fArray[i+1] = rhs[i];
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T* a, unsigned int len) {
+ assert(len == D);
+ for(unsigned int i=0; i<D; ++i)
+ fArray[i] = a[i];
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& rhs) {
+ for(unsigned int i=0; i<D; ++i)
+ fArray[i] = rhs;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2) {
+ fArray[0] = a1; fArray[1] = a2;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5; fArray[5] = a6;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5; fArray[5] = a6; fArray[6] = a7;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8,
+ const T& a9) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
+ fArray[8] = a9;
+}
+
+template <class T, unsigned int D>
+SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
+ const T& a5, const T& a6, const T& a7, const T& a8,
+ const T& a9, const T& a10) {
+ fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
+ fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
+ fArray[8] = a9; fArray[9] = a10;
+}
+
+//==============================================================================
+// operator=
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::operator=(const T& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] = rhs;
+ }
+ return *this;
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>& SVector<T,D>::operator=(const Expr<A,T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] = rhs.apply(i);
+ }
+ return *this;
+}
+
+//==============================================================================
+// operator==
+//==============================================================================
+template <class T, unsigned int D>
+bool SVector<T,D>::operator==(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] == rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+bool SVector<T,D>::operator==(const SVector<T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] == rhs.apply(i));
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+template <class A>
+bool SVector<T,D>::operator==(const Expr<A,T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] == rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// operator!=
+//==============================================================================
+template <class T, unsigned int D>
+inline bool SVector<T,D>::operator!=(const T& rhs) const {
+ return !operator==(rhs);
+}
+
+template <class T, unsigned int D>
+inline bool SVector<T,D>::operator!=(const SVector<T,D>& rhs) const {
+ return !operator==(rhs);
+}
+
+template <class T, unsigned int D>
+template <class A>
+inline bool SVector<T,D>::operator!=(const Expr<A,T,D>& rhs) const {
+ return !operator==(rhs);
+}
+
+//==============================================================================
+// operator>
+//==============================================================================
+template <class T, unsigned int D>
+bool SVector<T,D>::operator>(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] > rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+bool SVector<T,D>::operator>(const SVector<T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] > rhs.apply(i));
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+template <class A>
+bool SVector<T,D>::operator>(const Expr<A,T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] > rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// operator<
+//==============================================================================
+template <class T, unsigned int D>
+bool SVector<T,D>::operator<(const T& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] < rhs);
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+bool SVector<T,D>::operator<(const SVector<T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] < rhs.apply(i));
+ }
+ return rc;
+}
+
+template <class T, unsigned int D>
+template <class A>
+bool SVector<T,D>::operator<(const Expr<A,T,D>& rhs) const {
+ bool rc = true;
+ for(unsigned int i=0; i<D; ++i) {
+ rc = rc && (fArray[i] < rhs.apply(i));
+ }
+ return rc;
+}
+
+//==============================================================================
+// operator+=
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::operator+=(const SVector<T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] += rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>& SVector<T,D>::operator+=(const Expr<A,T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] += rhs.apply(i);
+ }
+ return *this;
+}
+
+//==============================================================================
+// operator-=
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::operator-=(const SVector<T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] -= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>& SVector<T,D>::operator-=(const Expr<A,T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] -= rhs.apply(i);
+ }
+ return *this;
+}
+
+//==============================================================================
+// operator*=
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::operator*=(const SVector<T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] *= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>& SVector<T,D>::operator*=(const Expr<A,T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] *= rhs.apply(i);
+ }
+ return *this;
+}
+
+//==============================================================================
+// operator/=
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::operator/=(const SVector<T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] /= rhs.apply(i);
+ }
+ return *this;
+}
+
+template <class T, unsigned int D>
+template <class A>
+SVector<T,D>& SVector<T,D>::operator/=(const Expr<A,T,D>& rhs) {
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] /= rhs.apply(i);
+ }
+ return *this;
+}
+
+
+//==============================================================================
+// unit
+//==============================================================================
+template <class T, unsigned int D>
+SVector<T,D>& SVector<T,D>::Unit() {
+ const double len = Mag(*this);
+ for(unsigned int i=0; i<D; ++i) {
+ fArray[i] /= len;
+ }
+ return *this;
+}
+
+//==============================================================================
+// place_at
+//==============================================================================
+template <class T, unsigned int D>
+template <unsigned int D2>
+SVector<T,D>& SVector<T,D>::Place_at(const SVector<T,D2>& rhs, const unsigned int row) {
+
+ assert(row+D2 <= D);
+ // Sassert(end <= D);
+
+ for(unsigned int i=row, j=0; j<D2; ++i,++j)
+ fArray[i] = rhs.apply(j);
+
+ return *this;
+}
+
+
+//==============================================================================
+// place_at
+//==============================================================================
+template <class T, unsigned int D>
+template <class A, unsigned int D2>
+SVector<T,D>& SVector<T,D>::Place_at(const Expr<A,T,D2>& rhs, const unsigned int row) {
+
+ assert(row+D2 <= D);
+
+ for(unsigned int i=row, j=0; j<D2; ++i,++j)
+ fArray[i] = rhs.apply(j);
+
+ return *this;
+}
+
+//==============================================================================
+// print
+//==============================================================================
+template <class T, unsigned int D>
+std::ostream& SVector<T,D>::Print(std::ostream& os) const {
+ os.setf(std::ios::right,std::ios::adjustfield);
+ // os.setf(ios::fixed);
+
+ unsigned int i=0;
+ for(; i<D-1; ++i)
+ os << fArray[i] << ", ";
+ os << fArray[i];
+ return os;
+}
+
+//==============================================================================
+// Access functions
+//==============================================================================
+template <class T, unsigned int D>
+inline T SVector<T,D>::apply(unsigned int i) const { return fArray[i]; }
+
+template <class T, unsigned int D>
+inline const T* SVector<T,D>::Array() const { return fArray; }
+
+template <class T, unsigned int D>
+inline T* SVector<T,D>::Array() { return fArray; }
+
+//==============================================================================
+// Operators
+//==============================================================================
+template <class T, unsigned int D>
+inline const T& SVector<T,D>::operator[](unsigned int i) const { return fArray[i]; }
+
+template <class T, unsigned int D>
+inline const T& SVector<T,D>::operator()(unsigned int i) const { return fArray[i]; }
+
+template <class T, unsigned int D>
+inline T& SVector<T,D>::operator[](unsigned int i) { return fArray[i]; }
+
+template <class T, unsigned int D>
+inline T& SVector<T,D>::operator()(unsigned int i) { return fArray[i]; }
+
+//==============================================================================
+// operator<<
+//==============================================================================
+template <class T, unsigned int D>
+inline std::ostream& operator<<(std::ostream& os, const SVector<T,D>& rhs) {
+ return rhs.Print(os);
+}
+
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+#endif
diff --git a/smatrix/inc/Math/UnaryOperators.h b/smatrix/inc/Math/UnaryOperators.h
new file mode 100644
index 0000000000000000000000000000000000000000..afc0c9efee90bd5dd30558a36649e6c425bb8474
--- /dev/null
+++ b/smatrix/inc/Math/UnaryOperators.h
@@ -0,0 +1,271 @@
+// @(#)root/smatrix:$Name: $:$Id: UnaryOperators.hv 1.0 2005/11/24 12:00:00 moneta Exp $
+// Authors: T. Glebe, L. Moneta 2005
+
+#ifndef ROOT_Math_UnaryOperators
+#define ROOT_Math_UnaryOperators
+//======================================================
+//
+// ATTENTION: This file was automatically generated,
+// do not edit!
+//
+// author: Thorsten Glebe
+// HERA-B Collaboration
+// Max-Planck-Institut fuer Kernphysik
+// Saupfercheckweg 1
+// 69117 Heidelberg
+// Germany
+// E-mail: T.Glebe@mpi-hd.mpg.de
+//
+//======================================================
+
+#include <cmath>
+
+namespace ROOT {
+
+ namespace Math {
+
+
+
+template <class T, unsigned int D> class SVector;
+template <class T, unsigned int D1, unsigned int D2> class SMatrix;
+
+
+//==============================================================================
+// Minus
+//==============================================================================
+template <class T>
+class Minus {
+public:
+ static inline T apply(const T& rhs) {
+ return -(rhs);
+ }
+};
+
+//==============================================================================
+// operator- (Expr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<UnaryOp<Minus<T>, Expr<A,T,D>, T>, T, D>
+ operator-(const Expr<A,T,D>& rhs) {
+ typedef UnaryOp<Minus<T>, Expr<A,T,D>, T> MinusUnaryOp;
+
+ return Expr<MinusUnaryOp,T,D>(MinusUnaryOp(Minus<T>(),rhs));
+}
+
+
+//==============================================================================
+// operator- (SVector, unary)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<UnaryOp<Minus<T>, SVector<T,D>, T>, T, D>
+ operator-(const SVector<T,D>& rhs) {
+ typedef UnaryOp<Minus<T>, SVector<T,D>, T> MinusUnaryOp;
+
+ return Expr<MinusUnaryOp,T,D>(MinusUnaryOp(Minus<T>(),rhs));
+}
+
+//==============================================================================
+// operator- (MatrixExpr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Minus<T>, Expr<A,T,D,D2>, T>, T, D, D2>
+ operator-(const Expr<A,T,D,D2>& rhs) {
+ typedef UnaryOp<Minus<T>, Expr<A,T,D,D2>, T> MinusUnaryOp;
+
+ return Expr<MinusUnaryOp,T,D,D2>(MinusUnaryOp(Minus<T>(),rhs));
+}
+
+
+//==============================================================================
+// operator- (SMatrix, unary)
+//==============================================================================
+template <class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Minus<T>, SMatrix<T,D,D2>, T>, T, D, D2>
+ operator-(const SMatrix<T,D,D2>& rhs) {
+ typedef UnaryOp<Minus<T>, SMatrix<T,D,D2>, T> MinusUnaryOp;
+
+ return Expr<MinusUnaryOp,T,D,D2>(MinusUnaryOp(Minus<T>(),rhs));
+}
+
+
+//==============================================================================
+// Fabs
+//==============================================================================
+template <class T>
+class Fabs {
+public:
+ static inline T apply(const T& rhs) {
+ return std::fabs(rhs);
+ }
+};
+
+//==============================================================================
+// fabs (Expr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<UnaryOp<Fabs<T>, Expr<A,T,D>, T>, T, D>
+ fabs(const Expr<A,T,D>& rhs) {
+ typedef UnaryOp<Fabs<T>, Expr<A,T,D>, T> FabsUnaryOp;
+
+ return Expr<FabsUnaryOp,T,D>(FabsUnaryOp(Fabs<T>(),rhs));
+}
+
+
+//==============================================================================
+// fabs (SVector, unary)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<UnaryOp<Fabs<T>, SVector<T,D>, T>, T, D>
+ fabs(const SVector<T,D>& rhs) {
+ typedef UnaryOp<Fabs<T>, SVector<T,D>, T> FabsUnaryOp;
+
+ return Expr<FabsUnaryOp,T,D>(FabsUnaryOp(Fabs<T>(),rhs));
+}
+
+//==============================================================================
+// fabs (MatrixExpr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Fabs<T>, Expr<A,T,D,D2>, T>, T, D, D2>
+ fabs(const Expr<A,T,D,D2>& rhs) {
+ typedef UnaryOp<Fabs<T>, Expr<A,T,D,D2>, T> FabsUnaryOp;
+
+ return Expr<FabsUnaryOp,T,D,D2>(FabsUnaryOp(Fabs<T>(),rhs));
+}
+
+
+//==============================================================================
+// fabs (SMatrix, unary)
+//==============================================================================
+template <class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Fabs<T>, SMatrix<T,D,D2>, T>, T, D, D2>
+ fabs(const SMatrix<T,D,D2>& rhs) {
+ typedef UnaryOp<Fabs<T>, SMatrix<T,D,D2>, T> FabsUnaryOp;
+
+ return Expr<FabsUnaryOp,T,D,D2>(FabsUnaryOp(Fabs<T>(),rhs));
+}
+
+
+//==============================================================================
+// Sqr
+//==============================================================================
+template <class T>
+class Sqr {
+public:
+ static inline T apply(const T& rhs) {
+ return square(rhs);
+ }
+};
+
+//==============================================================================
+// sqr (Expr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<UnaryOp<Sqr<T>, Expr<A,T,D>, T>, T, D>
+ sqr(const Expr<A,T,D>& rhs) {
+ typedef UnaryOp<Sqr<T>, Expr<A,T,D>, T> SqrUnaryOp;
+
+ return Expr<SqrUnaryOp,T,D>(SqrUnaryOp(Sqr<T>(),rhs));
+}
+
+
+//==============================================================================
+// sqr (SVector, unary)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<UnaryOp<Sqr<T>, SVector<T,D>, T>, T, D>
+ sqr(const SVector<T,D>& rhs) {
+ typedef UnaryOp<Sqr<T>, SVector<T,D>, T> SqrUnaryOp;
+
+ return Expr<SqrUnaryOp,T,D>(SqrUnaryOp(Sqr<T>(),rhs));
+}
+
+//==============================================================================
+// sqr (MatrixExpr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Sqr<T>, Expr<A,T,D,D2>, T>, T, D, D2>
+ sqr(const Expr<A,T,D,D2>& rhs) {
+ typedef UnaryOp<Sqr<T>, Expr<A,T,D,D2>, T> SqrUnaryOp;
+
+ return Expr<SqrUnaryOp,T,D,D2>(SqrUnaryOp(Sqr<T>(),rhs));
+}
+
+
+//==============================================================================
+// sqr (SMatrix, unary)
+//==============================================================================
+template <class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Sqr<T>, SMatrix<T,D,D2>, T>, T, D, D2>
+ sqr(const SMatrix<T,D,D2>& rhs) {
+ typedef UnaryOp<Sqr<T>, SMatrix<T,D,D2>, T> SqrUnaryOp;
+
+ return Expr<SqrUnaryOp,T,D,D2>(SqrUnaryOp(Sqr<T>(),rhs));
+}
+
+
+//==============================================================================
+// Sqrt
+//==============================================================================
+template <class T>
+class Sqrt {
+public:
+ static inline T apply(const T& rhs) {
+ return std::sqrt(rhs);
+ }
+};
+
+//==============================================================================
+// sqrt (Expr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D>
+inline Expr<UnaryOp<Sqrt<T>, Expr<A,T,D>, T>, T, D>
+ sqrt(const Expr<A,T,D>& rhs) {
+ typedef UnaryOp<Sqrt<T>, Expr<A,T,D>, T> SqrtUnaryOp;
+
+ return Expr<SqrtUnaryOp,T,D>(SqrtUnaryOp(Sqrt<T>(),rhs));
+}
+
+
+//==============================================================================
+// sqrt (SVector, unary)
+//==============================================================================
+template <class T, unsigned int D>
+inline Expr<UnaryOp<Sqrt<T>, SVector<T,D>, T>, T, D>
+ sqrt(const SVector<T,D>& rhs) {
+ typedef UnaryOp<Sqrt<T>, SVector<T,D>, T> SqrtUnaryOp;
+
+ return Expr<SqrtUnaryOp,T,D>(SqrtUnaryOp(Sqrt<T>(),rhs));
+}
+
+//==============================================================================
+// sqrt (MatrixExpr, unary)
+//==============================================================================
+template <class A, class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Sqrt<T>, Expr<A,T,D,D2>, T>, T, D, D2>
+ sqrt(const Expr<A,T,D,D2>& rhs) {
+ typedef UnaryOp<Sqrt<T>, Expr<A,T,D,D2>, T> SqrtUnaryOp;
+
+ return Expr<SqrtUnaryOp,T,D,D2>(SqrtUnaryOp(Sqrt<T>(),rhs));
+}
+
+
+//==============================================================================
+// sqrt (SMatrix, unary)
+//==============================================================================
+template <class T, unsigned int D, unsigned int D2>
+inline Expr<UnaryOp<Sqrt<T>, SMatrix<T,D,D2>, T>, T, D, D2>
+ sqrt(const SMatrix<T,D,D2>& rhs) {
+ typedef UnaryOp<Sqrt<T>, SMatrix<T,D,D2>, T> SqrtUnaryOp;
+
+ return Expr<SqrtUnaryOp,T,D,D2>(SqrtUnaryOp(Sqrt<T>(),rhs));
+}
+
+
+ } // namespace Math
+
+} // namespace ROOT
+
+
+
+#endif /* ROOT_Math_UnaryOperators */
diff --git a/smatrix/src/Dict.h b/smatrix/src/Dict.h
new file mode 100644
index 0000000000000000000000000000000000000000..ffb4ff9a2da680c5b0943cdc89dda0fb7fe71afe
--- /dev/null
+++ b/smatrix/src/Dict.h
@@ -0,0 +1,5 @@
+#include "SMatrix.hh"
+#include "SVector.hh"
+
+SMatrix<double,2,2> m;
+SVector<double,2> v;
diff --git a/smatrix/src/Selection.xml b/smatrix/src/Selection.xml
new file mode 100644
index 0000000000000000000000000000000000000000..411f2a0c39ec64c58f6eed1fc8a5219286634486
--- /dev/null
+++ b/smatrix/src/Selection.xml
@@ -0,0 +1,4 @@
+<lcgdict>
+ <class pattern ="SMatrix<double,*"/>
+ <class pattern ="SVector<double,*"/>
+</lcgdict>
\ No newline at end of file
diff --git a/smatrix/test/TestTimer.h b/smatrix/test/TestTimer.h
new file mode 100644
index 0000000000000000000000000000000000000000..312cbad43d361cd0b1a028c7994b2ae267f6f8db
--- /dev/null
+++ b/smatrix/test/TestTimer.h
@@ -0,0 +1,52 @@
+// simple class to measure time
+
+#include "TStopwatch.h"
+
+
+
+namespace ROOT {
+
+ namespace Math{
+
+ namespace test {
+
+
+ void printTime(TStopwatch & time, std::string s) {
+ int pr = std::cout.precision(8);
+ std::cout << s << "\t" << " time = " << time.RealTime() << "\t(sec)\t"
+ // << time.CpuTime()
+ << std::endl;
+ std::cout.precision(pr);
+ }
+
+ class Timer {
+
+ public:
+
+ Timer(const std::string & s = "") : fName(s), fTime(0)
+ {
+ fWatch.Start();
+ }
+ Timer(double & t, const std::string & s = "") : fName(s), fTime(&t)
+ {
+ fWatch.Start();
+ }
+
+ ~Timer() {
+ fWatch.Stop();
+ printTime(fWatch,fName);
+ if (fTime) *fTime += fWatch.RealTime();
+ }
+
+
+ private:
+
+ std::string fName;
+ double * fTime;
+ TStopwatch fWatch;
+
+ };
+ }
+
+ }
+}
diff --git a/smatrix/test/matrix_op.h b/smatrix/test/matrix_op.h
new file mode 100644
index 0000000000000000000000000000000000000000..2548045592cfaf3b3333e3f18828a3b5027b9bf5
--- /dev/null
+++ b/smatrix/test/matrix_op.h
@@ -0,0 +1,103 @@
+#include "TestTimer.h"
+
+// define funcitons for matrix operations
+
+//#define DEBUG
+#ifndef NLOOP
+#define NLOOP 1000000
+#endif
+
+
+// vector assignment
+template<class V>
+void testVeq(const V & v, double & time, V & result) {
+ test::Timer t(time,"V=V ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = v;
+ }
+}
+
+// matrix assignmnent
+template<class M>
+void testMeq(const M & m, double & time, M & result) {
+ test::Timer t(time,"M=M ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = m;
+ }
+}
+
+
+
+// vector sum
+template<class V>
+void testVad(const V & v1, const V & v2, double & time, V & result) {
+ test::Timer t(time,"V+V ");;
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = v1 + v2;
+ }
+}
+
+// matrix sum
+template<class M>
+void testMad(const M & m1, const M & m2, double & time, M & result) {
+ test::Timer t(time,"M+M ");;
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = m1 + m2;
+ }
+}
+
+// vector * constant
+template<class V>
+void testVscale(const V & v1, double a, double & time, V & result) {
+ test::Timer t(time,"a*V ");;
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = a * v1; // v1 * a does not exist in ROOT
+ }
+}
+
+
+// matrix * constant
+template<class M>
+void testMscale(const M & m1, double a, double & time, M & result) {
+ test::Timer t(time,"a*M ");;
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = m1 * a;
+ }
+}
+
+
+// simple Matrix vector op
+template<class M, class V>
+void testMV(const M & mat, const V & v, double & time, V & result) {
+ test::Timer t(time,"M*V ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = mat * v;
+ }
+}
+
+// general matrix vector op
+template<class M, class V>
+void testGMV(const M & mat, const V & v1, const V & v2, double & time, V & result) {
+ test::Timer t(time,"M*V+");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = mat * v1 + v2;
+ }
+}
+
+// general matrix matrix op
+template<class A, class B, class C>
+void testMM(const A & a, const B & b, const C & c, double & time, C & result) {
+ test::Timer t(time,"M*M ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = a * b + c;
+ }
+}
diff --git a/smatrix/test/testIO.C b/smatrix/test/testIO.C
new file mode 100644
index 0000000000000000000000000000000000000000..9cda39e44977d3bab80dd63991b3b8c75c81b871
--- /dev/null
+++ b/smatrix/test/testIO.C
@@ -0,0 +1,36 @@
+{
+gSystem->Load("libSmatrix");
+using namespace ROOT::Math;
+SMatrix<double,2,2> m;
+//cout << m << endl;
+
+double * d = m.Array();
+d[1] = 1;
+d[2] = 2;
+d[3] = 3;
+//m.Print();
+m.print(cout);
+
+TFile f("testSmatrix.root","RECREATE");
+
+f.WriteObjectAny(&m,"SMatrix<double,2,2>","m");
+
+f.Close();
+
+cout << "\nReading File..\n\n";
+
+TFile f2("testSmatrix.root");
+
+SMatrix<double,2,2> * p_m2 = (SMatrix<double,2,2> *) f2.GetObjectUnchecked("m");
+SMatrix<double,2,2> m2 = *p_m2;
+m2.print(cout);
+cout << endl;
+
+cout << " Test read matrix = original matrix ";
+if ( m2 == m )
+ cout << " OK " << endl;
+else
+cout << " Failed " << endl;
+
+f2.Close();
+}
diff --git a/smatrix/test/testKalman.cxx b/smatrix/test/testKalman.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9c0bcee215a753c38a078d77595dc859676e329a
--- /dev/null
+++ b/smatrix/test/testKalman.cxx
@@ -0,0 +1,283 @@
+
+
+#include "Math/SVector.h"
+#include "Math/SMatrix.h"
+
+#include "TMatrixD.h"
+#include "TVectorD.h"
+
+#include "TRandom3.h"
+
+// #include "SealUtil/SealTimer.h"
+// #include "SealUtil/SealHRRTChrono.h"
+// #include "SealUtil/TimingReport.h"
+
+#include <iostream>
+
+#ifndef NDIM1
+#define NDIM1 2
+#endif
+#ifndef NDIM2
+#define NDIM2 5
+#endif
+
+#define NITER 1 // number of iterations
+
+#define NLOOP 1000000 // number of time the test is repeted
+
+using namespace ROOT::Math;
+
+#include "TestTimer.h"
+
+int test_smatrix_kalman() {
+
+ // need to write explicitly the dimensions
+
+
+ typedef SMatrix<double, NDIM1, NDIM1> MnMatrixNN;
+ typedef SMatrix<double, NDIM2, NDIM2> MnMatrixMM;
+ typedef SMatrix<double, NDIM1, NDIM2> MnMatrixNM;
+ typedef SMatrix<double, NDIM2 , NDIM1> MnMatrixMN;
+ typedef SMatrix<double, NDIM1 > MnSymMatrixNN;
+ typedef SMatrix<double, NDIM2 > MnSymMatrixMM;
+ typedef SVector<double, NDIM1> MnVectorN;
+ typedef SVector<double, NDIM2> MnVectorM;
+
+
+
+ int first = NDIM1; //Can change the size of the matrices
+ int second = NDIM2;
+
+
+ std::cout << "************************************************\n";
+ std::cout << " SMatrix kalman test " << first << " x " << second << std::endl;
+ std::cout << "************************************************\n";
+
+
+
+
+ int npass = NITER;
+ TRandom3 r(111);
+ double x2sum = 0;
+ for (int k = 0; k < npass; k++) {
+
+
+
+ MnMatrixNM H;
+ MnMatrixMN K0;
+ MnSymMatrixMM Cp;
+ MnSymMatrixNN V;
+ MnVectorN m;
+ MnVectorM xp;
+
+
+ {
+
+ // fill matrices with ranodm data
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < second; j++)
+ H(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < first; j++)
+ K0(i,j) = r.Rndm() + 1.;
+ // sym matrices
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j <=i; j++) {
+ Cp(i,j) = r.Rndm() + 1.;
+ if (j != i) Cp(j,i) = Cp(i,j);
+ }
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j <=i; j++) {
+ V(i,j) = r.Rndm() + 1.;
+ if (j != i) V(j,i) = V(i,j);
+ }
+ // vectors
+ for(int i = 0; i < first; i++)
+ m(i) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ xp(i) = r.Rndm() + 1.;
+
+ }
+
+
+ MnSymMatrixMM I;
+ for(int i = 0; i < second; i++)
+ I(i,i) = 1;
+
+#ifdef DEBUG
+ std::cout << "pass " << k << std::endl;
+ if (k == 0) {
+ std::cout << " K0 = " << K0 << std::endl;
+ std::cout << " H = " << K0 << std::endl;
+ std::cout << " Cp = " << Cp << std::endl;
+ std::cout << " V = " << V << std::endl;
+ std::cout << " m = " << m << std::endl;
+ std::cout << " xp = " << xp << std::endl;
+ }
+#endif
+
+
+ {
+ double x2 = 0;
+ test::Timer t("SMatrix Kalman ");
+
+ MnVectorM x;
+ MnMatrixMN tmp;
+ MnSymMatrixNN Rinv;
+ MnMatrixMN K;
+ MnSymMatrixMM C;
+ MnVectorN vtmp;
+
+ for (int l = 0; l < NLOOP; l++)
+ {
+
+
+
+
+ x = xp + K0 * (m- H * xp);
+ tmp = Cp * transpose(H);
+ Rinv = V + H * tmp;
+
+ bool test = Rinv.Invert();
+ if(!test) {
+ std::cout<<"inversion failed" <<std::endl;
+ std::cout << Rinv << std::endl;
+ }
+
+ K = tmp * Rinv;
+ C = ( I - K * H ) * Cp;
+ //x2 = product(Rinv,m-H*xp); // this does not compile on WIN32
+ vtmp = m-H*xp;
+ x2 = Dot(vtmp, Rinv*vtmp);
+
+ }
+ //std::cout << k << " chi2 = " << x2 << std::endl;
+ x2sum += x2;
+ }
+ }
+ //tr.dump();
+
+ std::cout << "x2sum = " << x2sum << std::endl;
+
+ return 0;
+}
+
+// ROOT test
+
+
+int test_tmatrix_kalman() {
+
+
+
+
+ typedef TMatrixD MnMatrix;
+ typedef TVectorD MnVector;
+
+// typedef boost::numeric::ublas::matrix<double> MnMatrix;
+ //typedef HepSymMatrix MnSymMatrixHep;
+
+
+ int first = NDIM1; //Can change the size of the matrices
+ int second = NDIM2;
+
+
+ std::cout << "************************************************\n";
+ std::cout << " TMatrix Kalman test " << first << " x " << second << std::endl;
+ std::cout << "************************************************\n";
+
+
+
+ int npass = NITER;
+ TRandom3 r(111);
+ double x2sum = 0;
+
+ for (int k = 0; k < npass; k++)
+ {
+ MnMatrix H(first,second);
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < second; j++)
+ H(i,j) = r.Rndm() + 1.;
+ MnMatrix K0(second,first);
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < first; j++)
+ K0(i,j) = r.Rndm() + 1.;
+ //MnSymMatrix Cp(second);
+ MnMatrix Cp(second,second);
+ //in ROOT a sym matrix is like a normal matrix
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j <=i; j++){
+ Cp(i,j) = r.Rndm() + 1.;
+ if (j != i) Cp(j,i) = Cp(i,j);
+ }
+ //MnSymMatrix V(first);
+ MnMatrix V(first,first);
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j <=i; j++) {
+ V(i,j) = r.Rndm() + 1.;
+ if (j != i) V(j,i) = V(i,j);
+ }
+ MnVector m(first);
+ for(int j = 0; j < first; j++)
+ m(j) = r.Rndm() + 1.;
+ MnVector xp(second);
+ for(int j = 0; j < second; j++)
+ xp(j) = r.Rndm() + 1.;
+ //MnSymMatrix I(second);//Identity matrix
+ MnMatrix I(second,second);//Identity matrix
+ for (int i = 0; i < second; i++)
+ for(int j = 0; j <second; j++) {
+ I(i,j) = 0.0;
+ if (i==j) I(i,i) = 1.0;
+ }
+
+// if (k==0) {
+// std::cout << " Cp " << std::endl;
+// Cp.Print();
+// }
+
+ {
+ double x2 = 0;
+ MnVector x(second);
+ MnMatrix tmp(second,first);
+ MnMatrix Rinv(first,first);
+ MnMatrix K(second,first);
+ MnMatrix C(second,second);
+ MnVector vtmp(first);
+
+ test::Timer t("TMatrix Kalman ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ //MnVector x = xp + K0*(m-H*xp);
+ x = xp + K0*(m-H*xp);
+ MnMatrix HT = H;
+ //MnMatrix tmp = Cp*HT.Transpose(HT);
+ tmp = Cp*HT.Transpose(HT);
+ double det;
+ //MnMatrix Rinv =(V + H *tmp).InvertFast(&det);
+ Rinv =(V + H *tmp).InvertFast(&det);
+ //MnMatrix K = tmp * Rinv;
+ //MnMatrix C = (I - K * H) * Cp;
+ //K = tmp * Rinv;
+ C = (I - tmp* Rinv * H) * Cp;
+ vtmp = m-H*xp;
+ x2 = vtmp * (Rinv * vtmp);
+ }
+ //std::cout << k << " chi2 " << x2 << std::endl;
+ x2sum += x2;
+ }
+ // }
+ }
+ //tr.dump();
+ //print sum of x2 to check that result is same as other tests
+ std::cout << "x2sum = " << x2sum << std::endl;
+
+ return 0;
+}
+
+
+
+int main() {
+
+ test_smatrix_kalman();
+ test_tmatrix_kalman();
+}
diff --git a/smatrix/test/testOperations.cxx b/smatrix/test/testOperations.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8ff381ee27c495ea01eb9ffd946c6806a726338a
--- /dev/null
+++ b/smatrix/test/testOperations.cxx
@@ -0,0 +1,366 @@
+
+
+#include "Math/SVector.h"
+#include "Math/SMatrix.h"
+
+#include "TMatrixD.h"
+#include "TVectorD.h"
+
+#include "TRandom3.h"
+
+// #include "SealUtil/SealTimer.h"
+// #include "SealUtil/SealHRRTChrono.h"
+// #include "SealUtil/TimingReport.h"
+
+#include <iostream>
+
+#ifndef NDIM1
+#define NDIM1 5
+#endif
+#ifndef NDIM2
+#define NDIM2 5
+#endif
+
+#define NITER 1 // number of iterations
+
+using namespace ROOT::Math;
+
+#include "matrix_op.h"
+
+
+
+template<class V>
+double testDot_S(const V & v1, const V & v2, double & time) {
+ test::Timer t(time,"dot ");
+ double result=0;
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = Dot(v1,v2);
+ }
+ return result;
+}
+
+template<class M, class V>
+double testInnerProd_S(const M & a, const V & v, double & time) {
+ test::Timer t(time,"prod");
+ double result=0;
+ for (int l = 0; l < NLOOP; l++)
+ {
+#ifndef WIN32
+ result = Product(v,a);
+#else
+ // cannot instantiate on Windows (don't know why? )
+ V tmp = a*v;
+ result = Dot(v,tmp);
+#endif
+ }
+ return result;
+}
+
+//inversion
+template<class M>
+void testInv_S(const M & a, double & time, M& result){
+ test::Timer t(time,"inv ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = a;
+ result.Invert();
+ }
+}
+
+// for root
+
+
+template<class V>
+double testDot_T(const V & v1, const V & v2, double & time) {
+ test::Timer t(time,"dot ");
+ double result=0;
+ int size = v1.GetNrows();
+ assert(size == v2.GetNrows() );
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = 0;
+ for (int i = 0; i < size; ++i)
+ result += v1(i)*v2(i);
+ }
+ return result;
+}
+
+template<class M, class V>
+double testInnerProd_T(const M & a, const V & v, double & time) {
+ test::Timer t(time,"prod");
+ double result=0;
+ int size = v.GetNrows();
+ for (int l = 0; l < NLOOP; l++)
+ {
+ V tmp = a * v;
+ //result = v * tmp;
+ result = 0;
+ for (int i = 0; i < size; ++i)
+ result += v(i)*tmp(i);
+ }
+ return result;
+}
+
+//inversion
+template<class M>
+void testInv_T(const M & a, double & time, M& result){
+ test::Timer t(time,"inv ");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = a;
+ result.InvertFast();
+ }
+}
+
+template<class M>
+void testInv_T2(const M & a, double & time, M& result){
+ test::Timer t(time,"inv2");
+ for (int l = 0; l < NLOOP; l++)
+ {
+ result = a;
+ result.InvertFast();
+ }
+}
+
+
+
+int test_smatrix_op() {
+
+ // need to write explicitly the dimensions
+
+
+ typedef SMatrix<double, NDIM1, NDIM1> MnMatrixNN;
+ typedef SMatrix<double, NDIM2, NDIM2> MnMatrixMM;
+ typedef SMatrix<double, NDIM1, NDIM2> MnMatrixNM;
+ typedef SMatrix<double, NDIM2 , NDIM1> MnMatrixMN;
+ typedef SVector<double, NDIM1> MnVectorN;
+ typedef SVector<double, NDIM2> MnVectorM;
+
+
+
+ int first = NDIM1; //Can change the size of the matrices
+ int second = NDIM2;
+
+
+ std::cout << "************************************************\n";
+ std::cout << " SMatrix operations test " << first << " x " << second << std::endl;
+ std::cout << "************************************************\n";
+
+
+// seal::TimingReport tr;
+// seal::TimingItem & init = tr.item<seal::SealHRRTChrono>("init");
+
+// seal::TimingItem & t_veq = tr.item<seal::SealHRRTChrono>("smatrix veq");
+// seal::TimingItem & t_meq = tr.item<seal::SealHRRTChrono>("smatrix meq");
+// seal::TimingItem & t_vad = tr.item<seal::SealHRRTChrono>("smatrix vad");
+// seal::TimingItem & t_mad = tr.item<seal::SealHRRTChrono>("smatrix mad");
+// seal::TimingItem & t_dot = tr.item<seal::SealHRRTChrono>("smatrix dot");
+// seal::TimingItem & t_mv = tr.item<seal::SealHRRTChrono>("smatrix MV ");
+// seal::TimingItem & t_gmv = tr.item<seal::SealHRRTChrono>("smatrix GMV");
+// seal::TimingItem & t_mm = tr.item<seal::SealHRRTChrono>("smatrix GMM");
+// seal::TimingItem & t_prd = tr.item<seal::SealHRRTChrono>("smatrix prd");
+// seal::TimingItem & t_inv = tr.item<seal::SealHRRTChrono>("smatrix inv");
+
+ double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc = 0;
+
+
+
+ double r1,r2;
+ int npass = NITER;
+ TRandom3 r(111);
+ for (int k = 0; k < npass; k++) {
+
+
+ MnMatrixNM A;
+ MnMatrixMN B;
+ MnMatrixNN C;
+ MnMatrixMM D;
+ MnVectorN v;
+ MnVectorM p;
+
+ {
+ //seal::SealTimer t(init,false);
+ // fill matrices with ranodm data
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < second; j++)
+ A(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < first; j++)
+ B(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < first; j++)
+ C(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < second; j++)
+ D(i,j) = r.Rndm() + 1.;
+
+// // sym matrices
+// for(int i = 0; i < second; i++)
+// for(int j = 0; j <=i; j++) {
+// C(i,j) = r.Rndm() + 1.;
+// if (j != i) Cp(j,i) = Cp(i,j);
+// }
+ // vectors
+ for(int i = 0; i < first; i++)
+ v(i) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ p(i) = r.Rndm() + 1.;
+
+ }
+
+
+// MnSymMatrixMM I;
+// for(int i = 0; i < second; i++)
+// I(i,i) = 1;
+
+#ifdef DEBUG
+ std::cout << "pass " << k << std::endl;
+ if (k == 0) {
+ std::cout << " A = " << A << std::endl;
+ std::cout << " B = " << B << std::endl;
+ std::cout << " C = " << C << std::endl;
+ std::cout << " D = " << D << std::endl;
+ std::cout << " v = " << v << std::endl;
+ std::cout << " p = " << p << std::endl;
+ }
+#endif
+
+ MnVectorN v1; testMV(A,v,t_mv,v1);
+ MnVectorN v2; testGMV(A,v,v1,t_gmv,v2);
+ MnMatrixNN C1; testMM(A,B,C,t_mm,C1);
+ MnMatrixNN C2; testInv_S(C1,t_inv,C2);
+ MnVectorN v3; testVeq(v,t_veq,v3);
+ MnVectorN v4; testVad(v2,v3,t_vad,v4);
+ MnVectorN v5; testVscale(v4,2.0,t_vsc,v5);
+ MnMatrixNN C3; testMeq(C,t_meq,C3);
+ MnMatrixNN C4; testMad(C2,C3,t_mad,C4);
+ MnMatrixNN C5; testMscale(C4,0.5,t_msc,C5);
+
+ r1 = testDot_S(v3,v5,t_dot);
+ r2 = testInnerProd_S(C5,v5,t_prd);
+
+ }
+ //tr.dump();
+
+ std::cout << " r1 = " << r1 << " r2 = " << r2 << std::endl;
+
+ return 0;
+}
+
+// ROOT test
+
+
+int test_tmatrix_op() {
+
+
+
+
+ typedef TMatrixD MnMatrix;
+ typedef TVectorD MnVector;
+
+// typedef boost::numeric::ublas::matrix<double> MnMatrix;
+ //typedef HepSymMatrix MnSymMatrixHep;
+
+
+ int first = NDIM1; //Can change the size of the matrices
+ int second = NDIM2;
+
+
+ std::cout << "************************************************\n";
+ std::cout << " TMatrix operations test " << first << " x " << second << std::endl;
+ std::cout << "************************************************\n";
+
+ double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_inv2, t_vsc, t_msc = 0;
+
+
+ double r1,r2;
+ int npass = NITER;
+ TRandom3 r(111);
+
+ for (int k = 0; k < npass; k++) {
+
+
+ MnMatrix A(NDIM1,NDIM2);
+ MnMatrix B(NDIM2,NDIM1);
+ MnMatrix C(NDIM1,NDIM1);
+ MnMatrix D(NDIM2,NDIM2);
+ MnVector v(NDIM1);
+ MnVector p(NDIM2);
+
+ {
+ // seal::SealTimer t(init,false);
+ // fill matrices with ranodm data
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < second; j++)
+ A(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < first; j++)
+ B(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < first; i++)
+ for(int j = 0; j < first; j++)
+ C(i,j) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ for(int j = 0; j < second; j++)
+ D(i,j) = r.Rndm() + 1.;
+
+// // sym matrices
+// for(int i = 0; i < second; i++)
+// for(int j = 0; j <=i; j++) {
+// C(i,j) = r.Rndm() + 1.;
+// if (j != i) Cp(j,i) = Cp(i,j);
+// }
+ // vectors
+ for(int i = 0; i < first; i++)
+ v(i) = r.Rndm() + 1.;
+ for(int i = 0; i < second; i++)
+ p(i) = r.Rndm() + 1.;
+
+ }
+
+
+// MnSymMatrixMM I;
+// for(int i = 0; i < second; i++)
+// I(i,i) = 1;
+
+#ifdef DEBUG
+ std::cout << "pass " << k << std::endl;
+ if (k == 0) {
+ A.Print(); B.Print(); C.Print(); D.Print(); v.Print(); p.Print();
+ }
+#endif
+
+
+ MnVector v1(NDIM1); testMV(A,v,t_mv,v1);
+ MnVector v2(NDIM1); testGMV(A,v,v1,t_gmv,v2);
+ MnMatrix C1(NDIM1,NDIM1); testMM(A,B,C,t_mm,C1);
+ MnMatrix C2(NDIM1,NDIM1); testInv_T(C1,t_inv,C2);
+ MnVector v3(NDIM1); testVeq(v,t_veq,v3);
+ MnVector v4(NDIM1); testVad(v2,v3,t_vad,v4);
+ MnVector v5(NDIM1); testVscale(v4,2.0,t_vsc,v5);
+ MnMatrix C3(NDIM1,NDIM1); testMeq(C,t_meq,C3);
+ MnMatrix C4(NDIM1,NDIM1); testMad(C2,C3,t_mad,C4);
+ MnMatrix C5(NDIM1,NDIM1); testMscale(C4,0.5,t_msc,C5);
+
+
+ r1 = testDot_T(v3,v5,t_dot);
+ r2 = testInnerProd_T(C5,v5,t_prd);
+
+ MnMatrix C2b(NDIM1,NDIM1); testInv_T2(C1,t_inv2,C2b);
+
+
+ }
+ // tr.dump();
+
+ std::cout << " r1 = " << r1 << " r2 = " << r2 << std::endl;
+
+ return 0;
+}
+
+
+
+int main() {
+
+ test_smatrix_op();
+ test_tmatrix_op();
+}
diff --git a/smatrix/test/testSMatrix.cxx b/smatrix/test/testSMatrix.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f67f987de222a37f33752063fb3d1d4147f51dc3
--- /dev/null
+++ b/smatrix/test/testSMatrix.cxx
@@ -0,0 +1,235 @@
+#include <cmath>
+#include "Math/SVector.h"
+#include "Math/SMatrix.h"
+
+#include <iostream>
+
+using namespace ROOT::Math;
+
+using std::cout;
+using std::endl;
+
+
+#define XXX
+
+int test1() {
+
+ SVector<float,3> x(4,5,6);
+ SVector<float,2> y(2,3);
+ cout << "x: " << x << endl;
+ cout << "y: " << y << endl;
+
+ SMatrix<float,4,3> A;
+ SMatrix<float,2,2> B;
+
+ A.Place_in_row(y, 1, 1);
+ A.Place_in_col(x + 2, 1, 0);
+ A.Place_at(B , 2, 1);
+ cout << "A: " << endl << A << endl;
+
+ return 0;
+
+}
+
+int test2() {
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(1,0) = 1;
+ A(0,1) = 3;
+ A(1,1) = A(2,2) = 2;
+ cout << "A: " << endl << A << endl;
+
+ SVector<double,3> x = A.Row(0);
+ cout << "x: " << x << endl;
+
+ SVector<double,3> y = A.Col(1);
+ cout << "y: " << y << endl;
+
+ return 0;
+#endif
+}
+
+int test3() {
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 2;
+
+ SMatrix<double,3> B = A; // save A in B
+ cout << "A: " << endl << A << endl;
+
+ double det = 0.;
+ A.Sdet(det);
+ cout << "Determinant: " << det << endl;
+ // WARNING: A has changed!!
+ cout << "A again: " << endl << A << endl;
+
+ A.Sinvert();
+ cout << "A^-1: " << endl << A << endl;
+
+ // check if this is really the inverse:
+ cout << "A^-1 * B: " << endl << A * B << endl;
+
+ return 0;
+#endif
+}
+
+int test4() {
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 2;
+ cout << " A: " << endl << A << endl;
+
+ SVector<double,3> x(1,2,3);
+ cout << "x: " << x << endl;
+
+ // we add 1 to each component of x and A
+ cout << " (x+1)^T * (A+1) * (x+1): " << Product(x+1,A+1) << endl;
+
+ return 0;
+#endif
+}
+
+int test5() {
+#ifdef XXX
+ SMatrix<float,4,3> A;
+ A(0,0) = A(0,1) = A(1,1) = A(2,2) = 4.;
+ A(2,3) = 1.;
+ cout << "A: " << endl << A << endl;
+ SVector<float,4> x(1,2,3,4);
+ cout << " x: " << x << endl;
+ SVector<float,3> a(1,2,3);
+ cout << " a: " << a << endl;
+ SVector<float,4> y = x + A * a;
+ // SVector<float,4> y = A * a;
+ cout << " y: " << y << endl;
+
+ SVector<float,3> b = (x+1) * (A+1);
+ cout << " b: " << b << endl;
+
+ return 0;
+#endif
+}
+
+int test6() {
+#ifdef XXX
+ SMatrix<float,4,2> A;
+ A(0,0) = A(0,1) = A(1,1) = A(2,0) = A(3,1) = 4.;
+ cout << "A: " << endl << A << endl;
+
+ SMatrix<float,2,3> S;
+ S(0,0) = S(0,1) = S(1,1) = S(0,2) = 1.;
+ cout << " S: " << endl << S << endl;
+
+ SMatrix<float,4,3> C = A * S;
+ cout << " C: " << endl << C << endl;
+
+ return 0;
+#endif
+}
+
+int test7() {
+#ifdef XXX
+ SVector<float,3> xv(4,4,4);
+ SVector<float,3> yv(5,5,5);
+ SVector<float,2> zv(64,64);
+ SMatrix<float,2,3> x;
+ x.Place_in_row(xv,0,0);
+ x.Place_in_row(xv,1,0);
+ SMatrix<float,2,3> y;
+ y.Place_in_row(yv,0,0);
+ y.Place_in_row(yv,1,0);
+ SMatrix<float,2,3> z;
+ z.Place_in_col(zv,0,0);
+ z.Place_in_col(zv,0,1);
+ z.Place_in_col(zv,0,2);
+
+ // element wise multiplication
+ cout << "x * y: " << endl << times(x, -y) << endl;
+
+ x += z - y;
+ cout << "x += z - y: " << endl << x << endl;
+
+ // element wise square root
+ cout << "sqrt(z): " << endl << sqrt(z) << endl;
+
+ // element wise multiplication with constant
+ cout << "2 * y: " << endl << 2 * y << endl;
+
+ // a more complex expression
+ cout << "fabs(-z + 3*x): " << endl << fabs(-z + 3*x) << endl;
+
+ return 0;
+#endif
+}
+
+int test8() {
+#ifdef XXX
+ SMatrix<float,2,3> A;
+ SVector<float,3> av1(5.,15.,5.);
+ SVector<float,3> av2(15.,5.,15.);
+ A.Place_in_row(av1,0,0);
+ A.Place_in_row(av2,1,0);
+
+ cout << "A: " << endl << A << endl;
+
+ SVector<float,3> x(1,2,3);
+ SVector<float,3> y(4,5,6);
+
+ cout << "dot(x,y): " << Dot(x,y) << endl;
+
+ cout << "mag(x): " << Mag(x) << endl;
+
+ cout << "cross(x,y): " << Cross(x,y) << endl;
+
+ cout << "unit(x): " << Unit(x) << endl;
+
+ SVector<float,3> z(4,16,64);
+ cout << "x + y: " << x+y << endl;
+
+ cout << "x * y: " << x * -y << endl;
+ x += z - y;
+ cout << "x += z - y: " << x << endl;
+
+ // element wise square root
+ cout << "sqrt(z): " << sqrt(z) << endl;
+
+ // element wise multiplication with constant
+ cout << "2 * y: " << 2 * y << endl;
+
+ // a more complex expression
+ cout << "fabs(-z + 3*x): " << fabs(-z + 3*x) << endl;
+
+ SVector<float,4> a;
+ SVector<float,2> b(1,2);
+ a.Place_at(b,2);
+ cout << "a: " << a << endl;
+#endif
+
+
+ return 0;
+}
+
+#define TEST(N) \
+ itest = N; \
+ if (test##N() == 0) std::cout << " Test " << itest << " OK " << std::endl; \
+ else {std::cout << " Test " << itest << " Failed " << std::endl; \
+ return -1;}
+
+
+
+int main(void) {
+
+ int itest;
+ TEST(1);
+ TEST(2);
+ TEST(3);
+ TEST(4);
+ TEST(5);
+ TEST(6);
+ TEST(7);
+ TEST(8);
+
+ return 0;
+}
diff --git a/smatrix/test/test_smatrix.cxx b/smatrix/test/test_smatrix.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8434894032c33fa0f15d3656133c4c44f93d8c1c
--- /dev/null
+++ b/smatrix/test/test_smatrix.cxx
@@ -0,0 +1,516 @@
+#include <iostream>
+#include <cmath>
+#include "Math/SVector.h"
+#include "Math/SMatrix.h"
+
+using namespace ROOT::Math;
+
+using std::cout;
+using std::endl;
+
+int main(void) {
+
+ SVector<float,3> x(4,5,6);
+ SVector<float,2> y(2,3);
+ cout << "x: " << x << endl;
+ cout << "y: " << y << endl;
+
+ SMatrix<float,4,3> A;
+ SMatrix<float,2,2> B;
+
+ A.place_in_row(y, 1, 1);
+ A.place_in_col(x + 2, 1, 0);
+ A.place_at(B , 2, 1);
+ cout << "A: " << endl << A << endl;
+
+ return 0;
+
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(1,0) = 1;
+ A(0,1) = 3;
+ A(1,1) = A(2,2) = 2;
+ cout << "A: " << endl << A << endl;
+
+ SVector<double,3> x = A.row(0);
+ cout << "x: " << x << endl;
+
+ SVector<double,3> y = A.col(1);
+ cout << "y: " << y << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 2;
+
+ SMatrix<double,3> B = A; // save A in B
+ cout << "A: " << endl << A << endl;
+
+ double det = 0.;
+ A.sdet(det);
+ cout << "Determinant: " << det << endl;
+ // WARNING: A has changed!!
+ cout << "A again: " << endl << A << endl;
+
+ A.sinvert();
+ cout << "A^-1: " << endl << A << endl;
+
+ // check if this is really the inverse:
+ cout << "A^-1 * B: " << endl << A * B << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 2;
+ cout << " A: " << endl << A << endl;
+
+ SVector<double,3> x(1,2,3);
+ cout << "x: " << x << endl;
+
+ // we add 1 to each component of x and A
+ cout << " (x+1)^T * (A+1) * (x+1): " << product(x+1,A+1) << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<float,4,3> A;
+ A(0,0) = A(0,1) = A(1,1) = A(2,2) = 4.;
+ A(2,3) = 1.;
+ cout << "A: " << endl << A << endl;
+ SVector<float,4> x(1,2,3,4);
+ cout << " x: " << x << endl;
+ SVector<float,3> a(1,2,3);
+ cout << " a: " << a << endl;
+ SVector<float,4> y = x + A * a;
+ // SVector<float,4> y = A * a;
+ cout << " y: " << y << endl;
+
+ SVector<float,3> b = (x+1) * (A+1);
+ cout << " b: " << b << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<float,4,2> A;
+ A(0,0) = A(0,1) = A(1,1) = A(2,0) = A(3,1) = 4.;
+ cout << "A: " << endl << A << endl;
+
+ SMatrix<float,2,3> S;
+ S(0,0) = S(0,1) = S(1,1) = S(0,2) = 1.;
+ cout << " S: " << endl << S << endl;
+
+ SMatrix<float,4,3> C = A * S;
+ cout << " C: " << endl << C << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<float,2,3> x = 4.;
+ SMatrix<float,2,3> y = 5.;
+ SMatrix<float,2,3> z = 64.;
+
+ // element wise multiplication
+ cout << "x * y: " << endl << times(x, -y) << endl;
+
+ x += z - y;
+ cout << "x += z - y: " << endl << x << endl;
+
+ // element wise square root
+ cout << "sqrt(z): " << endl << sqrt(z) << endl;
+
+ // element wise multiplication with constant
+ cout << "2 * y: " << endl << 2 * y << endl;
+
+ // a more complex expression
+ cout << "fabs(-z + 3*x): " << endl << fabs(-z + 3*x) << endl;
+
+ return 0;
+#endif
+
+#ifdef XXX
+ SMatrix<float,2,3> A = 15.;
+ A(0,0) = A(1,1) = A(0,2) = 5.;
+
+ cout << "A: " << endl << A << endl;
+
+ SVector<float,3> x(1,2,3);
+ SVector<float,3> y(4,5,6);
+
+ cout << "dot(x,y): " << dot(x,y) << endl;
+
+ cout << "mag(x): " << mag(x) << endl;
+
+ cout << "cross(x,y): " << cross(x,y) << endl;
+
+ cout << "unit(x): " << unit(x) << endl;
+
+ SVector<float,3> z(4,16,64);
+ cout << "x + y: " << x+y << endl;
+
+ cout << "x * y: " << x * -y << endl;
+ x += z - y;
+ cout << "x += z - y: " << x << endl;
+
+ // element wise square root
+ cout << "sqrt(z): " << sqrt(z) << endl;
+
+ // element wise multiplication with constant
+ cout << "2 * y: " << 2 * y << endl;
+
+ // a more complex expression
+ cout << "fabs(-z + 3*x): " << fabs(-z + 3*x) << endl;
+
+ SVector<float,4> a;
+ SVector<float,2> b(1,2);
+ a.place_at(b,2);
+ cout << "a: " << a << endl;
+#endif
+
+ for(unsigned int i=0; i<1000000; ++i) {
+#ifdef XXX
+ VtVector a(1.,2.,3.);
+ VtVector b(4.,5.,6.);
+ VtVector c(8.,9.,10.);
+
+ VtVector d = a*(-1) + b;
+ VtVector d = a + b + c;
+#endif
+#ifdef XXX
+ VtMatrix A(4,3);
+ A(0,0) = A(1,1) = A(2,2) = 4.;
+ A(2,3) = 1.;
+ VtVector a(1,2,3);
+ // cout << " a: " << a << endl;
+ // cout << " A: " << A << endl;
+ // VtVector x(1,2,3,4);
+ // cout << " x: " << x << endl;
+ // VtVector y = x + A * a;
+ VtVector y = A * a;
+ cout << " y: " << y << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ // ==========================================
+ VtMatrix A(3,3);
+ A(0,0) = A(1,1) = A(2,2) = 4.;
+ // A(2,3) = 1.;
+ // cout << " A: " << -A << endl;
+ VtMatrix B(3,3);
+ B(0,1) = B(1,0) = B(1,1) = B(0,2) = 1.;
+ // cout << " B: " << B << endl;
+ // cout << " B: " << B << endl;
+ VtMatrix C(3,3);
+ C(0,2) = C(1,2) = C(2,2) = 2.;
+ // cout << " C: " << C << endl;
+ // VtMatrix D = B + C + (-A);
+ VtMatrix D = A + B + C + A + B + C;
+ cout << " D: " << D << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ VtSymMatrix VA(3,2.);
+ VA(0,0) = VA(0,1) = VA(1,0) = 1;
+ // cout << " A: " << VA << endl;
+ // cout << " det A: " << VA.det() << endl;
+ VA.det();
+ VA.VtDsinv();
+ cout << " A^-1: " << VA << endl;
+#endif
+#ifdef XXX
+ //===================================
+ VtSymMatrix A(3,3.);
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ // cout << " A: " << VA << endl;
+ VtSymMatrix B(3,3);
+ B(0,1) = B(1,0) = B(1,1) = B(2,2) = 1.;
+ // cout << " B: " << B << endl;
+ VtSymMatrix C(3,3);
+ C(0,1) = C(1,0) = C(2,1) = C(0,2) = 2.;
+ cout << " C: " << C << endl;
+ VtVector x(1,2,3);
+ // cout << " x: " << x << endl;
+ VtVector y(4,5,6);
+ // cout << " y: " << y << endl;
+ (A+B+C).product(x+y);
+ cout << " (x+y)^T * (A+B+C) * (x+y): " << (A+B+C).product(x+y) << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ // cout << " c: " << c << endl;
+ SVector<float,3> a;
+ a[0] = 1.; a[1] = 2.; a[2] = 3.;
+ SVector<float,3> b;
+ b[0] = 4.; b[1] = 5.; b[2] = 6.;
+ SVector<float,3> c;
+ c[0] = 8.; c[1] = 9.; c[2] = 10.;
+
+ SVector<float,3> d = a + b + c;
+ cout << "d: " << d << endl;
+ exit(0);
+
+ cout << "d: " << d << endl;
+ d -= b + c;
+ cout << "d -= b + c: " << d << endl;
+ cout << "dot(a,b): " << dot(a,b) << endl;
+ cout << "dot(a+b,c+d): " << dot(a+b,c+d) << endl;
+ cout << "dot(a*b,c+d): " << dot(a*b,c+d) << endl;
+ cout << "dot(a*b+c,c+d): " << dot(a*b+c,c+d) << endl;
+
+ cout << "mag2(a) " << mag2(a) << endl;
+ cout << "mag(a) " << mag(a) << endl;
+ cout << "mag2(a+b+c)" << mag2(a+b+c) << endl;
+ cout << "mag(a+b+c) " << mag(a+b+c) << endl;
+ // mag2(a);
+ // mag2(a+b+c);
+
+ // SVector<float,3> d = a + static_cast<float>(3.);
+ // SVector<float,3> d = sqr(a + b + -20);
+ // SVector<float,3> d = (-a+b) * 3;
+ SVector<float,3> d = a * (b + c);
+ d = unit(a + b + c);
+ // d.unit();
+ // cout << "d unit: " << d << endl;
+ // cout << "d = -a + b: " << d << endl;
+ // cout << "a: " << a << " b: " << b << endl;
+ // cout << "mag2(3*a): " << mag2(3*a) << endl;
+ // cout << "cross(3*a,b+c): " << cross(3*a,b+c) << endl;
+ // cout << "dot(-a,c+d+3): " << dot(-a,c+d+3) << endl;
+ // cout << "mag2(-a+d): " << mag2(-a+d) << " mag(d): " << mag(-a+d) << endl;
+
+ SMatrix<float,3,3> A = 15.;
+ SMatrix<float,3,3> B = A;
+ A(0,0) = A(1,1) = A(2,0) = A(2,2) = 5.;
+ SMatrix<float,3,3> C = fabs(A + -B + 2);
+// cout << "A: " << endl << A << endl;
+ cout << "C: " << endl << C << endl;
+#endif
+
+#ifdef XXX
+ SMatrix<float,4,3> A;
+ A(0,0) = A(1,1) = A(2,2) = 4.;
+ A(2,3) = 1.;
+ // cout << "A: " << endl << A << endl;
+ SVector<float,4> x(1,2,3,4);
+ // cout << " x: " << x << endl;
+ SVector<float,3> a(1,2,3);
+ // cout << " a: " << a << endl;
+ SVector<float,4> y = x + A * a;
+ // SVector<float,4> y = A * a;
+ cout << " y: " << y << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ // =======================================
+ SMatrix<float,3,3> A;
+ A(0,0) = A(1,1) = A(2,2) = 4.;
+ // cout << "A: " << endl << A << endl;
+ // cout << "B: " << endl << B << endl;
+ SMatrix<float,3,3> B;
+ B(0,1) = B(1,0) = B(1,1) = B(0,2) = 1.;
+ // cout << "B: " << endl << B << endl;
+ // SMatrix<float,3,3> C = times(A,B); // component wise multiplication
+ SMatrix<float,3,3> C;
+ C(0,2) = C(1,2) = C(2,2) = 2.;
+ // cout << "C: " << endl << C << endl;
+ // SMatrix<float,3,3> D = -A + B + C;
+ SMatrix<float,3,3> D = A + B + C + A + B + C;
+ cout << "D: " << endl << D << endl;
+ exit(0);
+ SMatrix<float,3,2> E = 4.;
+ cout << "E: " << endl << E << endl;
+ // D = A + E;
+#endif
+
+#ifdef XXX
+ SVector<float,3> b(4,5,6);
+ cout << " x: " << x << endl;
+
+ // SVector<float,4> y = x + (A * 2) * (a + 1);
+ SVector<float,3> y = (x+1) * (A+1);
+ cout << " y: " << y << endl;
+
+ SMatrix<float,3> S;
+ S(0,0) = S(1,0) = S(2,0) = 1.;
+ cout << " S: " << endl << S << endl;
+ SMatrix<float,4,3> C = A * S;
+ cout << " C: " << endl << C << endl;
+#endif
+#ifdef XXX
+ SMatrix<double,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 2;
+ // cout << "A: " << endl << A << endl;
+ double det = 0.;
+ A.sdet(det);
+ cout << "Determinant: " << det << endl;
+ cout << "A again: " << endl << A << endl;
+ exit(0);
+ A.sinvert();
+ cout << "A^-1: " << endl << A << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SVector<double,3> x(1,2,3);
+ cout << "x: " << x << endl;
+ cout << " x^T * A * x: " << product(x+1,A+1) << endl;
+ // product(A,x);
+
+ SMatrix<double,3> B = 1;
+ cout << "B: " << endl << B << endl;
+ A /= B + 1;
+ cout << "A/=B: " << endl << A << endl;
+
+ SVector<double,3> y(4,5,6);
+ cout << "y: " << y << endl;
+ y /= x;
+ cout << "y/=x: " << y << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ //===================================
+ SMatrix<float,3> A;
+ A(0,0) = A(0,1) = A(1,0) = 1;
+ A(1,1) = A(2,2) = 3.;
+ // cout << " A: " << endl << VA << endl;
+ SMatrix<float,3> B;
+ B(0,1) = B(1,0) = B(1,1) = B(2,2) = 1.;
+ B(0,0) = 3;
+ // cout << " B: " << endl << B << endl;
+ SMatrix<float,3> C;
+ C(0,1) = C(1,0) = C(2,1) = C(0,2) = 2.;
+ C(0,0) = C(1,1) = C(2,2) = 3;
+ // cout << " C: " << endl << C << endl;
+ SVector<float,3> x(1,2,3);
+ // cout << " x: " << x << endl;
+ SVector<float,3> y(4,5,6);
+ // cout << " y: " << y << endl;
+ product(A+B+C,x+y);
+ cout << " (x+y)^T * (A+B+C) * (x+y): " << product(A+B+C,x+y) << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SVector<float,4> x;
+ SVector<float,2> y(1,2);
+ // cout << "y: " << y << endl;
+ x.place_at(y,2);
+ cout << "x: " << x << endl;
+ SMatrix<float,4,3> A;
+ SMatrix<float,2,2> B = 1;
+ // SVector<float,2> x(1,2);
+ // A.place_in_row(x+2,1,1);
+ // A.place_in_col(x,1,1);
+ A.place_at(B,0,1);
+ cout << "A: " << endl << A << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ VtVector x(4);
+ VtVector y(1,2);
+ // cout << "y: " << y << endl;
+ x.place_at(y,2);
+ cout << "x: " << x << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,4,3> A;
+ SMatrix<float,2,2> B = 1;
+ // SVector<float,2> x(1,2);
+ // A.place_in_row(x+2,1,1);
+ // A.place_in_col(x,1,1);
+ A.place_at(B,0,1);
+ cout << "A: " << endl << A << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ VtMatrix A(4,3);
+ VtMatrix B(2,2);
+ B(0,0) = B(1,1) = 2.;
+ A.place_at(B,0,1);
+ cout << "A: " << A << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,3> C;
+ C(0,1) = C(1,0) = C(2,1) = C(0,2) = 2.;
+ C(0,0) = C(1,1) = C(2,2) = 3;
+ // cout << " C: " << endl << C << endl;
+ float det = 0.;
+ // Dfact<SMatrix<float,3>,3,3>(C,det);
+ C.det(det);
+ cout << "Dfact(C): " << det << endl;
+ cout << "C after: " << endl << C << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,4> C;
+ C(0,1) = C(1,0) = C(2,1) = C(0,2) = C(3,1) = C(2,3) = 2.;
+ C(0,0) = C(1,1) = C(2,2) = C(3,3) = 3;
+// SMatrix<float,4> B(C);
+// cout << " C: " << endl << C << endl;
+ Dinv<SMatrix<float,4>,4,4>(C);
+ cout << "C after: " << endl << C << endl;
+ SMatrix<float,4> D = B * C;
+ // cout.setf(ios::fixed);
+ cout << "D = B * C: " << endl << D << endl;
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,2> C;
+ C(0,0) = C(1,1) = C(0,1) = 5.;
+ C(1,0) = 1;
+ SMatrix<float,2> B(C);
+ cout << " C: " << endl << C << endl;
+ if(Invert<0>::Dinv(C) ) {
+ cout << "C after: " << endl << C << endl;
+ SMatrix<float,2> D = C * B;
+ cout << " D: " << endl << D << endl;
+ } else { cerr << " inversion failed! " << endl; }
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,3> C;
+ C(0,1) = C(1,1) = C(0,2) = C(2,2) = 5.;
+ C(0,0) = 10;
+ SMatrix<float,3> B(C);
+ cout << " C: " << endl << C << endl;
+ if(Invert<3>::Dinv(C)) {
+ cout << "C after: " << endl << C << endl;
+ SMatrix<float,3> D = B * C;
+ cout << "D: " << endl << D << endl;
+ } else { cerr << " inversion failed! " << endl; }
+ exit(0);
+#endif
+#ifdef XXX
+ SMatrix<float,4> C;
+ C(0,1) = C(1,0) = C(2,1) = C(0,2) = 2;
+ C(3,1) = 2;
+ C(2,3) = 2.;
+ C(0,0) = C(1,1) = C(2,2) = C(3,3) = 3;
+ // Invert<4>::Dinv<SMatrix<float,4>, 4>(C);
+ // C.invert();
+ SMatrix<float,4> B(C);
+ cout << " C: " << endl << C << endl;
+ // if(Invert<4>::Dinv(C)) {
+ if(C.invert()) {
+ cout << "C after: " << endl << C << endl;
+ SMatrix<float,4> D = B * C;
+ cout.setf(ios::fixed);
+ cout << "D: " << endl << D << endl;
+ } else { cerr << " inversion failed! " << endl; }
+
+ cout << "C+B: " << endl << C+B << endl;
+ exit(0);
+#endif
+ }
+
+ return 0;
+}