root/Math/Dsinv.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef  ROOT_Math_Dsinv
0005 #define  ROOT_Math_Dsinv
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: Inversion of a symmetric, positive definite matrix.
0023 //              Code was taken from CERNLIB::kernlib dsinv function, translated
0024 //              from FORTRAN to C++ and optimized.
0025 //
0026 // changes:
0027 // 22 Mar 2001 (TG) creation
0028 //
0029 // ********************************************************************
0030
0031 #include "Math/SMatrixDfwd.h"
0032
0033 namespace ROOT {
0034
0035   namespace Math {
0036
0037 /** Dsinv.
0038     Compute inverse of a symmetric, positive definite matrix of dimension
0039     \f$idim\f$ and order \f$n\f$.
0040
0041     @author T. Glebe
0042 */
0043 template <class T, int n, int idim>
0044 class SInverter
0045 {
0046
0047 public:
0048   template <class MatrixRep>
0049   inline static bool Dsinv(MatrixRep& rhs) {
0050
0051     /* Local variables */
0052     int i, j, k, l;
0053     T s31, s32;
0054     int jm1, jp1;
0055
0056     /* Parameter adjustments */
0057     const int arrayOffset = -1*(idim + 1);
0058
0059
0060     /* Function Body */
0061     if (idim < n || n <= 1) {
0062       return false;
0063     }
0064
0065     /* sfact.inc */
0066     for (j = 1; j <= n; ++j) {
0067       const int ja  = j * idim;
0068       const int jj  = j + ja;
0069       const int ja1 = ja + idim;
0070
0071
0072       if (rhs[jj + arrayOffset] <= 0.) { return false; }
0073       rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
0074       if (j == n) { break; }
0075
0076       for (l = j + 1; l <= n; ++l) {
0077         rhs[j + (l * idim) + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ja + arrayOffset];
0078         const int lj = l + ja1;
0079         for (i = 1; i <= j; ++i) {
0080           rhs[lj + arrayOffset] -= rhs[l + (i * idim)  + arrayOffset] * rhs[i + ja1 + arrayOffset];
0081         }
0082       }
0083     }
0084
0085     /* sfinv.inc */
0086     // idim << 1 is equal to idim * 2
0087     // compiler will compute the arguments!
0088     rhs[((idim << 1) + 1) + arrayOffset] = -rhs[((idim << 1) + 1) + arrayOffset];
0089     rhs[idim + 2 + arrayOffset] = rhs[((idim << 1)) + 1 + arrayOffset] * rhs[((idim << 1)) + 2 + arrayOffset];
0090
0091     if(n > 2) {
0092
0093       for (j = 3; j <= n; ++j) {
0094         const int jm2 = j - 2;
0095         const int ja = j * idim;
0096         const int jj = j + ja;
0097         const int j1 = j - 1 + ja;
0098
0099         for (k = 1; k <= jm2; ++k) {
0100           s31 = rhs[k + ja + arrayOffset];
0101
0102           for (i = k; i <= jm2; ++i) {
0103             s31 += rhs[k + ((i + 1) * idim) + arrayOffset] * rhs[i + 1 + ja + arrayOffset];
0104           } // for i
0105           rhs[k + ja + arrayOffset] = -s31;
0106           rhs[j + (k * idim) + arrayOffset] = -s31 * rhs[jj + arrayOffset];
0107         } // for k
0108         rhs[j1 + arrayOffset] *= -1;
0109         //      rhs[j1] = -rhs[j1];
0110         rhs[jj - idim + arrayOffset] = rhs[j1 + arrayOffset] * rhs[jj + arrayOffset];
0111       } // for j
0112     } // if (n>2)
0113
0114     j = 1;
0115     do {
0116       const int jad = j * idim;
0117       const int jj = j + jad;
0118
0119       jp1 = j + 1;
0120       for (i = jp1; i <= n; ++i) {
0121         rhs[jj + arrayOffset] += rhs[j + (i * idim) + arrayOffset] * rhs[i + jad + arrayOffset];
0122       } // for i
0123
0124       jm1 = j;
0125       j = jp1;
0126       const int ja = j * idim;
0127
0128       for (k = 1; k <= jm1; ++k) {
0129         s32 = 0.;
0130         for (i = j; i <= n; ++i) {
0131           s32 += rhs[k + (i * idim) + arrayOffset] * rhs[i + ja + arrayOffset];
0132         } // for i
0133         //rhs[k + ja + arrayOffset] = rhs[j + (k * idim) + arrayOffset] = s32;
0134         rhs[k + ja + arrayOffset] = s32;
0135       } // for k
0136     } while(j < n);
0137
0138     return true;
0139   }
0140
0141
0142     // for symmetric matrices
0143
0144   static bool Dsinv(MatRepSym<T,n> & rhs) {
0145     // not very efficient but need to re-do Dsinv for new storage of
0146     // symmetric matrices
0147     MatRepStd<T,n,n> tmp;
0148     for (int i = 0; i< n*n; ++i)
0149       tmp[i] = rhs[i];
0150     // call dsinv
0151     if (! SInverter<T,n,n>::Dsinv(tmp) ) return false;
0152     //if (! Inverter<n>::Dinv(tmp) ) return false;
0153     // recopy the data
0154     for (int i = 0; i< n*n; ++i)
0155       rhs[i] = tmp[i];
0156
0157     return true;
0158
0159   }
0160
0161 }; // end of Dsinv
0162
0163
0164
0165   }  // namespace Math
0166
0167 }  // namespace ROOT
0168
0169
0170 #endif  /* ROOT_Math_Dsinv */