Garfield++ v2r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
inverse.cpp
Go to the documentation of this file.
1#include "wcpplib/matrix/abs_inverse.h"
2#include "wcpplib/matrix/inverse.h"
3/*
4Copyright (c) 2001 Igor B. Smirnov
5
6The file can be used, copied, modified, and distributed
7according to the terms of GNU Lesser General Public License version 2.1
8as published by the Free Software Foundation,
9and provided that the above copyright notice, this permission notice,
10and notices about any modifications of the original text
11appear in all copies and in supporting documentation.
12The file is provided "as is" without express or implied warranty.
13*/
14
15namespace Heed {
16
17void inverse_DynArr_prot(const DynArr<DoubleAc>& mi, DynArr<DoubleAc>& mr,
18 int& szero, int& serr, int s_stop) {
19 mfunname(
20 "void inverse_DynArr_prot(const DynArr<DoubleAc>& mi, "
21 "DynArr<DoubleAc>& mr, int& s_zero, int& serr, int s_stop)");
22 // mcout<<"inverse_DynArr_prot:\n";
23 const DynLinArr<long>& miqel(mi.get_qel());
24 check_econd11(miqel.get_qel(), != 2, mcerr);
25 check_econd11(miqel[0], <= 0, mcerr);
26 check_econd12(miqel[0], !=, miqel[1], mcerr);
27 serr = 0;
28 szero = 0;
29 long q = miqel[0];
30 mr = DynArr<DoubleAc>(q, q);
31 if (q == 1) {
32 if (mi.ac(0, 0).get() == 0.0) {
33 szero = 1;
34 serr = 1;
35 return;
36 }
37 mr.ac(0, 0) = 1.0 / mi.ac(0, 0);
38 if (fabs(mr.ac(0, 0)).left_limit() == 0) {
39 serr = 1;
40 }
41 return;
42 }
43 DynArr<DoubleAc> mii(mi);
44 mr.assignAll(0.0);
45 for (long n = 0; n < q; n++) mr.ac(n, n) = DoubleAc(1.0);
46
47 for (long nr = 0; nr < q; nr++) {
48 // Iprintn(mcout, nr);
49 long nmax = 0;
50 DoubleAc d(0.0);
51 for (long nr1 = nr; nr1 < q; nr1++) {
52 if (fabs(mii.ac(nr1, nr)) > d) {
53 d = fabs(mii.ac(nr1, nr));
54 nmax = nr1;
55 }
56 }
57 // mcout<<"d="<<d<<'\n';
58 // mcout<<"nmax="<<nmax<<'\n';
59 if (d.get() == 0.0) {
60 szero = 1;
61 serr = 1;
62 return;
63 }
64 if (d.left_limit() == 0) {
65 serr = 1;
66 if (s_stop == 1) return;
67 }
68 if (nmax > nr) {
69 for (long nc = nr; nc < q; nc++) {
70 DoubleAc t(mii.ac(nr, nc));
71 mii.ac(nr, nc) = mii.ac(nmax, nc);
72 mii.ac(nmax, nc) = t;
73 }
74 for (long nc = 0; nc < q; nc++) {
75 DoubleAc t(mr.ac(nr, nc));
76 mr.ac(nr, nc) = mr.ac(nmax, nc);
77 mr.ac(nmax, nc) = t;
78 }
79 }
80 DoubleAc t = mii.ac(nr, nr);
81 for (long nr1 = 0; nr1 < q; nr1++) {
82 if (nr1 != nr) {
83 DoubleAc k(mii.ac(nr1, nr) / t);
84 // mcout<<"nr1="<<nr1<<" nr="<<nr<<'\n';
85 // mcout<<"k="<<k<<'\n';
86 for (long nc = nr; nc < q; nc++) {
87 mii.ac(nr1, nc) -= k * mii.ac(nr, nc);
88 }
89 for (long nc = 0; nc < q; nc++) {
90 mr.ac(nr1, nc) -= k * mr.ac(nr, nc);
91 }
92 }
93 }
94 for (long nc = nr; nc < q; nc++) mii.ac(nr, nc) /= t;
95 for (long nc = 0; nc < q; nc++) mr.ac(nr, nc) /= t;
96 }
97}
98
99void inverse_DynArr(const DynArr<double>& mi, DynArr<double>& mr, int& serr) {
100 mfunname(
101 "void inverse_DynArr(const DynArr<double>& mi, DynArr<double>& mr, "
102 "int& serr)");
103 const DynLinArr<long>& miqel(mi.get_qel());
104 check_econd11(miqel.get_qel(), != 2, mcerr);
105 check_econd11(miqel[0], <= 0, mcerr);
106 check_econd12(miqel[0], !=, miqel[1], mcerr);
107 serr = 0;
108 long q = miqel[0];
109 mr = DynArr<double>(q, q);
110 if (q == 1) {
111 if (mi.ac(0, 0) == 0.0) {
112 serr = 1;
113 return;
114 }
115 mr.ac(0, 0) = 1.0 / mi.ac(0, 0);
116 return;
117 }
118 DynArr<DoubleAc> mii;
119 DynArr<DoubleAc> mrr(q, q);
120 copy_DynArr(mi, mii);
121 mrr.assignAll(0.0);
122 for (long n = 0; n < miqel[0]; n++) mrr.ac(n, n) = 1.0;
123 int szero;
124 inverse_DynArr_prot(mii, mrr, szero, serr);
125 copy_DynArr(mrr, mr);
126}
127
128void inverse_DynArr(const DynArr<DoubleAc>& mi, DynArr<DoubleAc>& mr1,
129 int& szero, int& serr1, DynArr<DoubleAc>& mr2, int& serr2) {
130 mfunname(
131 "void inverse_DynArr(const DynArr<DoubleAc>& mi, DynArr<DoubleAc>& "
132 "mr1, int& szero, int& serr1, DynArr<DoubleAc>& mr2, int& serr2)");
133 const DynLinArr<long>& miqel(mi.get_qel());
134 check_econd11(miqel.get_qel(), != 2, mcerr);
135 check_econd11(miqel[0], <= 0, mcerr);
136 check_econd12(miqel[0], !=, miqel[1], mcerr);
137 serr1 = 0;
138 serr2 = 0;
139 long q = miqel[0];
140 // mr = DynArr<DoubleAc>(miqel[0], miqel[0]);
141 if (q == 1) {
142 if (mi.ac(0, 0).get() == 0.0) {
143 serr1 = 1;
144 return;
145 }
146 mr1 = DynArr<DoubleAc>(q, q);
147 mr2 = DynArr<DoubleAc>(q, q);
148 mr1.ac(0, 0) = 1.0 / mi.ac(0, 0).get();
149 mr2.ac(0, 0) = DoubleAc(1.0) / mi.ac(0, 0);
150 if (fabs(mr2.ac(0, 0)).left_limit() == 0.0) serr2 = 1;
151 return;
152 }
153
154 DynArr<DoubleAc> mii(q, q);
155 for (long n1 = 0; n1 < q; n1++) {
156 for (long n2 = 0; n2 < q; n2++) {
157 mii.ac(n1, n2) = double(mi.ac(n1, n2));
158 }
159 }
160 DynArr<DoubleAc> mrr(q, q);
161 mrr.assignAll(0.0);
162 for (long n = 0; n < q; n++) mrr.ac(n, n) = 1.0;
163 inverse_DynArr_prot(mii, mrr, szero, serr1, 0);
164 copy_DynArr(mrr, mr1);
165 if (szero != 0) return;
166 mii = mi;
167 mrr.assignAll(0.0);
168 for (long n = 0; n < q; n++) mrr.ac(n, n) = DoubleAc(1.0);
169 inverse_DynArr_prot(mii, mrr, szero, serr2, 0);
170 check_econd11(szero, != 0, mcerr);
171 copy_DynArr(mrr, mr2);
172}
173
174void inverse_DynArr(const DynArr<double>& mi,
175 const DynLinArr<int>& s_var, // 1 if variable
176 DynArr<double>& mr, int& serr) {
177 mfunname(
178 "void inverse_DynArr(const DynArr<double>& mi, const "
179 "DynLinArr<long>& s_var, DynArr<double>& mr, int& serr)");
180 const DynLinArr<long>& miqel(mi.get_qel());
181 check_econd11(miqel.get_qel(), != 2, mcerr);
182 check_econd11(miqel[0], <= 0, mcerr);
183 check_econd12(miqel[0], !=, miqel[1], mcerr);
184 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
185
186 long q = s_var.get_qel();
187 long qvar = 0;
188 int s = 1; // sign that all parameters are varied
189 for (long n = 0; n < q; n++) {
190 if (s_var[n] == 0)
191 s = 0;
192 else
193 qvar++;
194 }
195 if (s == 1) {
196 inverse_DynArr(mi, mr, serr);
197 return;
198 }
199 check_econd11(qvar, <= 0, mcerr);
200 DynArr<double> mi1(qvar, qvar);
201 int nv1 = 0;
202 for (int n1 = 0; n1 < q; n1++) {
203 if (s_var[n1] != 1) continue;
204 int nv2 = 0;
205 for (int n2 = 0; n2 < q; n2++) {
206 if (s_var[n2] == 1) {
207 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
208 nv2++;
209 }
210 }
211 nv1++;
212 }
213 DynArr<double> mr1;
214 inverse_DynArr(mi1, mr1, serr);
215 mr = DynArr<double>(q, q);
216 mr.assignAll(0.0);
217 nv1 = 0;
218 for (int n1 = 0; n1 < q; n1++) {
219 if (s_var[n1] != 1) continue;
220 int nv2 = 0;
221 for (int n2 = 0; n2 < q; n2++) {
222 if (s_var[n2] == 1) {
223 mr.ac(n1, n2) = mr1.ac(nv1, nv2);
224 nv2++;
225 }
226 }
227 nv1++;
228 }
229}
230
231void inverse_DynArr_prot(const DynArr<DoubleAc>& mi,
232 const DynLinArr<int>& s_var, // 1 if variable
233 DynArr<DoubleAc>& mr, int& szero, int& serr,
234 int s_stop) {
235 mfunname(
236 "void inverse_DynArr_prot(const DynArr<DoubleAc>& mi, const "
237 "DynLinArr<long>& s_var, DynArr<DoubleAc>& mr, int& szero, int& "
238 "serr, int s_stop=1)");
239 const DynLinArr<long>& miqel(mi.get_qel());
240 check_econd11(miqel.get_qel(), != 2, mcerr);
241 check_econd11(miqel[0], <= 0, mcerr);
242 check_econd12(miqel[0], !=, miqel[1], mcerr);
243 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
244
245 long q = s_var.get_qel();
246
247 long qvar = 0;
248 int s = 1; // sign that all parameters are varied
249 for (long n = 0; n < q; n++) {
250 if (s_var[n] == 0)
251 s = 0;
252 else
253 qvar++;
254 }
255 if (s == 1) {
256 inverse_DynArr_prot(mi, mr, szero, serr, s_stop);
257 return;
258 }
259 check_econd11(qvar, <= 0, mcerr);
260 DynArr<DoubleAc> mi1(qvar, qvar);
261 int nv1 = 0;
262 for (int n1 = 0; n1 < q; n1++) {
263 if (s_var[n1] != 1) continue;
264 int nv2 = 0;
265 for (int n2 = 0; n2 < q; n2++) {
266 if (s_var[n2] == 1) {
267 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
268 nv2++;
269 }
270 }
271 nv1++;
272 }
273 DynArr<DoubleAc> mr1;
274 inverse_DynArr_prot(mi1, mr1, szero, serr, s_stop);
275 mr = DynArr<DoubleAc>(q, q);
276 mr.assignAll(0.0);
277 nv1 = 0;
278 for (int n1 = 0; n1 < q; n1++) {
279 if (s_var[n1] != 1) continue;
280 int nv2 = 0;
281 for (int n2 = 0; n2 < q; n2++) {
282 if (s_var[n2] == 1) {
283 mr.ac(n1, n2) = mr1.ac(nv1, nv2);
284 nv2++;
285 }
286 }
287 nv1++;
288 }
289}
290
291void inverse_DynArr(const DynArr<DoubleAc>& mi,
292 const DynLinArr<int>& s_var, // 1 if variable
293 DynArr<DoubleAc>& mr1, int& szero, int& serr1,
294 DynArr<DoubleAc>& mr2, int& serr2) {
295 mfunname(
296 "void inverse_DynArr(const DynArr<DoubleAc>& mi, const "
297 "DynLinArr<long>& s_var, DynArr<DoubleAc>& mr1, int& szero, int& "
298 "serr1, DynArr<DoubleAc>& mr2, int& serr2 )");
299 const DynLinArr<long>& miqel(mi.get_qel());
300 check_econd11(miqel.get_qel(), != 2, mcerr);
301 check_econd11(miqel[0], <= 0, mcerr);
302 check_econd12(miqel[0], !=, miqel[1], mcerr);
303 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
304
305 long q = s_var.get_qel();
306 long qvar = 0;
307 int s = 1; // sign that all parameters are varied
308 for (long n = 0; n < q; n++) {
309 if (s_var[n] == 0)
310 s = 0;
311 else
312 qvar++;
313 }
314 if (s == 1) {
315 inverse_DynArr(mi, mr1, szero, serr1, mr2, serr2);
316 return;
317 }
318 check_econd11(qvar, <= 0, mcerr);
319 DynArr<DoubleAc> mi1(qvar, qvar);
320 int nv1 = 0;
321 for (int n1 = 0; n1 < q; n1++) {
322 if (s_var[n1] != 1) continue;
323 int nv2 = 0;
324 for (int n2 = 0; n2 < q; n2++) {
325 if (s_var[n2] == 1) {
326 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
327 nv2++;
328 }
329 }
330 nv1++;
331 }
332 DynArr<DoubleAc> mrr1, mrr2;
333 inverse_DynArr(mi1, mrr1, szero, serr1, mrr2, serr2);
334 mr1 = DynArr<DoubleAc>(q, q);
335 mr1.assignAll(0.0);
336 if (serr1 != 1) {
337 mr2 = DynArr<DoubleAc>(q, q);
338 mr2.assignAll(0.0);
339 }
340 nv1 = 0;
341 for (int n1 = 0; n1 < q; n1++) {
342 if (s_var[n1] == 1) {
343 int nv2 = 0;
344 for (int n2 = 0; n2 < q; n2++) {
345 if (s_var[n2] == 1) {
346 mr1.ac(n1, n2) = mrr1.ac(nv1, nv2);
347 if (serr1 != 1) {
348 mr2.ac(n1, n2) = mrr2.ac(nv1, nv2);
349 }
350 nv2++;
351 }
352 }
353 nv1++;
354 }
355 }
356}
357
358DoubleAc determinant_DynArr(const DynArr<DoubleAc>& mi, long q) {
359 mfunname("DoubleAc determinant_DynArr(const DynArr<DoubleAc>& mi, long q)");
360 const DynLinArr<long>& miqel(mi.get_qel());
361 check_econd11(miqel.get_qel(), != 2, mcerr);
362 check_econd11(miqel[0], <= 0, mcerr);
363 if (q == 0) {
364 check_econd12(miqel[0], !=, miqel[1], mcerr);
365 q = miqel[0];
366 } else {
367 check_econd11(miqel[0], < q, mcerr);
368 check_econd11(miqel[1], < q, mcerr);
369 }
370 // serr=0;
371 DynArr<DoubleAc> mii(mi);
372#ifndef ALWAYS_USE_TEMPLATE_PAR_AS_FUN_PAR
373 return abstract_determinant<DynArr<DoubleAc>, DoubleAc>(mii, q);
374#else
375 DoubleAc fict;
376 return abstract_determinant(mii, q, fict);
377#endif
378}
379
380DoubleAc determinant_DynArr(const DynArr<DoubleAc>& mi,
381 const DynLinArr<int>& s_var, // 1 if variable
382 long q) {
383 mfunname("DoubleAc determinant_DynArr(...)");
384 const DynLinArr<long>& miqel(mi.get_qel());
385 check_econd11(miqel.get_qel(), != 2, mcerr);
386 check_econd11(miqel[0], <= 0, mcerr);
387 // check_econd12(miqel[0] , != , miqel[1] , mcerr);
388 if (q == 0) {
389 check_econd12(miqel[0], !=, miqel[1], mcerr);
390 q = miqel[0];
391 } else {
392 check_econd11(miqel[0], < q, mcerr);
393 check_econd11(miqel[1], < q, mcerr);
394 }
395 check_econd12(q, >, s_var.get_qel(), mcerr);
396 long miq = std::min(miqel[0], miqel[1]);
397 long qvar = 0;
398 int s = 1; // sign that all parameters are varied
399 for (long n = 0; n < s_var.get_qel(); n++) {
400 if (s_var[n] == 0) {
401 s = 0;
402 } else {
403 qvar++;
404 if (qvar == q) break;
405 }
406 }
407 if (s == 1) return determinant_DynArr(mi, q);
408
409 check_econd11(qvar, <= 0, mcerr);
410 check_econd11(qvar, < q, mcerr);
411 DynArr<DoubleAc> mi1(q, q);
412 int nv1 = 0;
413 for (int n1 = 0; n1 < miq; n1++) {
414 if (s_var[n1] != 1) continue;
415 int nv2 = 0;
416 for (int n2 = 0; n2 < miq; n2++) {
417 if (s_var[n2] == 1) {
418 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
419 nv2++;
420 if (nv2 >= q) break;
421 }
422 }
423 nv1++;
424 if (nv1 >= q) break;
425 }
426 return determinant_DynArr(mi1, q);
427}
428}
check_econd11
#define check_econd11(a, signb, stream)
Definition: FunNameStack.h:155
check_econd12
#define check_econd12(a, sign, b, stream)
Definition: FunNameStack.h:163
mfunname
#define mfunname(string)
Definition: FunNameStack.h:45
abs_inverse.h
Heed::DoubleAc::get
double get(void) const
Definition: DoubleAc.h:76
Heed::DoubleAc::left_limit
double left_limit(void) const
Definition: DoubleAc.h:79
Heed::DynArr::ac
T & ac(long i)
Definition: AbsArr.h:1717
Heed::DynArr::assignAll
void assignAll(const T &val)
Definition: AbsArr.h:2479
Heed::DynArr::get_qel
const DynLinArr< long > & get_qel(void) const
Definition: AbsArr.h:2169
Heed::DynLinArr::get_qel
long get_qel(void) const
Definition: AbsArr.h:283
Heed
Definition: BGMesh.cpp:5
Heed::determinant_DynArr
DoubleAc determinant_DynArr(const DynArr< DoubleAc > &mi, long q)
Definition: inverse.cpp:358
Heed::abstract_determinant
X abstract_determinant(M &mi, long q, X)
Definition: abs_inverse.h:32
Heed::inverse_DynArr_prot
void inverse_DynArr_prot(const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
Definition: inverse.cpp:17
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::inverse_DynArr
void inverse_DynArr(const DynArr< double > &mi, DynArr< double > &mr, int &serr)
Definition: inverse.cpp:99
Heed::copy_DynArr
void copy_DynArr(const DynArr< T > &s, DynArr< X > &d)
Definition: AbsArr.h:2494
mcerr
#define mcerr
Definition: prstream.h:128