FEDRA emulsion software from the OPERA Collaboration
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Dfinv.hh File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

template<class Matrix , unsigned int n, unsigned int idim>
bool Dfinv (Matrix &rhs, unsigned int *ir)
 

Function Documentation

◆ Dfinv()

template<class Matrix , unsigned int n, unsigned int idim>
bool Dfinv ( Matrix &  rhs,
unsigned int *  ir 
)

Dfinv. Function to compute the inverse of a square matrix ($A^{-1}$) of dimension $idim$ and order $n$. The routine Dfactir must be called before Dfinv!

Author
T. Glebe
36 {
37#ifdef XXX
38 if (idim < n || n <= 0 || n==1) {
39 return false;
40 }
41#endif
42
43 typename Matrix::value_type* a = rhs.Array();
44
45 /* Local variables */
46 static unsigned int nxch, i, j, k, m, ij;
47 static unsigned int im2, nm1, nmi;
48 static typename Matrix::value_type s31, s34, ti;
49
50 /* Parameter adjustments */
51 a -= idim + 1;
52 --ir;
53
54 /* Function Body */
55
56 /* finv.inc */
57
58 a[idim + 2] = -a[(idim << 1) + 2] * (a[idim + 1] * a[idim + 2]);
59 a[(idim << 1) + 1] = -a[(idim << 1) + 1];
60
61 if (n != 2) {
62 for (i = 3; i <= n; ++i) {
63 const unsigned int ii = i * idim;
64 const unsigned int iii = i + ii;
65 const unsigned int imi = ii - idim;
66 const unsigned int iimi = i + imi;
67 im2 = i - 2;
68 for (j = 1; j <= im2; ++j) {
69 const unsigned int ji = j * idim;
70 const unsigned int jii = j + ii;
71 s31 = 0.;
72 for (k = j; k <= im2; ++k) {
73 s31 += a[k + ji] * a[i + k * idim];
74 a[jii] += a[j + (k + 1) * idim] * a[k + 1 + ii];
75 } // for k
76 a[i + ji] = -a[iii] * (a[i - 1 + ji] * a[iimi] + s31);
77 a[jii] *= -1;
78 } // for j
79 a[iimi] = -a[iii] * (a[i - 1 + imi] * a[iimi]);
80 a[i - 1 + ii] *= -1;
81 } // for i
82 } // if n!=2
83
84 nm1 = n - 1;
85 for (i = 1; i <= nm1; ++i) {
86 const unsigned int ii = i * idim;
87 nmi = n - i;
88 for (j = 1; j <= i; ++j) {
89 const unsigned int ji = j * idim;
90 const unsigned int iji = i + ji;
91 for (k = 1; k <= nmi; ++k) {
92 a[iji] += a[i + k + ji] * a[i + (i + k) * idim];
93 } // for k
94 } // for j
95
96 for (j = 1; j <= nmi; ++j) {
97 const unsigned int ji = j * idim;
98 s34 = 0.;
99 for (k = j; k <= nmi; ++k) {
100 s34 += a[i + k + ii + ji] * a[i + (i + k) * idim];
101 } // for k
102 a[i + ii + ji] = s34;
103 } // for j
104 } // for i
105
106 nxch = ir[n];
107 if (nxch == 0) {
108 return true;
109 }
110
111 for (m = 1; m <= nxch; ++m) {
112 k = nxch - m + 1;
113 ij = ir[k];
114 i = ij / 4096;
115 j = ij % 4096;
116 const unsigned int ii = i * idim;
117 const unsigned int ji = j * idim;
118 for (k = 1; k <= n; ++k) {
119 ti = a[k + ii];
120 a[k + ii] = a[k + ji];
121 a[k + ji] = ti;
122 } // for k
123 } // for m
124
125 return true;
126} // Dfinv
EdbTrackP * ti[10]
Definition: RecDispMC_Profiles.C:54
void a()
Definition: check_aligned.C:59