Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
HelperOps.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: J. Palacios 2006
3 
4 #ifndef ROOT_Math_HelperOps
5 #define ROOT_Math_HelperOps 1
6 
7 // Include files
8 
9 /** @class HelperOps HelperOps.h Math/HelperOps.h
10  *
11  *
12  * @author Juan PALACIOS
13  * @date 2006-01-11
14  *
15  * Specialised helper classes for binary operators =, +=, -=
16  * between SMatrices and Expressions with arbitrary representations.
17  * Specialisations at the moment only for Symmetric LHS and Generic RHS
18  * and used to throw static assert.
19  */
20 #include "Math/StaticCheck.h"
21 #include <algorithm> // required by std::copy
22 #include <assert.h>
23 
24 namespace ROOT {
25 
26 namespace Math {
27 
28  template <class T, unsigned int D1, unsigned int D2, class R>
29  class SMatrix;
30 
31  template <class A, class T, unsigned int D1, unsigned int D2, class R>
32  class Expr;
33 
34  template <class T, unsigned int D>
35  class MatRepSym;
36 
37  template <class T, unsigned int D1, unsigned int D2>
38  class MatRepStd;
39 
40  //=========================================================================
41  /**
42  Structure to assign from an expression based to general matrix to general matrix
43  */
44  template <class T,
45  unsigned int D1, unsigned int D2,
46  class A, class R1, class R2>
47 
48  struct Assign
49  {
50  /**
51  Evaluate the expression from general to general matrices.
52  If the matrix to assign the value is in use in the expression,
53  a temporary object is created to store the value (case A = B * A)
54  */
55  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const Expr<A,T,D1,D2,R2>& rhs)
56  {
57  if (! rhs.IsInUse(lhs.begin() ) ) {
58  unsigned int l = 0;
59  for(unsigned int i=0; i<D1; ++i)
60  for(unsigned int j=0; j<D2; ++j) {
61  lhs.fRep[l] = rhs(i,j);
62  l++;
63  }
64  }
65  // lhs is in use in expression, need to create a temporary with the result
66  else {
67  // std::cout << "create temp for " << typeid(rhs).name() << std::endl;
68  T tmp[D1*D2];
69  unsigned int l = 0;
70  for(unsigned int i=0; i<D1; ++i)
71  for(unsigned int j=0; j<D2; ++j) {
72  tmp[l] = rhs(i,j);
73  l++;
74  }
75  // copy now the temp object
76  for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = tmp[i];
77  }
78 
79  }
80 
81  };
82 
83  /**
84  Structure to assign from an expression based to symmetric matrix to symmetric matrix
85  */
86  template <class T,
87  unsigned int D1, unsigned int D2,
88  class A>
89 
90  struct Assign<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
91  {
92  /**
93  Evaluate the expression from symmetric to symmetric matrices.
94  If the matrix to assign the value is in use in the expression,
95  a temporary object is created to store the value (case A = B * A)
96  */
97  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
98  const Expr<A,T,D1,D2,MatRepSym<T,D1> >& rhs)
99  {
100  if (! rhs.IsInUse(lhs.begin() ) ) {
101  unsigned int l = 0;
102  for(unsigned int i=0; i<D1; ++i)
103  // storage of symmetric matrix is in lower block
104  for(unsigned int j=0; j<=i; ++j) {
105  lhs.fRep.Array()[l] = rhs(i,j);
106  l++;
107  }
108  }
109  // create a temporary object to store result
110  else {
111  T tmp[MatRepSym<T,D1>::kSize];
112  unsigned int l = 0;
113  for(unsigned int i=0; i<D1; ++i)
114  for(unsigned int j=0; j<=i; ++j) {
115  tmp[l] = rhs(i,j);
116  l++;
117  }
118  // copy now the temp object
119  for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] = tmp[i];
120  }
121  }
122  };
123 
124 
125 
126  /**
127  Dummy Structure which flags an error to avoid assigment from expression based on a
128  general matrix to a symmetric matrix
129  */
130  template <class T, unsigned int D1, unsigned int D2, class A>
131  struct Assign<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
132  {
133  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
134  const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
135  {
136  STATIC_CHECK(0==1, Cannot_assign_general_to_symmetric_matrix);
137  }
138 
139  }; // struct Assign
140 
141 
142  /**
143  Force Expression evaluation from general to symmetric.
144  To be used when is known (like in similarity products) that the result
145  is symmetric
146  Note this is function used in the simmilarity product: no check for temporary is
147  done since in that case is not needed
148  */
149  struct AssignSym
150  {
151  /// assign a symmetric matrix from an expression
152  template <class T,
153  unsigned int D,
154  class A,
155  class R>
156  static void Evaluate(SMatrix<T,D,D,MatRepSym<T,D> >& lhs, const Expr<A,T,D,D,R>& rhs)
157  {
158  //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
159  unsigned int l = 0;
160  for(unsigned int i=0; i<D; ++i)
161  // storage of symmetric matrix is in lower block
162  for(unsigned int j=0; j<=i; ++j) {
163  lhs.fRep.Array()[l] = rhs(i,j);
164  l++;
165  }
166  }
167  /// assign the symmetric matric from a general matrix
168  template <class T,
169  unsigned int D,
170  class R>
171  static void Evaluate(SMatrix<T,D,D,MatRepSym<T,D> >& lhs, const SMatrix<T,D,D,R>& rhs)
172  {
173  //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
174  unsigned int l = 0;
175  for(unsigned int i=0; i<D; ++i)
176  // storage of symmetric matrix is in lower block
177  for(unsigned int j=0; j<=i; ++j) {
178  lhs.fRep.Array()[l] = rhs(i,j);
179  l++;
180  }
181  }
182 
183 
184  }; // struct AssignSym
185 
186 
187  //=========================================================================
188  /**
189  Evaluate the expression performing a += operation
190  Need to check whether creating a temporary object with the expression result
191  (like in op: A += A * B )
192  */
193  template <class T, unsigned int D1, unsigned int D2, class A,
194  class R1, class R2>
195  struct PlusEquals
196  {
197  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const Expr<A,T,D1,D2,R2>& rhs)
198  {
199  if (! rhs.IsInUse(lhs.begin() ) ) {
200  unsigned int l = 0;
201  for(unsigned int i=0; i<D1; ++i)
202  for(unsigned int j=0; j<D2; ++j) {
203  lhs.fRep[l] += rhs(i,j);
204  l++;
205  }
206  }
207  else {
208  T tmp[D1*D2];
209  unsigned int l = 0;
210  for(unsigned int i=0; i<D1; ++i)
211  for(unsigned int j=0; j<D2; ++j) {
212  tmp[l] = rhs(i,j);
213  l++;
214  }
215  // += now using the temp object
216  for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] += tmp[i];
217  }
218  }
219  };
220 
221  /**
222  Specialization for symmetric matrices
223  Evaluate the expression performing a += operation for symmetric matrices
224  Need to have a separate functions to avoid to modify two times the off-diagonal
225  elements (i.e applying two times the expression)
226  Need to check whether creating a temporary object with the expression result
227  (like in op: A += A * B )
228  */
229  template <class T,
230  unsigned int D1, unsigned int D2,
231  class A>
232  struct PlusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
233  {
234  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs, const Expr<A,T,D1,D2, MatRepSym<T,D1> >& rhs)
235  {
236  if (! rhs.IsInUse(lhs.begin() ) ) {
237  unsigned int l = 0; // l span storage of sym matrices
238  for(unsigned int i=0; i<D1; ++i)
239  for(unsigned int j=0; j<=i; ++j) {
240  lhs.fRep.Array()[l] += rhs(i,j);
241  l++;
242  }
243  }
244  else {
245  T tmp[MatRepSym<T,D1>::kSize];
246  unsigned int l = 0;
247  for(unsigned int i=0; i<D1; ++i)
248  for(unsigned int j=0; j<=i; ++j) {
249  tmp[l] = rhs(i,j);
250  l++;
251  }
252  // += now using the temp object
253  for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] += tmp[i];
254  }
255  }
256  };
257  /**
258  Specialization for symmetrix += general : NOT Allowed operation
259  */
260  template <class T, unsigned int D1, unsigned int D2, class A>
261  struct PlusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
262  {
263  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
264  const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
265  {
266  STATIC_CHECK(0==1, Cannot_plusEqual_general_to_symmetric_matrix);
267  }
268  }; // struct PlusEquals
269 
270  //=========================================================================
271 
272  /**
273  Evaluate the expression performing a -= operation
274  Need to check whether creating a temporary object with the expression result
275  (like in op: A -= A * B )
276  */
277  template <class T, unsigned int D1, unsigned int D2, class A,
278  class R1, class R2>
279  struct MinusEquals
280  {
281  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const Expr<A,T,D1,D2,R2>& rhs)
282  {
283  if (! rhs.IsInUse(lhs.begin() ) ) {
284  unsigned int l = 0;
285  for(unsigned int i=0; i<D1; ++i)
286  for(unsigned int j=0; j<D2; ++j) {
287  lhs.fRep[l] -= rhs(i,j);
288  l++;
289  }
290  }
291  else {
292  T tmp[D1*D2];
293  unsigned int l = 0;
294  for(unsigned int i=0; i<D1; ++i)
295  for(unsigned int j=0; j<D2; ++j) {
296  tmp[l] = rhs(i,j);
297  l++;
298  }
299  // -= now using the temp object
300  for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] -= tmp[i];
301  }
302  }
303  };
304  /**
305  Specialization for symmetric matrices.
306  Evaluate the expression performing a -= operation for symmetric matrices
307  Need to have a separate functions to avoid to modify two times the off-diagonal
308  elements (i.e applying two times the expression)
309  Need to check whether creating a temporary object with the expression result
310  (like in op: A -= A + B )
311  */
312  template <class T,
313  unsigned int D1, unsigned int D2,
314  class A>
315  struct MinusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
316  {
317  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs, const Expr<A,T,D1,D2, MatRepSym<T,D1> >& rhs)
318  {
319  if (! rhs.IsInUse(lhs.begin() ) ) {
320  unsigned int l = 0; // l span storage of sym matrices
321  for(unsigned int i=0; i<D1; ++i)
322  for(unsigned int j=0; j<=i; ++j) {
323  lhs.fRep.Array()[l] -= rhs(i,j);
324  l++;
325  }
326  }
327  else {
328  T tmp[MatRepSym<T,D1>::kSize];
329  unsigned int l = 0;
330  for(unsigned int i=0; i<D1; ++i)
331  for(unsigned int j=0; j<=i; ++j) {
332  tmp[l] = rhs(i,j);
333  l++;
334  }
335  // -= now using the temp object
336  for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] -= tmp[i];
337  }
338  }
339  };
340 
341  /**
342  Specialization for symmetrix -= general : NOT Allowed operation
343  */
344  template <class T, unsigned int D1, unsigned int D2, class A>
345  struct MinusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
346  {
347  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
348  const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
349  {
350  STATIC_CHECK(0==1, Cannot_minusEqual_general_to_symmetric_matrix);
351  }
352  }; // struct MinusEquals
353 
354 
355  /** Structure to deal when a submatrix is placed in a matrix.
356  We have different cases according to the matrix representation
357  */
358  template <class T, unsigned int D1, unsigned int D2,
359  unsigned int D3, unsigned int D4,
360  class R1, class R2>
361  struct PlaceMatrix
362  {
363  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const SMatrix<T,D3,D4,R2>& rhs,
364  unsigned int row, unsigned int col) {
365 
366  assert(row+D3 <= D1 && col+D4 <= D2);
367  const unsigned int offset = row*D2+col;
368 
369  for(unsigned int i=0; i<D3*D4; ++i) {
370  lhs.fRep[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
371  }
372 
373  }
374  }; // struct PlaceMatrix
375 
376  template <class T, unsigned int D1, unsigned int D2,
377  unsigned int D3, unsigned int D4,
378  class A, class R1, class R2>
379  struct PlaceExpr {
380  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const Expr<A,T,D3,D4,R2>& rhs,
381  unsigned int row, unsigned int col) {
382 
383  assert(row+D3 <= D1 && col+D4 <= D2);
384  const unsigned int offset = row*D2+col;
385 
386  for(unsigned int i=0; i<D3*D4; ++i) {
387  lhs.fRep[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
388  }
389  }
390  }; // struct PlaceExpr
391 
392  // specialization for general matrix in symmetric matrices
393  template <class T, unsigned int D1, unsigned int D2,
394  unsigned int D3, unsigned int D4 >
395  struct PlaceMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
396  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
397  const SMatrix<T,D3,D4,MatRepStd<T,D3,D4> >& ,
398  unsigned int , unsigned int )
399  {
400  STATIC_CHECK(0==1, Cannot_Place_Matrix_general_in_symmetric_matrix);
401  }
402  }; // struct PlaceMatrix
403 
404  // specialization for general expression in symmetric matrices
405  template <class T, unsigned int D1, unsigned int D2,
406  unsigned int D3, unsigned int D4, class A >
407  struct PlaceExpr<T, D1, D2, D3, D4, A, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
408  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
409  const Expr<A,T,D3,D4,MatRepStd<T,D3,D4> >& ,
410  unsigned int , unsigned int )
411  {
412  STATIC_CHECK(0==1, Cannot_Place_Matrix_general_in_symmetric_matrix);
413  }
414  }; // struct PlaceExpr
415 
416  // specialization for symmetric matrix in symmetric matrices
417 
418  template <class T, unsigned int D1, unsigned int D2,
419  unsigned int D3, unsigned int D4 >
420  struct PlaceMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepSym<T,D3> > {
421  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
422  const SMatrix<T,D3,D4,MatRepSym<T,D3> >& rhs,
423  unsigned int row, unsigned int col )
424  {
425  // can work only if placed on the diagonal
426  assert(row == col);
427 
428  for(unsigned int i=0; i<D3; ++i) {
429  for(unsigned int j=0; j<=i; ++j)
430  lhs.fRep(row+i,col+j) = rhs(i,j);
431  }
432  }
433  }; // struct PlaceMatrix
434 
435  // specialization for symmetric expression in symmetric matrices
436  template <class T, unsigned int D1, unsigned int D2,
437  unsigned int D3, unsigned int D4, class A >
438  struct PlaceExpr<T, D1, D2, D3, D4, A, MatRepSym<T,D1>, MatRepSym<T,D3> > {
439  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
440  const Expr<A,T,D3,D4,MatRepSym<T,D3> >& rhs,
441  unsigned int row, unsigned int col )
442  {
443  // can work only if placed on the diagonal
444  assert(row == col);
445 
446  for(unsigned int i=0; i<D3; ++i) {
447  for(unsigned int j=0; j<=i; ++j)
448  lhs.fRep(row+i,col+j) = rhs(i,j);
449  }
450  }
451  }; // struct PlaceExpr
452 
453 
454 
455  /** Structure for getting sub matrices
456  We have different cases according to the matrix representations
457  */
458  template <class T, unsigned int D1, unsigned int D2,
459  unsigned int D3, unsigned int D4,
460  class R1, class R2>
461  struct RetrieveMatrix
462  {
463  static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const SMatrix<T,D3,D4,R2>& rhs,
464  unsigned int row, unsigned int col) {
465  STATIC_CHECK( D1 <= D3,Smatrix_nrows_too_small);
466  STATIC_CHECK( D2 <= D4,Smatrix_ncols_too_small);
467 
468  assert(row + D1 <= D3);
469  assert(col + D2 <= D4);
470 
471  for(unsigned int i=0; i<D1; ++i) {
472  for(unsigned int j=0; j<D2; ++j)
473  lhs(i,j) = rhs(i+row,j+col);
474  }
475  }
476  }; // struct RetrieveMatrix
477 
478  // specialization for getting symmetric matrices from general matrices (MUST fail)
479  template <class T, unsigned int D1, unsigned int D2,
480  unsigned int D3, unsigned int D4 >
481  struct RetrieveMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
482  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
483  const SMatrix<T,D3,D4,MatRepStd<T,D3,D4> >& ,
484  unsigned int , unsigned int )
485  {
486  STATIC_CHECK(0==1, Cannot_Sub_Matrix_symmetric_in_general_matrix);
487  }
488  }; // struct RetrieveMatrix
489 
490  // specialization for getting symmetric matrices from symmetric matrices (OK if row == col)
491  template <class T, unsigned int D1, unsigned int D2,
492  unsigned int D3, unsigned int D4 >
493  struct RetrieveMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepSym<T,D3> > {
494  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
495  const SMatrix<T,D3,D4,MatRepSym<T,D3> >& rhs,
496  unsigned int row, unsigned int col )
497  {
498  STATIC_CHECK( D1 <= D3,Smatrix_dimension1_too_small);
499  // can work only if placed on the diagonal
500  assert(row == col);
501  assert(row + D1 <= D3);
502 
503  for(unsigned int i=0; i<D1; ++i) {
504  for(unsigned int j=0; j<=i; ++j)
505  lhs(i,j) = rhs(i+row,j+col );
506  }
507  }
508 
509  }; // struct RetrieveMatrix
510 
511  /**
512  Structure for assignment to a general matrix from iterator.
513  Optionally a check is done that iterator size
514  is not larger than matrix size
515  */
516  template <class T, unsigned int D1, unsigned int D2, class R>
517  struct AssignItr {
518  template<class Iterator>
519  static void Evaluate(SMatrix<T,D1,D2,R>& lhs, Iterator begin, Iterator end,
520  bool triang, bool lower,bool check=true) {
521  // require size match exactly (better)
522 
523  if (triang) {
524  Iterator itr = begin;
525  if (lower) {
526  for (unsigned int i = 0; i < D1; ++i)
527  for (unsigned int j =0; j <= i; ++j) {
528  // we assume iterator is well bounded within matrix
529  lhs.fRep[i*D2+j] = *itr++;
530  }
531 
532  }
533  else { // upper
534  for (unsigned int i = 0; i < D1; ++i)
535  for (unsigned int j = i; j <D2; ++j) {
536  if (itr != end)
537  lhs.fRep[i*D2+j] = *itr++;
538  else
539  return;
540  }
541 
542  }
543  }
544  // case of filling the full matrix
545  else {
546  if (check) assert( begin + R::kSize == end);
547  // copy directly the elements
548  std::copy(begin, end, lhs.fRep.Array() );
549  }
550  }
551 
552  }; // struct AssignItr
553 
554  /**
555  Specialized structure for assignment to a symmetrix matrix from iterator.
556  Optionally a check is done that iterator size
557  is the same as the matrix size
558  */
559  template <class T, unsigned int D1, unsigned int D2>
560  struct AssignItr<T, D1, D2, MatRepSym<T,D1> > {
561  template<class Iterator>
562  static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs, Iterator begin, Iterator end, bool , bool lower, bool check = true) {
563 
564  if (lower) {
565  if (check) {
566  assert(begin+ static_cast< int>( MatRepSym<T,D1>::kSize) == end);
567  }
568  std::copy(begin, end, lhs.fRep.Array() );
569  }
570  else {
571  Iterator itr = begin;
572  for (unsigned int i = 0; i < D1; ++i)
573  for (unsigned int j = i; j <D2; ++j) {
574  if (itr != end)
575  lhs(i,j) = *itr++;
576  else
577  return;
578  }
579  }
580  }
581 
582  }; // struct AssignItr
583 
584 
585 } // namespace Math
586 
587 } // namespace ROOT
588 
589 #endif // MATH_HELPEROPS_H