Garfield++ v1r0
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.
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
16 int& szero, int& serr, int s_stop) {
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}
112
113void inverse_DynArr(const DynArr<double>& mi, DynArr<double>& mr, int& serr) {
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 }
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}
142
144 int& szero, int& serr1, DynArr<DoubleAc>& mr2, int& serr2) {
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}
203
205 const DynLinArr<int>& s_var, // 1 if variable
206 DynArr<double>& mr, int& serr) {
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}
266
268 const DynLinArr<int>& s_var, // 1 if variable
269 DynArr<DoubleAc>& mr, int& szero, int& serr,
270 int s_stop) {
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 }
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}
331
333 const DynLinArr<int>& s_var, // 1 if variable
334 DynArr<DoubleAc>& mr1, int& szero, int& serr1,
335 DynArr<DoubleAc>& mr2, int& serr2) {
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}
403
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}
425
427 const DynLinArr<int>& s_var, // 1 if variable
428 long q) {
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}
481
482/*
483DoubleAc determinant_DynArr_prot(const DynArr<DoubleAc>& mi,
484 long q, // dimension of minor
485 int& szero, // sign that the calculations are
486 // terminated owing to attempt to divide by 0.
487 // The final determinant is not correct.
488 int& serr, // sign that the interval precision
489 // is broken
490 // (but the final matrix may be provided if
491 // szero=0)
492 int s_stop=1 // directive to stop if
493 // the interval precision is broken
494 )
495{
496 mfunname("DoubleAc determinant_prot(...)");
497 check_econd11(mi.get_qdim() , != 2 , mcerr);
498 check_econd12(mi.get_qel()[0] , < , q , mcerr);
499 check_econd12(mi.get_qel()[1] , < , q , mcerr);
500 check_econd11(q , < 1 , mcerr);
501 serr=0;
502 szero=0;
503 if(q == 1)
504 {
505 return mi.ac(0,0);
506 }
507 else if(q == 2)
508 {
509 return mi.ac(0,0)*mi.ac(1,1) - mi.ac(0,1)*mi.ac(1,0);
510 }
511 else if(q == 3)
512 {
513 return
514 mi.ac(0,0)*mi.ac(1,1)*mi.ac(2,2) +
515 mi.ac(0,2)*mi.ac(1,0)*mi.ac(2,1) +
516 mi.ac(0,1)*mi.ac(1,2)*mi.ac(2,0) -
517 mi.ac(0,2)*mi.ac(1,1)*mi.ac(2,0) -
518 mi.ac(0,0)*mi.ac(1,2)*mi.ac(2,1) -
519 mi.ac(0,1)*mi.ac(1,0)*mi.ac(2,2);
520 }
521 else
522 {
523
524 long nr, nr1, nc;
525 DoubleAc koef=1;
526 DynArr<DoubleAc> mii(mi);
527 for(nr=0; nr<q; nr++)
528 {
529 long nmax=0;
530 DoubleAc d=0;
531 for(nr1=nr; nr1<q; nr1++)
532 {
533 if(fabs(mii.ac(nr1,nr)) > d)
534 {
535 d = fabs(mii.ac(nr1,nr));
536 nmax = nr1;
537 }
538 }
539 //mcout<<"d="<<d<<'\n';
540 //mcout<<"nmax="<<nmax<<'\n';
541 if(d.get() == 0.0)
542 {
543 szero=1;
544 serr = 1;
545 return koef;
546 }
547 if(d.left_limit() == 0)
548 {
549 serr = 1;
550 if(s_stop == 1)
551 return koef;
552 }
553 if(nmax > nr)
554 {
555 for(nc=nr; nc<q; nc++)
556 {
557 DoubleAc t(mii.ac(nr,nc));
558 mii.ac(nr,nc) = mii.ac(nmax,nc);
559 mii.ac(nmax,nc) = t;
560 }
561 koef *= -1; // trancposition of rows: determinant changes sign
562 }
563 DoubleAc t=mii.ac(nr,nr);
564 for(nr1=nr+1; nr1<q; nr1++)
565 {
566 DoubleAc k(mii.ac(nr1,nr)/t);
567 //mcout<<"nr1="<<nr1<<" nr="<<nr<<'\n';
568 //mcout<<"k="<<k<<'\n';
569 for(nc=nr; nc<q; nc++)
570 {
571 mii.ac(nr1,nc) -= k * mii.ac(nr,nc);
572 } // add elements of another row: the main value of
573 // determinant is not affected (proven in linear algebra)
574 // But the resolution gets worser.
575 }
576 //for(nc=nr; nc<q; nc++)
577 //{
578 // mii.ac(nr,nc) /= t;
579 //}
580 koef *= t;
581 }
582 return koef;
583 }
584}
585*/
void copy_DynArr(const DynArr< T > &s, DynArr< X > &d)
Definition: AbsArr.h:2894
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
DoubleAc find_min(const DoubleAc &a, const DoubleAc &b)
Definition: DoubleAc.h:627
#define check_econd11(a, signb, stream)
Definition: FunNameStack.h:366
#define check_econd12(a, sign, b, stream)
Definition: FunNameStack.h:380
#define mfunname(string)
Definition: FunNameStack.h:67
X abstract_determinant(M &mi, long q, X)
Definition: abs_inverse.h:104
double left_limit(void) const
Definition: DoubleAc.h:79
double get(void) const
Definition: DoubleAc.h:76
void assignAll(const T &val)
Definition: AbsArr.h:2880
T & ac(long i)
Definition: AbsArr.h:2057
const DynLinArr< long > & get_qel(void) const
Definition: AbsArr.h:2548
long get_qel(void) const
Definition: AbsArr.h:420
void inverse_DynArr_prot(const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
Definition: inverse.cpp:15
void inverse_DynArr(const DynArr< double > &mi, DynArr< double > &mr, int &serr)
Definition: inverse.cpp:113
DoubleAc determinant_DynArr(const DynArr< DoubleAc > &mi, long q)
Definition: inverse.cpp:404
#define mcerr
Definition: prstream.h:135