Logo ROOT  
Reference Guide
Dsinv.h
Go to the documentation of this file.
1// @(#)root/smatrix:$Id$
2// Authors: T. Glebe, L. Moneta 2005
3
4#ifndef ROOT_Math_Dsinv
5#define ROOT_Math_Dsinv
6// ********************************************************************
7//
8// source:
9//
10// type: source code
11//
12// created: 22. Mar 2001
13//
14// author: Thorsten Glebe
15// HERA-B Collaboration
16// Max-Planck-Institut fuer Kernphysik
17// Saupfercheckweg 1
18// 69117 Heidelberg
19// Germany
20// E-mail: T.Glebe@mpi-hd.mpg.de
21//
22// Description: Inversion of a symmetric, positive definite matrix.
23// Code was taken from CERNLIB::kernlib dsinv function, translated
24// from FORTRAN to C++ and optimized.
25//
26// changes:
27// 22 Mar 2001 (TG) creation
28//
29// ********************************************************************
30
31#include "Math/SMatrixDfwd.h"
32
33namespace ROOT {
34
35 namespace Math {
36
37/** Dsinv.
38 Compute inverse of a symmetric, positive definite matrix of dimension
39 \f$idim\f$ and order \f$n\f$.
40
41 @author T. Glebe
42*/
43template <class T, int n, int idim>
45{
46
47public:
48 template <class MatrixRep>
49 inline static bool Dsinv(MatrixRep& rhs) {
50
51 /* Local variables */
52 int i, j, k, l;
53 T s31, s32;
54 int jm1, jp1;
55
56 /* Parameter adjustments */
57 const int arrayOffset = -1*(idim + 1);
58
59
60 /* Function Body */
61 if (idim < n || n <= 1) {
62 return false;
63 }
64
65 /* sfact.inc */
66 for (j = 1; j <= n; ++j) {
67 const int ja = j * idim;
68 const int jj = j + ja;
69 const int ja1 = ja + idim;
70
71
72 if (rhs[jj + arrayOffset] <= 0.) { return false; }
73 rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
74 if (j == n) { break; }
75
76 for (l = j + 1; l <= n; ++l) {
77 rhs[j + (l * idim) + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ja + arrayOffset];
78 const int lj = l + ja1;
79 for (i = 1; i <= j; ++i) {
80 rhs[lj + arrayOffset] -= rhs[l + (i * idim) + arrayOffset] * rhs[i + ja1 + arrayOffset];
81 }
82 }
83 }
84
85 /* sfinv.inc */
86 // idim << 1 is equal to idim * 2
87 // compiler will compute the arguments!
88 rhs[((idim << 1) + 1) + arrayOffset] = -rhs[((idim << 1) + 1) + arrayOffset];
89 rhs[idim + 2 + arrayOffset] = rhs[((idim << 1)) + 1 + arrayOffset] * rhs[((idim << 1)) + 2 + arrayOffset];
90
91 if(n > 2) {
92
93 for (j = 3; j <= n; ++j) {
94 const int jm2 = j - 2;
95 const int ja = j * idim;
96 const int jj = j + ja;
97 const int j1 = j - 1 + ja;
98
99 for (k = 1; k <= jm2; ++k) {
100 s31 = rhs[k + ja + arrayOffset];
101
102 for (i = k; i <= jm2; ++i) {
103 s31 += rhs[k + ((i + 1) * idim) + arrayOffset] * rhs[i + 1 + ja + arrayOffset];
104 } // for i
105 rhs[k + ja + arrayOffset] = -s31;
106 rhs[j + (k * idim) + arrayOffset] = -s31 * rhs[jj + arrayOffset];
107 } // for k
108 rhs[j1 + arrayOffset] *= -1;
109 // rhs[j1] = -rhs[j1];
110 rhs[jj - idim + arrayOffset] = rhs[j1 + arrayOffset] * rhs[jj + arrayOffset];
111 } // for j
112 } // if (n>2)
113
114 j = 1;
115 do {
116 const int jad = j * idim;
117 const int jj = j + jad;
118
119 jp1 = j + 1;
120 for (i = jp1; i <= n; ++i) {
121 rhs[jj + arrayOffset] += rhs[j + (i * idim) + arrayOffset] * rhs[i + jad + arrayOffset];
122 } // for i
123
124 jm1 = j;
125 j = jp1;
126 const int ja = j * idim;
127
128 for (k = 1; k <= jm1; ++k) {
129 s32 = 0.;
130 for (i = j; i <= n; ++i) {
131 s32 += rhs[k + (i * idim) + arrayOffset] * rhs[i + ja + arrayOffset];
132 } // for i
133 //rhs[k + ja + arrayOffset] = rhs[j + (k * idim) + arrayOffset] = s32;
134 rhs[k + ja + arrayOffset] = s32;
135 } // for k
136 } while(j < n);
137
138 return true;
139 }
140
141
142 // for symmetric matrices
143
144 static bool Dsinv(MatRepSym<T,n> & rhs) {
145 // not very efficient but need to re-do Dsinv for new storage of
146 // symmetric matrices
148 for (int i = 0; i< n*n; ++i)
149 tmp[i] = rhs[i];
150 // call dsinv
151 if (! SInverter<T,n,n>::Dsinv(tmp) ) return false;
152 //if (! Inverter<n>::Dinv(tmp) ) return false;
153 // recopy the data
154 for (int i = 0; i< n*n; ++i)
155 rhs[i] = tmp[i];
156
157 return true;
158
159 }
160
161}; // end of Dsinv
162
163
164
165 } // namespace Math
166
167} // namespace ROOT
168
169
170#endif /* ROOT_Math_Dsinv */
Expression wrapper class for Matrix objects.
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
static bool Dsinv(MatRepSym< T, n > &rhs)
Definition: Dsinv.h:144
static bool Dsinv(MatrixRep &rhs)
Definition: Dsinv.h:49
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
TLine l
Definition: textangle.C:4