add methods for retrieve slice and blocks of matrices. Capitaliza name of...

add methods for retrieve slice and blocks of matrices. Capitaliza name of transpose function. Add as well tests for slices


git-svn-id: http://root.cern.ch/svn/root/trunk@13525 27541ba8-7e3a-0410-8455-c3a389f83636

smatrix/inc/Math/MatrixFunctions.h

-// @(#)root/smatrix:$Name:  $:$Id: MatrixFunctions.h,v 1.1 2005/11/24 16:03:42 brun Exp $
+// @(#)root/smatrix:$Name:  $:$Id: MatrixFunctions.h,v 1.2 2005/12/05 16:33:47 moneta Exp $
 // Authors: T. Glebe, L. Moneta    2005  
 
 #ifndef ROOT_Math_MatrixFunctions
@@ -438,7 +438,7 @@ protected:
 //==============================================================================
 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) {
+ 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));
@@ -449,7 +449,7 @@ inline Expr<TransposeOp<SMatrix<T,D1,D2>,T,D1,D2>, T, D2, D1>
 //==============================================================================
 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) {
+ 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));

smatrix/inc/Math/SMatrix.h

-// @(#)root/smatrix:$Name:  $:$Id: SMatrix.h,v 1.1 2005/11/24 16:03:42 brun Exp $
+// @(#)root/smatrix:$Name:  $:$Id: SMatrix.h,v 1.2 2005/12/05 16:33:47 moneta Exp $
 // Authors: T. Glebe, L. Moneta    2005  
 
 #ifndef ROOT_Math_SMatrix
@@ -318,11 +318,73 @@ public:
   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
+
+  /** 
+      return a full Matrix row as a vector (copy the content in a new vector)
+  */ 
   SVector<T,D2> Row(const unsigned int therow) const;
-  /// return a Matrix column as a vector
+
+  /** 
+      return a full Matrix column as a vector (copy the content in a new vector)
+  */
   SVector<T,D1> Col(const unsigned int thecol) const;
-  /// used by operator<<()
+
+  /**
+     return a slice of therow as a vector starting at the colum value col0 until col0+N. 
+     Condition  col0+N <= D2
+   */ 
+  template <unsigned int N>  
+  SVector<T,N> SubRow(const unsigned int therow, const unsigned int col0 = 0 ) const;
+
+  /**
+     return a slice of the column as a vector starting at the row value row0 until row0+Dsub.
+     Condition  row0+N <= D1
+   */ 
+  template <unsigned int N>  
+  SVector<T,N> SubCol(const unsigned int thecol, const unsigned int row0 = 0) const;
+
+  /**
+     return a submatrix with the upper left corner at the values (row0, col0) and with sizes N1, N2
+     Condition  row0+N1 <= D1 && col0+N2 <=D2
+   */ 
+  template <unsigned int N1, unsigned int N2 >  
+  SMatrix<T,N1,N2> SubMatrix(const unsigned int row0, const unsigned int col0) const;
+
+  /**
+     return diagonal elements of a matrix as a Vector.  
+     It works only for squared matrices D1 == D2, otherwise it will produce a compile error
+   */ 
+  SVector<T,D1> Diagonal() const;
+
+  /**
+     return the upper Triangular block of the matrices (including the diagonal) as 
+     a vector of sizes N = D1 * (D1 + 1)/2. 
+     It works only for square matrices with D1==D2, otherwise it will produce a compile error
+   */ 
+#ifndef UNSUPPORTED
+  SVector<T, D1 * (D2 +1)/2> UpperBlock() const;
+#else
+  template<unsigned int N>
+  SVector<T,N> UpperBlock() const;
+#endif
+  /**
+     return the lower Triangular block of the matrices (including the diagonal) as 
+     a vector of sizes N = D1 * (D1 + 1)/2. 
+     It works only for square matrices with D1==D2, otherwise it will produce a compile error
+   */ 
+#ifndef UNSUPPORTED
+  SVector<T, D1 * (D2 +1)/2> LowerBlock() const;
+#else
+  template<unsigned int N>
+  SVector<T,N> LowerBlock() const;
+#endif
+
+
+
+
+  // submatrices 
+
+  /// Print: used by operator<<()
   std::ostream& Print(std::ostream& os) const;
 
 private:

smatrix/inc/Math/SMatrix.icc

