ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Dfinv.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dfinv
5 #define ROOT_Math_Dfinv
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 03. Apr 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: Matrix inversion
23 // Code was taken from CERNLIB::kernlib dfinv function, translated
24 // from FORTRAN to C++ and optimized.
25 //
26 // changes:
27 // 03 Apr 2001 (TG) creation
28 //
29 // ********************************************************************
30 
31 
32 namespace ROOT {
33 
34  namespace Math {
35 
36 
37 
38 
39 /** Dfinv.
40  Function to compute the inverse of a square matrix ($A^{-1}$) of
41  dimension $idim$ and order $n$. The routine Dfactir must be called
42  before Dfinv!
43 
44  @author T. Glebe
45 */
46 template <class Matrix, unsigned int n, unsigned int idim>
47 bool Dfinv(Matrix& rhs, unsigned int* ir) {
48 #ifdef XXX
49  if (idim < n || n <= 0 || n==1) {
50  return false;
51  }
52 #endif
53 
54  typename Matrix::value_type* a = rhs.Array();
55 
56  /* Local variables */
57  unsigned int nxch, i, j, k, m, ij;
58  unsigned int im2, nm1, nmi;
59  typename Matrix::value_type s31, s34, ti;
60 
61  /* Parameter adjustments */
62  a -= idim + 1;
63  --ir;
64 
65  /* Function Body */
66 
67  /* finv.inc */
68 
69  a[idim + 2] = -a[(idim << 1) + 2] * (a[idim + 1] * a[idim + 2]);
70  a[(idim << 1) + 1] = -a[(idim << 1) + 1];
71 
72  if (n != 2) {
73  for (i = 3; i <= n; ++i) {
74  const unsigned int ii = i * idim;
75  const unsigned int iii = i + ii;
76  const unsigned int imi = ii - idim;
77  const unsigned int iimi = i + imi;
78  im2 = i - 2;
79  for (j = 1; j <= im2; ++j) {
80  const unsigned int ji = j * idim;
81  const unsigned int jii = j + ii;
82  s31 = 0.;
83  for (k = j; k <= im2; ++k) {
84  s31 += a[k + ji] * a[i + k * idim];
85  a[jii] += a[j + (k + 1) * idim] * a[k + 1 + ii];
86  } // for k
87  a[i + ji] = -a[iii] * (a[i - 1 + ji] * a[iimi] + s31);
88  a[jii] *= -1;
89  } // for j
90  a[iimi] = -a[iii] * (a[i - 1 + imi] * a[iimi]);
91  a[i - 1 + ii] *= -1;
92  } // for i
93  } // if n!=2
94 
95  nm1 = n - 1;
96  for (i = 1; i <= nm1; ++i) {
97  const unsigned int ii = i * idim;
98  nmi = n - i;
99  for (j = 1; j <= i; ++j) {
100  const unsigned int ji = j * idim;
101  const unsigned int iji = i + ji;
102  for (k = 1; k <= nmi; ++k) {
103  a[iji] += a[i + k + ji] * a[i + (i + k) * idim];
104  } // for k
105  } // for j
106 
107  for (j = 1; j <= nmi; ++j) {
108  const unsigned int ji = j * idim;
109  s34 = 0.;
110  for (k = j; k <= nmi; ++k) {
111  s34 += a[i + k + ii + ji] * a[i + (i + k) * idim];
112  } // for k
113  a[i + ii + ji] = s34;
114  } // for j
115  } // for i
116 
117  nxch = ir[n];
118  if (nxch == 0) {
119  return true;
120  }
121 
122  for (m = 1; m <= nxch; ++m) {
123  k = nxch - m + 1;
124  ij = ir[k];
125  i = ij / 4096;
126  j = ij % 4096;
127  const unsigned int ii = i * idim;
128  const unsigned int ji = j * idim;
129  for (k = 1; k <= n; ++k) {
130  ti = a[k + ii];
131  a[k + ii] = a[k + ji];
132  a[k + ji] = ti;
133  } // for k
134  } // for m
135 
136  return true;
137 } // Dfinv
138 
139 
140  } // namespace Math
141 
142 } // namespace ROOT
143 
144 
145 
146 #endif /* ROOT_Math_Dfinv */
TArc * a
Definition: textangle.C:12
bool Dfinv(Matrix &rhs, unsigned int *ir)
Dfinv.
Definition: Dfinv.h:47
TMarker * m
Definition: textangle.C:8
const Int_t n
Definition: legend1.C:16
int ii
Definition: hprod.C:34