root/Math/Dsfact.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_Dsfact
0005 #define ROOT_Math_Dsfact
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   22. Mar 2001
0013 //
0014 // author:    Thorsten Glebe
0015 //            HERA-B Collaboration
0016 //            Max-Planck-Institut fuer Kernphysik
0017 //            Saupfercheckweg 1
0018 //            69117 Heidelberg
0019 //            Germany
0020 //            E-mail: T.Glebe@mpi-hd.mpg.de
0021 //
0022 // Description: Determinant of a symmetric, positive definite matrix.
0023 //              Code was taken from CERNLIB::kernlib dsfact function, translated
0024 //              from FORTRAN to C++ and optimized.
0025 //
0026 // changes:
0027 // 22 Mar 2001 (TG) creation
0028 // 18 Apr 2001 (TG) removed internal copying of array.
0029 //
0030 // ********************************************************************
0031
0032 #include "Math/MatrixRepresentationsStatic.h"
0033
0034 namespace ROOT {
0035
0036   namespace Math {
0037
0038
0039
0040
0041 /** Dsfact.
0042     Compute determinant of a symmetric, positive definite matrix of dimension
0043     \f$idim\f$ and order \f$n\f$.
0044
0045     @author T. Glebe
0046 */
0047
0048 template <unsigned int n, unsigned int idim =n>
0049 class SDeterminant {
0050
0051 public:
0052 template <class T>
0053 static bool Dsfact(MatRepStd<T,n,idim>& rhs, T& det) {
0054
0055 #ifdef XXX
0056   /* Function Body */
0057   if (idim < n || n <= 0) {
0058     return false;
0059   }
0060 #endif
0061
0062 #ifdef OLD_IMPL
0063   typename MatrixRep::value_type* a = rhs.Array();
0064 #endif
0065
0066 #ifdef XXX
0067   const typename MatrixRep::value_type* A = rhs.Array();
0068   typename MatrixRep::value_type array[MatrixRep::kSize];
0069   typename MatrixRep::value_type* a = array;
0070
0071   // copy contents of matrix to working place
0072   for(unsigned int i=0; i<MatrixRep::kSize; ++i) {
0073     array[i] = A[i];
0074   }
0075 #endif
0076
0077   /* Local variables */
0078   unsigned int i, j, l;
0079
0080   /* Parameter adjustments */
0081   //  a -= idim + 1;
0082   const int arrayOffset = -(idim+1);
0083   /* sfactd.inc */
0084   det = 1.;
0085   for (j = 1; j <= n; ++j) {
0086     const unsigned int ji = j * idim;
0087     const unsigned int jj = j + ji;
0088
0089     if (rhs[jj + arrayOffset] <= 0.) {
0090       det = 0.;
0091       return false;
0092     }
0093
0094     const unsigned int jp1 = j + 1;
0095     const unsigned int jpi = jp1 * idim;
0096
0097     det *= rhs[jj + arrayOffset];
0098     rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
0099
0100     for (l = jp1; l <= n; ++l) {
0101       rhs[j + l * idim + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ji + arrayOffset];
0102
0103       const unsigned int lj = l + jpi;
0104
0105       for (i = 1; i <= j; ++i) {
0106          rhs[lj + arrayOffset] -= rhs[l + i * idim + arrayOffset] * rhs[i + jpi + arrayOffset];
0107       } // for i
0108     } // for l
0109   } // for j
0110
0111   return true;
0112 } // end of Dsfact
0113
0114
0115    // t.b.d re-implement methods for symmetric
0116   // symmetric function (copy in a general  one)
0117   template <class T>
0118   static bool Dsfact(MatRepSym<T,n> & rhs,  T & det) {
0119     // not very efficient but need to re-do Dsinv for new storage of
0120     // symmetric matrices
0121     MatRepStd<T,n> tmp;
0122     for (unsigned int i = 0; i< n*n; ++i)
0123       tmp[i] = rhs[i];
0124     if (!  SDeterminant<n>::Dsfact(tmp,det) ) return false;
0125 //     // recopy the data
0126 //     for (int i = 0; i< idim*n; ++i)
0127 //       rhs[i] = tmp[i];
0128
0129     return true;
0130   }
0131
0132
0133 };  // end of class Sdeterminant
0134
0135   }  // namespace Math
0136
0137 }  // namespace ROOT
0138
0139 #endif  /* ROOT_Math_Dsfact */
0140