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

 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
EdbSegP * s
Definition: tlg2pattern.C:32