FEDRA emulsion software from the OPERA Collaboration
Namespaces | Functions
VtSqMatrix.C File Reference
#include "Riostream.h"
#include <string.h>
#include "vt++/VtSqMatrix.hh"
#include "vt++/VtVector.hh"
Include dependency graph for VtSqMatrix.C:

Namespaces

namespace  MATRIX
 

Functions

void Dfact1 (int *n, double *a, int *idim, double *ir, int *ifail, double *det, int *jfail)
 
void Dfactir1 (int *n, double *a, int *idim, int *ir, int *ifail, double *det, int *jfail)
 
void Dfinv1 (int *n, double *a, int *idim, int *ir, int *ifail)
 
void Dinv1 (int *n, double *a, int *idim, double *ir, int *ifail)
 

Function Documentation

◆ Dfact1()

void Dfact1 ( int *  n,
double *  a,
int *  idim,
double *  ir,
int *  ifail,
double *  det,
int *  jfail 
)
225{
226 /* Local variables */
227 static int nxch, i, j, k, l;
228 static double p, q, tf;
229
230 *ifail = 0;
231 *jfail = 0;
232
233 if (*idim < *n || *n <= 0) {
234 *ifail =-1;
235 return;
236 }
237
238 /* Parameter adjustments */
239 a -= *idim + 1;
240
241 /* Function Body */
242
243 // fact.inc
244
245 nxch = 0;
246 *det = 1.;
247 for (j = 1; j <= *n; ++j) {
248 const int ji = j * (*idim);
249 const int jj = j + ji;
250
251 k = j;
252 p = fabs(a[jj]);
253
254 if (j != *n) {
255 for (i = j + 1; i <= *n; ++i) {
256 q = fabs(a[i + ji]);
257 if (q > p) {
258 k = i;
259 p = q;
260 }
261 } // for i
262 if (k != j) {
263 for (l = 1; l <= *n; ++l) {
264 const int li = l*(*idim);
265 const int jli = j + li;
266 const int kli = k + li;
267 tf = a[jli];
268 a[jli] = a[kli];
269 a[kli] = tf;
270 } // for l
271 ++nxch;
272 } // if k != j
273 } // if j!=n
274
275 if (p <= 0.) {
276 *det = 0;
277 *ifail = -1;
278 return;
279 }
280
281 *det *= a[jj];
282
283// t = fabs(*det);
284// if (t < 1e-19 || t > 1e19) {
285// *det = 0;
286// *ifail = -1;
287// return;
288// }
289
290 a[jj] = 1. / a[jj];
291 if (j == *n) {
292 continue;
293 }
294
295 const int jm1 = j - 1;
296 const int jpi = (j + 1) * (*idim);
297 const int jjpi = j + jpi;
298
299 for (k = j + 1; k <= *n; ++k) {
300 const int ki = k * (*idim);
301 const int jki = j + ki;
302 const int kji = k + jpi;
303 if (j != 1) {
304 for (i = 1; i <= jm1; ++i) {
305 const int ii = i * (*idim);
306 a[jki] -= a[i + ki] * a[j + ii];
307 a[kji] -= a[i + jpi] * a[k + ii];
308 } // for i
309 }
310 a[jki] *= a[jj];
311 a[kji] -= a[jjpi] * a[k + ji];
312 } // for k
313 } // for j
314
315 if (nxch % 2 != 0) {
316 *det = -(*det);
317 }
318 return;
319} // end of Dfact1
fabs
Expr< UnaryOp< Fabs< T >, Expr< A, T, D >, T >, T, D > fabs(const Expr< A, T, D > &rhs)
Definition: UnaryOperators.hh:96
a
void a()
Definition: check_aligned.C:59
q
q
Definition: testBGReduction_AllMethods.C:55
p
p
Definition: testBGReduction_AllMethods.C:8

◆ Dfactir1()

