ROOT   Reference Guide
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
32namespace ROOT {
33
34 namespace Math {
35
36
37
38
39/** Dfinv.
40 Function to compute the inverse of a square matrix (\f$A^{-1}\f$) of
41 dimension \f$idim\f$ and order \f$n\f$. The routine Dfactir must be called
42 before Dfinv!
43
44 @author T. Glebe
45*/
46template <class Matrix, unsigned int n, unsigned int idim>
47bool 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 */
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
bool Dfinv(Matrix &rhs, unsigned int *ir)
Dfinv.
Definition: Dfinv.h:47
VSD Structures.
Definition: StringConv.hxx:21
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12