Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Functions
inverse.cpp File Reference
#include "wcpplib/matrix/abs_inverse.h"
#include "wcpplib/matrix/inverse.h"

Go to the source code of this file.

Functions

void inverse_DynArr_prot (const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
 
void inverse_DynArr (const DynArr< double > &mi, DynArr< double > &mr, int &serr)
 
void inverse_DynArr (const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr1, int &szero, int &serr1, DynArr< DoubleAc > &mr2, int &serr2)
 
void inverse_DynArr (const DynArr< double > &mi, const DynLinArr< int > &s_var, DynArr< double > &mr, int &serr)
 
void inverse_DynArr_prot (const DynArr< DoubleAc > &mi, const DynLinArr< int > &s_var, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
 
void inverse_DynArr (const DynArr< DoubleAc > &mi, const DynLinArr< int > &s_var, DynArr< DoubleAc > &mr1, int &szero, int &serr1, DynArr< DoubleAc > &mr2, int &serr2)
 
DoubleAc determinant_DynArr (const DynArr< DoubleAc > &mi, long q)
 
DoubleAc determinant_DynArr (const DynArr< DoubleAc > &mi, const DynLinArr< int > &s_var, long q)
 

Function Documentation

◆ determinant_DynArr() [1/2]

DoubleAc determinant_DynArr ( const DynArr< DoubleAc > &  mi,
const DynLinArr< int > &  s_var,
long  q 
)

Definition at line 426 of file inverse.cpp.

428 {
429 mfunname("DoubleAc determinant_DynArr(...)");
430 const DynLinArr<long>& miqel(mi.get_qel());
431 check_econd11(miqel.get_qel(), != 2, mcerr);
432 check_econd11(miqel[0], <= 0, mcerr);
433 //check_econd12(miqel[0] , != , miqel[1] , mcerr);
434 if (q == 0) {
435 check_econd12(miqel[0], !=, miqel[1], mcerr);
436 q = miqel[0];
437 } else {
438 check_econd11(miqel[0], < q, mcerr);
439 check_econd11(miqel[1], < q, mcerr);
440 }
441 check_econd12(q, >, s_var.get_qel(), mcerr);
442 //DynArr<DoubleAc> mi(der[1]);
443 long miq = find_min(miqel[0], miqel[1]);
444 long n;
445 long qvar = 0;
446 int s = 1; // sign that all parameters are varied
447 for (n = 0; n < s_var.get_qel(); n++) {
448 if (s_var[n] == 0)
449 s = 0;
450 else {
451 qvar++;
452 if (qvar == q) break;
453 }
454 }
455 if (s == 1) {
456 return determinant_DynArr(mi, q);
457 } else {
458 check_econd11(qvar, <= 0, mcerr);
459 check_econd11(qvar, < q, mcerr);
460 DynArr<DoubleAc> mi1(q, q);
461 //HepMatrix tempmat(qvar, qvar, 0);
462 int n1, n2, nv1 = 0;
463 int nv2; //DynLinArr<long>& indexes(2);
464 for (n1 = 0; n1 < miq; n1++) {
465 if (s_var[n1] == 1) {
466 nv2 = 0;
467 for (n2 = 0; n2 < miq; n2++) {
468 if (s_var[n2] == 1) {
469 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
470 nv2++;
471 if (nv2 >= q) break;
472 }
473 }
474 nv1++;
475 if (nv1 >= q) break;
476 }
477 }
478 return determinant_DynArr(mi1, q);
479 }
480}
find_min
DoubleAc find_min(const DoubleAc &a, const DoubleAc &b)
Definition: DoubleAc.h:627
check_econd11
#define check_econd11(a, signb, stream)
Definition: FunNameStack.h:366
check_econd12
#define check_econd12(a, sign, b, stream)
Definition: FunNameStack.h:380
mfunname
#define mfunname(string)
Definition: FunNameStack.h:67
DynArr::ac
T & ac(long i)
Definition: AbsArr.h:2057
DynArr::get_qel
const DynLinArr< long > & get_qel(void) const
Definition: AbsArr.h:2548
DynLinArr::get_qel
long get_qel(void) const
Definition: AbsArr.h:420
determinant_DynArr
DoubleAc determinant_DynArr(const DynArr< DoubleAc > &mi, long q)
Definition: inverse.cpp:404
mcerr
#define mcerr
Definition: prstream.h:135

◆ determinant_DynArr() [2/2]

DoubleAc determinant_DynArr ( const DynArr< DoubleAc > &  mi,
long  q 
)

Definition at line 404 of file inverse.cpp.

404 {
405 mfunname("DoubleAc determinant_DynArr(const DynArr<DoubleAc>& mi, long q)");
406 const DynLinArr<long>& miqel(mi.get_qel());
407 check_econd11(miqel.get_qel(), != 2, mcerr);
408 check_econd11(miqel[0], <= 0, mcerr);
409 if (q == 0) {
410 check_econd12(miqel[0], !=, miqel[1], mcerr);
411 q = miqel[0];
412 } else {
413 check_econd11(miqel[0], < q, mcerr);
414 check_econd11(miqel[1], < q, mcerr);
415 }
416 //serr=0;
417 DynArr<DoubleAc> mii(mi);
418#ifndef ALWAYS_USE_TEMPLATE_PAR_AS_FUN_PAR
419 return abstract_determinant<DynArr<DoubleAc>, DoubleAc>(mii, q);
420#else
421 DoubleAc fict;
422 return abstract_determinant(mii, q, fict);
423#endif
424}
abstract_determinant
X abstract_determinant(M &mi, long q, X)
Definition: abs_inverse.h:104

Referenced by determinant_DynArr().

◆ inverse_DynArr() [1/4]

void inverse_DynArr ( const DynArr< double > &  mi,
const DynLinArr< int > &  s_var,
DynArr< double > &  mr,
int &  serr 
)

Definition at line 204 of file inverse.cpp.

206 {
207 mfunname("void inverse_DynArr(const DynArr<double>& mi, const "
208 "DynLinArr<long>& s_var, DynArr<double>& mr, int& serr)");
209 const DynLinArr<long>& miqel(mi.get_qel());
210 check_econd11(miqel.get_qel(), != 2, mcerr);
211 check_econd11(miqel[0], <= 0, mcerr);
212 check_econd12(miqel[0], !=, miqel[1], mcerr);
213 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
214
215 long q = s_var.get_qel();
216 //DynArr<DoubleAc> mi(der[1]);
217
218 long n;
219 long qvar = 0;
220 int s = 1; // sign that all parameters are varied
221 for (n = 0; n < q; n++) {
222 if (s_var[n] == 0)
223 s = 0;
224 else
225 qvar++;
226 }
227 if (s == 1) {
228 inverse_DynArr(mi, mr, serr);
229 } else {
230 check_econd11(qvar, <= 0, mcerr);
231 DynArr<double> mi1(qvar, qvar);
232 //HepMatrix tempmat(qvar, qvar, 0);
233 int n1, n2, nv1 = 0;
234 int nv2; //DynLinArr<long>& indexes(2);
235 for (n1 = 0; n1 < q; n1++) {
236 if (s_var[n1] == 1) {
237 nv2 = 0;
238 for (n2 = 0; n2 < q; n2++) {
239 if (s_var[n2] == 1) {
240 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
241 nv2++;
242 }
243 }
244 nv1++;
245 }
246 }
247 DynArr<double> mr1;
248 inverse_DynArr(mi1, mr1, serr);
249 mr = DynArr<double>(q, q);
250 mr.assignAll(0.0);
251 nv1 = 0;
252 for (n1 = 0; n1 < q; n1++) {
253 if (s_var[n1] == 1) {
254 nv2 = 0;
255 for (n2 = 0; n2 < q; n2++) {
256 if (s_var[n2] == 1) {
257 mr.ac(n1, n2) = mr1.ac(nv1, nv2);
258 nv2++;
259 }
260 }
261 nv1++;
262 }
263 }
264 }
265}
DynArr::assignAll
void assignAll(const T &val)
Definition: AbsArr.h:2880
inverse_DynArr
void inverse_DynArr(const DynArr< double > &mi, DynArr< double > &mr, int &serr)
Definition: inverse.cpp:113

◆ inverse_DynArr() [2/4]

void inverse_DynArr ( const DynArr< double > &  mi,
DynArr< double > &  mr,
int &  serr 
)

Definition at line 113 of file inverse.cpp.

113 {
114 mfunname("void inverse_DynArr(const DynArr<double>& mi, DynArr<double>& mr, "
115 "int& serr)");
116 const DynLinArr<long>& miqel(mi.get_qel());
117 check_econd11(miqel.get_qel(), != 2, mcerr);
118 check_econd11(miqel[0], <= 0, mcerr);
119 check_econd12(miqel[0], !=, miqel[1], mcerr);
120 serr = 0;
121 long q = miqel[0];
122 mr = DynArr<double>(q, q);
123 if (q == 1) {
124 if (mi.ac(0, 0) == 0.0) {
125 serr = 1;
126 return;
127 }
128 mr.ac(0, 0) = 1.0 / mi.ac(0, 0);
129 return;
130 }
131 DynArr<DoubleAc> mii;
132 DynArr<DoubleAc> mrr(q, q);
133 copy_DynArr(mi, mii);
134 mrr.assignAll(0.0);
135 long n;
136 for (n = 0; n < miqel[0]; n++)
137 mrr.ac(n, n) = 1.0;
138 int szero;
139 inverse_DynArr_prot(mii, mrr, szero, serr);
140 copy_DynArr(mrr, mr);
141}
copy_DynArr
void copy_DynArr(const DynArr< T > &s, DynArr< X > &d)
Definition: AbsArr.h:2894
inverse_DynArr_prot
void inverse_DynArr_prot(const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
Definition: inverse.cpp:15

Referenced by inverse_DynArr().

◆ inverse_DynArr() [3/4]

void inverse_DynArr ( const DynArr< DoubleAc > &  mi,
const DynLinArr< int > &  s_var,
DynArr< DoubleAc > &  mr1,
int &  szero,
int &  serr1,
DynArr< DoubleAc > &  mr2,
int &  serr2 
)

Definition at line 332 of file inverse.cpp.

335 {
336 mfunname("void inverse_DynArr(const DynArr<DoubleAc>& mi, const "
337 "DynLinArr<long>& s_var, DynArr<DoubleAc>& mr1, int& szero, int& "
338 "serr1, DynArr<DoubleAc>& mr2, int& serr2 )");
339 const DynLinArr<long>& miqel(mi.get_qel());
340 check_econd11(miqel.get_qel(), != 2, mcerr);
341 check_econd11(miqel[0], <= 0, mcerr);
342 check_econd12(miqel[0], !=, miqel[1], mcerr);
343 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
344
345 long q = s_var.get_qel();
346 //DynArr<DoubleAc> mi(der[1]);
347
348 long n;
349 long qvar = 0;
350 int s = 1; // sign that all parameters are varied
351 for (n = 0; n < q; n++) {
352 if (s_var[n] == 0)
353 s = 0;
354 else
355 qvar++;
356 }
357 if (s == 1) {
358 inverse_DynArr(mi, mr1, szero, serr1, mr2, serr2);
359 } else {
360 check_econd11(qvar, <= 0, mcerr);
361 DynArr<DoubleAc> mi1(qvar, qvar);
362 //HepMatrix tempmat(qvar, qvar, 0);
363 int n1, n2, nv1 = 0;
364 int nv2; //DynLinArr<long>& indexes(2);
365 for (n1 = 0; n1 < q; n1++) {
366 if (s_var[n1] == 1) {
367 nv2 = 0;
368 for (n2 = 0; n2 < q; n2++) {
369 if (s_var[n2] == 1) {
370 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
371 nv2++;
372 }
373 }
374 nv1++;
375 }
376 }
377 DynArr<DoubleAc> mrr1, mrr2;
378 inverse_DynArr(mi1, mrr1, szero, serr1, mrr2, serr2);
379 mr1 = DynArr<DoubleAc>(q, q);
380 mr1.assignAll(0.0);
381 if (serr1 != 1) {
382 mr2 = DynArr<DoubleAc>(q, q);
383 mr2.assignAll(0.0);
384 }
385 nv1 = 0;
386 for (n1 = 0; n1 < q; n1++) {
387 if (s_var[n1] == 1) {
388 nv2 = 0;
389 for (n2 = 0; n2 < q; n2++) {
390 if (s_var[n2] == 1) {
391 mr1.ac(n1, n2) = mrr1.ac(nv1, nv2);
392 if (serr1 != 1) {
393 mr2.ac(n1, n2) = mrr2.ac(nv1, nv2);
394 }
395 nv2++;
396 }
397 }
398 nv1++;
399 }
400 }
401 }
402}

◆ inverse_DynArr() [4/4]

void inverse_DynArr ( const DynArr< DoubleAc > &  mi,
DynArr< DoubleAc > &  mr1,
int &  szero,
int &  serr1,
DynArr< DoubleAc > &  mr2,
int &  serr2 
)

Definition at line 143 of file inverse.cpp.

144 {
145 mfunname("void inverse_DynArr(const DynArr<DoubleAc>& mi, DynArr<DoubleAc>& "
146 "mr1, int& szero, int& serr1, DynArr<DoubleAc>& mr2, int& serr2)");
147 const DynLinArr<long>& miqel(mi.get_qel());
148 check_econd11(miqel.get_qel(), != 2, mcerr);
149 check_econd11(miqel[0], <= 0, mcerr);
150 check_econd12(miqel[0], !=, miqel[1], mcerr);
151 serr1 = 0;
152 serr2 = 0;
153 long q = miqel[0];
154 //mr = DynArr<DoubleAc>(miqel[0], miqel[0]);
155 if (q == 1) {
156 if (mi.ac(0, 0).get() == 0.0) {
157 serr1 = 1;
158 return;
159 }
160 mr1 = DynArr<DoubleAc>(q, q);
161 mr2 = DynArr<DoubleAc>(q, q);
162 mr1.ac(0, 0) = 1.0 / mi.ac(0, 0).get();
163 mr2.ac(0, 0) = DoubleAc(1.0) / mi.ac(0, 0);
164 //Iprintdan(mcout, mr2.ac(0,0));
165 if (fabs(mr2.ac(0, 0)).left_limit() == 0.0) serr2 = 1;
166 return;
167 }
168 //DynArr<DoubleAc> mii(mi);
169 //mr.assignAll(0.0);
170 //long n;
171 //for( n=0; n<miqel[0]; n++)
172 // mr.ac(n,n)=DoubleAc(1.0);
173 //copy_DynArr(mi, mii);
174
175 DynArr<DoubleAc> mii(q, q);
176 long n1, n2;
177 for (n1 = 0; n1 < q; n1++) {
178 for (n2 = 0; n2 < q; n2++) {
179 mii.ac(n1, n2) = double(mi.ac(n1, n2));
180 }
181 }
182 DynArr<DoubleAc> mrr(q, q);
183 mrr.assignAll(0.0);
184 long n;
185 for (n = 0; n < q; n++)
186 mrr.ac(n, n) = 1.0;
187 //int szero;
188 inverse_DynArr_prot(mii, mrr, szero, serr1, 0);
189 copy_DynArr(mrr, mr1);
190 if (szero != 0) return;
191 //if(serr1 != 0) return;
192 mii = mi;
193 mrr.assignAll(0.0);
194 for (n = 0; n < q; n++)
195 mrr.ac(n, n) = DoubleAc(1.0);
196 inverse_DynArr_prot(mii, mrr, szero, serr2, 0);
197 check_econd11(szero, != 0, mcerr);
198 copy_DynArr(mrr, mr2);
199
200 //abstract_inverse<DynArr<DoubleAc> , DoubleAc>
201 // (mii, mr, miqel[0], serr);
202}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
DoubleAc::left_limit
double left_limit(void) const
Definition: DoubleAc.h:79

◆ inverse_DynArr_prot() [1/2]

void inverse_DynArr_prot ( const DynArr< DoubleAc > &  mi,
const DynLinArr< int > &  s_var,
DynArr< DoubleAc > &  mr,
int &  szero,
int &  serr,
int  s_stop 
)

Definition at line 267 of file inverse.cpp.

270 {
271 mfunname("void inverse_DynArr_prot(const DynArr<DoubleAc>& mi, const "
272 "DynLinArr<long>& s_var, DynArr<DoubleAc>& mr, int& szero, int& "
273 "serr, int s_stop=1)");
274 const DynLinArr<long>& miqel(mi.get_qel());
275 check_econd11(miqel.get_qel(), != 2, mcerr);
276 check_econd11(miqel[0], <= 0, mcerr);
277 check_econd12(miqel[0], !=, miqel[1], mcerr);
278 check_econd12(s_var.get_qel(), !=, miqel[0], mcerr);
279
280 long q = s_var.get_qel();
281 //DynArr<DoubleAc> mi(der[1]);
282
283 long n;
284 long qvar = 0;
285 int s = 1; // sign that all parameters are varied
286 for (n = 0; n < q; n++) {
287 if (s_var[n] == 0)
288 s = 0;
289 else
290 qvar++;
291 }
292 if (s == 1) {
293 inverse_DynArr_prot(mi, mr, szero, serr, s_stop);
294 } else {
295 check_econd11(qvar, <= 0, mcerr);
296 DynArr<DoubleAc> mi1(qvar, qvar);
297 //HepMatrix tempmat(qvar, qvar, 0);
298 int n1, n2, nv1 = 0;
299 int nv2; //DynLinArr<long>& indexes(2);
300 for (n1 = 0; n1 < q; n1++) {
301 if (s_var[n1] == 1) {
302 nv2 = 0;
303 for (n2 = 0; n2 < q; n2++) {
304 if (s_var[n2] == 1) {
305 mi1.ac(nv1, nv2) = mi.ac(n1, n2);
306 nv2++;
307 }
308 }
309 nv1++;
310 }
311 }
312 DynArr<DoubleAc> mr1;
313 inverse_DynArr_prot(mi1, mr1, szero, serr, s_stop);
314 mr = DynArr<DoubleAc>(q, q);
315 mr.assignAll(0.0);
316 nv1 = 0;
317 for (n1 = 0; n1 < q; n1++) {
318 if (s_var[n1] == 1) {
319 nv2 = 0;
320 for (n2 = 0; n2 < q; n2++) {
321 if (s_var[n2] == 1) {
322 mr.ac(n1, n2) = mr1.ac(nv1, nv2);
323 nv2++;
324 }
325 }
326 nv1++;
327 }
328 }
329 }
330}

◆ inverse_DynArr_prot() [2/2]

void inverse_DynArr_prot ( const DynArr< DoubleAc > &  mi,
DynArr< DoubleAc > &  mr,
int &  szero,
int &  serr,
int  s_stop 
)

Definition at line 15 of file inverse.cpp.

16 {
17 mfunname("void inverse_DynArr_prot(const DynArr<DoubleAc>& mi, "
18 "DynArr<DoubleAc>& mr, int& s_zero, int& serr, int s_stop)");
19 //mcout<<"inverse_DynArr_prot:\n";
20 //Iprintda_DoubleAc(mcout, mi, 3);
21 const DynLinArr<long>& miqel(mi.get_qel());
22 check_econd11(miqel.get_qel(), != 2, mcerr);
23 check_econd11(miqel[0], <= 0, mcerr);
24 check_econd12(miqel[0], !=, miqel[1], mcerr);
25 serr = 0;
26 szero = 0;
27 long q = miqel[0];
28 mr = DynArr<DoubleAc>(q, q);
29 if (q == 1) {
30 if (mi.ac(0, 0).get() == 0.0) {
31 szero = 1;
32 serr = 1;
33 return;
34 }
35 mr.ac(0, 0) = 1.0 / mi.ac(0, 0);
36 if (fabs(mr.ac(0, 0)).left_limit() == 0) {
37 serr = 1;
38 }
39 return;
40 }
41 DynArr<DoubleAc> mii(mi);
42 //Iprintda_DoubleAc(mcout, mii, 3);
43 mr.assignAll(0.0);
44 long n;
45 for (n = 0; n < q; n++)
46 mr.ac(n, n) = DoubleAc(1.0);
47 //Iprintda_DoubleAc(mcout, mr, 3);
48
49 long nr, nr1, nc;
50 for (nr = 0; nr < q; nr++) {
51 //Iprintn(mcout, nr);
52 long nmax = 0;
53 DoubleAc d(0.0);
54 for (nr1 = nr; nr1 < q; nr1++) {
55 if (fabs(mii.ac(nr1, nr)) > d) {
56 d = fabs(mii.ac(nr1, nr));
57 nmax = nr1;
58 }
59 }
60 //Iprintdan(mcout, d);
61 //mcout<<"d="<<d<<'\n';
62 //mcout<<"nmax="<<nmax<<'\n';
63 if (d.get() == 0.0) {
64 szero = 1;
65 serr = 1;
66 return;
67 }
68 if (d.left_limit() == 0) {
69 serr = 1;
70 if (s_stop == 1) return;
71 }
72 if (nmax > nr) {
73 for (nc = nr; nc < q; nc++) {
74 DoubleAc t(mii.ac(nr, nc));
75 mii.ac(nr, nc) = mii.ac(nmax, nc);
76 mii.ac(nmax, nc) = t;
77 }
78 for (nc = 0; nc < q; nc++) {
79 DoubleAc t(mr.ac(nr, nc));
80 mr.ac(nr, nc) = mr.ac(nmax, nc);
81 mr.ac(nmax, nc) = t;
82 }
83 //long tl=order[nr];
84 //order[nr] = order[nmax];
85 //order[nmax] = tl;
86 }
87 DoubleAc t = mii.ac(nr, nr);
88 for (nr1 = 0; nr1 < q; nr1++) {
89 if (nr1 != nr) {
90 DoubleAc k(mii.ac(nr1, nr) / t);
91 //mcout<<"nr1="<<nr1<<" nr="<<nr<<'\n';
92 //mcout<<"k="<<k<<'\n';
93 for (nc = nr; nc < q; nc++) {
94 mii.ac(nr1, nc) -= k * mii.ac(nr, nc);
95 }
96 for (nc = 0; nc < q; nc++) {
97 mr.ac(nr1, nc) -= k * mr.ac(nr, nc);
98 }
99 }
100 }
101 for (nc = nr; nc < q; nc++) {
102 mii.ac(nr, nc) /= t;
103 }
104 for (nc = 0; nc < q; nc++) {
105 mr.ac(nr, nc) /= t;
106 }
107 //Iprintda_DoubleAc(mcout, mii, 3);
108 //Iprintda_DoubleAc(mcout, mr, 3);
109 }
110 //Iprintda_DoubleAc(mcout, mr, 3);
111}

Referenced by inverse_DynArr(), inverse_DynArr_prot(), and Parabol::Parabol().