void Dfactir1 ( int *  n,
double *  a,
int *  idim,
int *  ir,
int *  ifail,
double *  det,
int *  jfail 
)
32{
33 /* Local variables */
34 static int nxch, i, j, k, l;
35 static double p, q, tf;
36
37 *ifail = 0;
38 *jfail = 0;
39
40 if (*idim < *n || *n <= 0) {
41 *ifail = -1;
42 return;
43 }
44
45
46 /* Parameter adjustments */
47 a -= *idim + 1;
48 --ir;
49
50 /* Function Body */
51
52 // fact.inc
53 nxch = 0;
54 *det = 1.;
55 for (j = 1; j <= *n; ++j) {
56 const int ji = j * (*idim);
57 const int jj = j + ji;
58
59 k = j;
60 p = fabs(a[jj]);
61
62 if (j != *n) {
63 for (i = j + 1; i <= *n; ++i) {
64 q = fabs(a[i + ji]);
65 if (q > p) {
66 k = i;
67 p = q;
68 }
69 } // for i
70
71 if (k != j) {
72 for (l = 1; l <= *n; ++l) {
73 const int li = l*(*idim);
74 const int jli = j + li;
75 const int kli = k + li;
76 tf = a[jli];
77 a[jli] = a[kli];
78 a[kli] = tf;
79 } // for l
80 ++nxch;
81 ir[nxch] = (j << 12) + k;
82 } // if k != j
83 } // if j!=n
84
85 if (p <= 0.) {
86 *det = 0;
87 *ifail = -1;
88 return;
89 }
90
91 *det *= a[jj];
92// t = fabs(*det);
93// if (t < 1e-19 || t > 1e19) {
94// *det = 0;
95// *ifail = -1;
96// return;
97// }
98
99 a[jj] = 1. / a[jj];
100 if (j == *n) {
101 continue;
102 }
103
104 const int jm1 = j - 1;
105 const int jpi = (j + 1) * (*idim);
106 const int jjpi = j + jpi;
107
108 for (k = j + 1; k <= *n; ++k) {
109 const int ki = k * (*idim);
110 const int jki = j + ki;
111 const int kji = k + jpi;
112 if (j != 1) {
113 for (i = 1; i <= jm1; ++i) {
114 const int ii = i * (*idim);
115 a[jki] -= a[i + ki] * a[j + ii];
116 a[kji] -= a[i + jpi] * a[k + ii];
117 } // for i
118 }
119 a[jki] *= a[jj];
120 a[kji] -= a[jjpi] * a[k + ji];
121 } // for k
122 } // for j
123
124 if (nxch % 2 != 0) {
125 *det = -(*det);
126 }
127 ir[*n] = nxch;
128 return;
129} // end of Dfactir1

◆ Dfinv1()

void Dfinv1 ( int *  n,
double *  a,
int *  idim,
int *  ir,
int *  ifail 
)
132{
133 /* Local variables */
134 static int nxch, i, j, k, m, ij;
135 static int im2, nm1, nmi;
136 static double s31, s34, ti;
137
138 *ifail = 0;
139
140 if (*idim < *n || *n <= 0 || *n==1) {
141 *ifail = -1;
142 return;
143 }
144
145
146 /* Parameter adjustments */
147 a -= *idim + 1;
148 --ir;
149
150 /* Function Body */
151
152 /* finv.inc */
153
154 a[*idim + 2] = -a[(*idim << 1) + 2] * (a[*idim + 1] * a[*idim + 2]);
155 a[(*idim << 1) + 1] = -a[(*idim << 1) + 1];
156
157 if (*n != 2) {
158 for (i = 3; i <= *n; ++i) {
159 const int ii = i * (*idim);
160 const int iii = i + ii;
161 const int imi = ii - *idim;
162 const int iimi = i + imi;
163 im2 = i - 2;
164 for (j = 1; j <= im2; ++j) {
165 const int ji = j * (*idim);
166 const int jii = j + ii;
167 s31 = 0.;
168 for (k = j; k <= im2; ++k) {
169 s31 += a[k + ji] * a[i + k * (*idim)];
170 a[jii] += a[j + (k + 1) * (*idim)] * a[k + 1 + ii];
171 } // for k
172 a[i + ji] = -a[iii] * (a[i - 1 + ji] * a[iimi] + s31);
173 a[jii] *= -1;
174 } // for j
175 a[iimi] = -a[iii] * (a[i - 1 + imi] * a[iimi]);
176 a[i - 1 + ii] *= -1;
177 } // for i
178 } // if n!=2
179
180 nm1 = *n - 1;
181 for (i = 1; i <= nm1; ++i) {
182 const int ii = i * (*idim);
183 nmi = *n - i;
184 for (j = 1; j <= i; ++j) {
185 const int ji = j * (*idim);
186 const int iji = i + ji;
187 for (k = 1; k <= nmi; ++k) {
188 a[iji] += a[i + k + ji] * a[i + (i + k) * (*idim)];
189 } // for k
190 } // for j
191
192 for (j = 1; j <= nmi; ++j) {
193 const int ji = j * (*idim);
194 s34 = 0.;
195 for (k = j; k <= nmi; ++k) {
196 s34 += a[i + k + ii + ji] * a[i + (i + k) * (*idim)];
197 } // for k
198 a[i + ii + ji] = s34;
199 } // for j
200 } // for i
201
202 nxch = ir[*n];
203 if (nxch == 0) {
204 return;
205 }
206
207 for (m = 1; m <= nxch; ++m) {
208 k = nxch - m + 1;
209 ij = ir[k];
210 i = ij / 4096;
211 j = ij % 4096;
212 const int ii = i * (*idim);
213 const int ji = j * (*idim);
214 for (k = 1; k <= *n; ++k) {
215 ti = a[k + ii];
216 a[k + ii] = a[k + ji];
217 a[k + ji] = ti;
218 } // for k
219 } // for m
220
221 return;
222} // Dfinv1
ti
EdbTrackP * ti[10]
Definition: RecDispMC_Profiles.C:54