-// @(#)root/smatrix:$Name:  $:$Id: SMatrix.icc,v 1.1 2005/11/24 16:03:42 brun Exp $
+// @(#)root/smatrix:$Name:  $:$Id: SMatrix.icc,v 1.2 2005/12/05 16:33:47 moneta Exp $
 // Authors: T. Glebe, L. Moneta    2005  
 
 #ifndef ROOT_Math_SMatrix_icc
@@ -575,7 +575,7 @@ SVector<T,D2> SMatrix<T,D1,D2>::Row(const unsigned int therow) const {
 
   const unsigned int offset = therow*D2;
 
-  static SVector<T,D2> tmp;
+  /*static*/  SVector<T,D2> tmp;
   for(unsigned int i=0; i<D2; ++i) {
     tmp[i] = fArray[offset+i];
   }
@@ -588,7 +588,7 @@ SVector<T,D2> SMatrix<T,D1,D2>::Row(const unsigned int therow) const {
 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;
+  /*static */ SVector<T,D1> tmp;
   for(unsigned int i=0; i<D1; ++i) {
     tmp[i] = fArray[thecol+i*D2];
   }
@@ -668,6 +668,136 @@ inline const T * SMatrix<T,D1,D2>::end() const {
 
 
 
+//==============================================================================
+// SubMatrices and slices of columns and rows
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int N>  
+SVector<T,N> SMatrix<T,D1,D2>::SubRow(const unsigned int therow, const unsigned int col0 ) const { 
+
+  const unsigned int offset = therow*D2 + col0;
+
+  STATIC_CHECK( N <= D2,SVector_dimension_too_small); 
+  assert(col0 + N <= D2);
+
+  SVector<T,N> tmp;
+  for(unsigned int i=0; i<N; ++i) {
+    tmp[i] = fArray[offset+i];
+  }
+  return tmp;
+}
+
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int N>  
+SVector<T,N> SMatrix<T,D1,D2>::SubCol(const unsigned int thecol, const unsigned int row0 ) const { 
+
+  const unsigned int offset = thecol + row0*D1;
+
+  STATIC_CHECK( N <= D1,SVector_dimension_too_small); 
+  assert(row0 + N <= D1);
+
+  SVector<T,N> tmp;
+  for(unsigned int i=0; i<N; ++i) {
+    tmp[i] = fArray[offset+i*D1];
+  }
+  return tmp;
+}
+
+// submatrix
+template <class T, unsigned int D1, unsigned int D2>
+template <unsigned int N1, unsigned int N2>  
+SMatrix<T,N1,N2> SMatrix<T,D1,D2>::SubMatrix(const unsigned int row0, const unsigned int col0) const { 
+
+  STATIC_CHECK( N1 <= D1,Smatrix_dimension1_too_small); 
+  STATIC_CHECK( N2 <= D2,Smatrix_dimension2_too_small); 
+  assert(row0 + N1 <= D1);
+  assert(col0 + N2 <= D2);
+
+  SMatrix<T,N1,N2> tmp;
+  for(unsigned int i=0; i<N1; ++i) {
+    for(unsigned int j=0; j<N2; ++j) {
+      tmp(i,j) = fArray[(i+row0)*D2+j+col0];
+    }
+  }
+  return tmp;
+}
+
+//diagonal
+template <class T, unsigned int D1, unsigned int D2>
+SVector<T,D1> SMatrix<T,D1,D2>::Diagonal( ) const { 
+
+  // only for squared matrices
+  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
+
+  SVector<T,D1> tmp;
+  for(unsigned int i=0; i<D1; ++i) {
+    tmp[i] = fArray[ i*D2 + i];
+  }
+  return tmp;
+}
+
+//upper block
+template <class T, unsigned int D1, unsigned int D2>
+#ifndef UNSUPPORTED
+SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2>::UpperBlock( ) const { 
+#else
+template <unsigned int N>  
+SVector<T,N> SMatrix<T,D1,D2>::UpperBlock( ) const { 
+#endif
+  // only for squared matrices
+  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
+
+#ifndef UNSUPPORTED
+  SVector<T, D1 * (D2 +1)/2 >  tmp;
+#else
+  // N must be equal D1 *(D1 +1)/2
+  STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
+  SVector<T,N> tmp;
+#endif
+
+  int k = 0;
+  for(unsigned int i=0; i<D1; ++i) {
+    for(unsigned int j=i; j<D2; ++j) {
+      tmp[k] = fArray[ i*D2 + j];
+      k++;
+    }
+  }
+  return tmp;
+}
+
+//lower block
+template <class T, unsigned int D1, unsigned int D2>
+#ifndef UNSUPPORTED
+SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2>::LowerBlock( ) const { 
+#else
+template <unsigned int N>  
+SVector<T,N> SMatrix<T,D1,D2>::LowerBlock( ) const { 
+#endif
+
+  // only for squared matrices
+  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
+
+#ifndef UNSUPPORTED
+  SVector<T, D1 * (D2 +1)/2 >  tmp;
+#else
+  // N must be equal D1 *(D1 +1)/2
+  STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
+  SVector<T,N> tmp;
+#endif
+
+  int k = 0;
+  for(unsigned int i=0; i<D1; ++i) {
+    for(unsigned int j=0; j<=i; ++j) {
+      tmp[k] = fArray[ i*D2 + j];
+      k++;
+    }
+  }
+  return tmp;
+}
+
+
+
+
   }  // namespace Math
 
 }  // namespace ROOT

smatrix/test/testKalman.cxx

@@ -135,7 +135,7 @@ int test_smatrix_kalman() {
 
 
 	  x = xp + K0 * (m- H * xp);
-	  tmp = Cp * transpose(H);
+	  tmp = Cp * Transpose(H);
 	  Rinv = V + H * tmp;
 
 	  bool test = Rinv.Invert();

smatrix/test/testSMatrix.cxx

@@ -15,6 +15,12 @@ using std::endl;
 
 #define XXX
 
+int compare( double a, double b) { 
+  if (a == b) return 0; 
+  std::cout << "\nFailure " << a << " diffent than " << b << std::endl;
+  return 1;
+}
+
 int test1() { 
 
   SVector<float,3> x(4,5,6);
@@ -96,8 +102,9 @@ int test3() {
   cout << "Determinant: " << det << endl;
   // WARNING: A has changed!!
   cout << "A again: " << endl << A << endl;
+  A = B; 
 
-  A.Sinvert();
+  A.Invert();
   cout << "A^-1: " << endl << A << endl;
 
   // check if this is really the inverse:
@@ -244,8 +251,64 @@ int test8() {
   return 0;
 }
 
+int test9() { 
+  // test non mutating inversions
+  SMatrix<double,3> A;
+  A(0,0) = A(0,1) = A(1,0) = 1;
+  A(1,1) = A(2,2) = 2;
 
 
+  double det = 0.;
+  A.Det2(det);
+  cout << "Determinant: " << det << endl;
+
+  SMatrix<double,3> Ainv = A.Inverse();
+  cout << "A^-1: " << endl << Ainv << endl;
+
+  // check if this is really the inverse:
+  cout << "A^-1 * A: " << endl << Ainv * A << endl;
+
+  return 0;
+}
+
+int test10() { 
+  // test slices
+  int iret = 0;
+  double d[9] = { 1,2,3,4,5,6,7,8,9};
+  SMatrix<double,3> A( d,d+9);
+
+  cout << "A: " << A << endl;
+
+  SVector<double,2> v23 = A.SubRow<2>( 0,1);    
+  SVector<double,2> v69 = A.SubCol<2>( 2,1);    
+
+  std::cout << " v23 =  " << v23 << " \tv69 = " << v69 << std::endl;
+  iret |= compare( Dot(v23,v69),(2*6+3*9) ); 
+  
+  SMatrix<double,2,2> subA1 = A.SubMatrix<2,2>( 1,0);
+  SMatrix<double,2,3> subA2 = A.SubMatrix<2,3>( 0,0);
+  std::cout << " subA1 =  " << subA1 << " \nsubA2 = " << subA2 << std::endl;
+  iret |= compare ( subA1(0,0), subA2(1,0)); 
+  iret |= compare ( subA1(0,1), subA2(1,1)); 
+
+
+
+  SVector<double,3> diag = A.Diagonal();
+  std::cout << " diagonal =  " << diag << std::endl; 
+  iret |= compare( Mag2(diag) , 1+5*5+9*9 ); 
+
+  SMatrix<double,3> B = Transpose(A);
+  std::cout << " B = " << B << std::endl;
+
+  SVector<double,6> vU = A.UpperBlock();
+  SVector<double,6> vL = B.LowerBlock();
+  std::cout << " vU =  " << vU << " \tvL = " << vL << std::endl;
+  // need to test mag since order can change
+  iret |= compare( Mag(vU), Mag(vL) ); 
+  
+ 
+  return iret;
+}
 
 #define TEST(N)                                                                 \
   itest = N;                                                                    \
@@ -266,6 +329,8 @@ int main(void) {
   TEST(6);
   TEST(7);
   TEST(8);
+  TEST(9);
+  TEST(10);
 
   return 0;
 }