◆ Dinv1()

void Dinv1 ( int *  n,
double *  a,
int *  idim,
double *  ir,
int *  ifail 
)
322{
323 static double det = 0;
324
325 *ifail = 0;
326
327 if (*n < 1 || *n > *idim) {
328 *ifail = -1;
329 return;
330 }
331
332 if ( *n == 1)
333 {
334 if (a[0] == 0.) {
335 *ifail = -1;
336 return;
337 }
338 a[0] = 1. / a[0];
339 return;
340 }
341
342 if ( *n == 2)
343 {
344 det = a[0] * a[3] - a[2] * a[1];
345
346 if (det == 0.) { *ifail = -1; return; }
347
348 double s = 1. / det;
349 double c11 = s * a[3];
350
351 a[2] = -s * a[2];
352 a[1] = -s * a[1];
353 a[3] = s * a[0];
354 a[0] = c11;
355 return;
356 }
357
358
359
360 if ( *n == 3 )
361 {
362 static double t1, t2, t3, temp, s;
363 static double c11, c12, c13, c21, c22, c23, c31, c32, c33;
364
365 /* COMPUTE COFACTORS. */
366 c11 = a[4] * a[8] - a[7] * a[5];
367 c12 = a[7] * a[2] - a[1] * a[8];
368 c13 = a[1] * a[5] - a[4] * a[2];
369 c21 = a[5] * a[6] - a[8] * a[3];
370 c22 = a[8] * a[0] - a[2] * a[6];
371 c23 = a[2] * a[3] - a[5] * a[0];
372 c31 = a[3] * a[7] - a[6] * a[4];
373 c32 = a[6] * a[1] - a[0] * a[7];
374 c33 = a[0] * a[4] - a[3] * a[1];
375
376 t1 = fabs(a[0]);
377 t2 = fabs(a[1]);
378 t3 = fabs(a[2]);
379
380 /* (SET TEMP=PIVOT AND DET=PIVOT*DET.) */
381 if(t1 < t2) {
382 if (t3 < t2) {
383 /* (PIVOT IS A21) */
384 temp = a[1];
385 det = c13 * c32 - c12 * c33;
386 } else {
387 /* (PIVOT IS A31) */
388 temp = a[2];
389 det = c23 * c12 - c22 * c13;
390 }
391 } else {
392 if(t3 < t1) {
393 /* (PIVOT IS A11) */
394 temp = a[0];
395 det = c22 * c33 - c23 * c32;
396 } else {
397 /* (PIVOT IS A31) */
398 temp = a[2];
399 det = c23 * c12 - c22 * c13;
400 }
401 }
402
403 if (det == 0.) {
404 *ifail = -1;
405 return;
406 }
407
408 /* SET ELEMENTS OF INVERSE IN A. */
409 s = temp / det;
410 a[0] = s * c11;
411 a[3] = s * c21;
412 a[6] = s * c31;
413 a[1] = s * c12;
414 a[4] = s * c22;
415 a[7] = s * c32;
416 a[2] = s * c13;
417 a[5] = s * c23;
418 a[8] = s * c33;
419 return;
420 }
421
422 /* Initialized data */
423
424 static int work[100], jfail;
425 for(int i=0; i<*n; ++i) work[i] = 0;
426
427 /* Function Body */
428
429 /* N.GT.3 CASES. FACTORIZE MATRIX AND INVERT. */
430 Dfactir1(n, a, idim, work, ifail, &det, &jfail);
431 if (*ifail != 0)
432 {
433 cerr << "Dfactir failed!!" << endl;
434 return;
435 }
436 Dfinv1(n, a, idim, work, ifail);
437 return;
438} // Dinv1
Dfinv1
void Dfinv1(int *n, double *a, int *idim, int *ir, int *ifail)
Definition: VtSqMatrix.C:131
Dfactir1
void Dfactir1(int *n, double *a, int *idim, int *ir, int *ifail, double *det, int *jfail)
Definition: VtSqMatrix.C:31
s
s
Definition: check_shower.C:55