Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Macros | Functions
NR.h File Reference

Go to the source code of this file.

Macros

#define NRGLOBAL   extern
 
#define _NR_H_
 

Functions

NRGLOBAL void nrerror ()
 
NRGLOBAL float * vector ()
 
NRGLOBAL void gaussj (double **a, int n, double *b, int m)
 
NRGLOBAL void ludcmp (double **a, int matsize, int *b, double *c)
 
NRGLOBAL void lubksb (double **a, int matsize, int *b, double *c)
 
NRGLOBAL void svdcmp (double **a, int matrow, int matcol, double *w, double **v)
 
NRGLOBAL void svbksb (double **a, double *w, double **v, int matrow, int matcol, double *V, double *ChDen)
 

Macro Definition Documentation

◆ _NR_H_

#define _NR_H_

Definition at line 140 of file NR.h.

◆ NRGLOBAL

#define NRGLOBAL   extern

Definition at line 7 of file NR.h.

Function Documentation

◆ gaussj()

NRGLOBAL void gaussj ( double **  a,
int  n,
double *  b,
int  m 
)

◆ lubksb()

NRGLOBAL void lubksb ( double **  a,
int  matsize,
int *  b,
double *  c 
)

Definition at line 90 of file luc.c.

90 {
91 int i, ii = 0, ip, j;
92 double sum;
93
94 for (i = 1; i <= n; i++) /* forward substitution. */
95 {
96 ip = index[i];
97 sum = b[ip];
98 b[ip] = b[i];
99 if (ii) {
100#ifdef _OPENMP
101#pragma omp parallel for private(j) reduction(- : sum)
102#endif
103 for (j = ii; j <= i; j++) sum = sum - a[i][j] * b[j];
104 } else if (sum)
105 ii = i;
106 b[i] = sum;
107 }
108
109 for (i = n; i >= 1; i--) /* back substitution. */
110 {
111 sum = b[i];
112#ifdef _OPENMP
113#pragma omp parallel for private(j) reduction(- : sum)
114#endif
115 for (j = i + 1; j <= n; j++) sum = sum - a[i][j] * b[j];
116 b[i] = sum / a[i][i];
117 }
118}

Referenced by InvertMatrix().

◆ ludcmp()

NRGLOBAL void ludcmp ( double **  a,
int  matsize,
int *  b,
double *  c 
)

Definition at line 12 of file luc.c.

12 {
13 int i, j, k;
14 int imax = 1;
15 double big, sum, dum;
16 double *vv; /* vv stores the implicit scaling of each row */
17
18 vv = dvector(1, n); /* No looping yet. */
19 *d = 1.0; /* loop over rows to get implicit scaling info. */
20#ifdef _OPENMP
21#pragma omp parallel for private(j, big)
22#endif
23 for (i = 1; i <= n; i++) {
24 big = 0.0;
25 for (j = 1; j <= n; j++)
26 if (fabs(a[i][j]) > big) big = fabs(a[i][j]);
27 if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
28 vv[i] = 1.0 / big;
29 }
30
31 for (j = 1; j <= n; j++) /* outer loop of Crout's method */
32 {
33 for (i = 1; i < j; i++) {
34 sum = a[i][j];
35#ifdef _OPENMP
36#pragma omp parallel for private(k) reduction(- : sum)
37#endif
38 for (k = 1; k < i; k++) // OMPChk: parallelization may not help much
39 sum = sum - a[i][k] * a[k][j];
40 a[i][j] = sum;
41 }
42
43 big = 0.0; /* search for largest pivotal element. */
44 for (i = j; i <= n; i++) {
45 sum = a[i][j];
46#ifdef _OPENMP
47#pragma omp parallel for private(k) reduction(- : sum)
48#endif
49 for (k = 1; k < j; k++) // OMPChk: parallelization may not help much
50 sum = sum - a[i][k] * a[k][j];
51 a[i][j] = sum;
52
53 dum = vv[i] * fabs(sum);
54 if (dum >= big) { /* is the figure of merit better */
55 big = dum; /* than the best? */
56 imax = i;
57 }
58 }
59
60 if (j != imax) {
61#ifdef _OPENMP
62#pragma omp parallel for private(k, dum)
63#endif
64 for (k = 1; k <= n; k++) {
65 dum = a[imax][k];
66 a[imax][k] = a[j][k];
67 a[j][k] = dum;
68 }
69 *d = -(*d); /* change the parity of d. */
70 vv[imax] = vv[j]; /* interchange the scale factor. */
71 }
72 index[j] = imax;
73 if (a[j][j] == 0.0) /* for some applications on singular */
74 a[j][j] = TINY; /* matrices, it is desirable to */
75 /* substitute TINY for zero. */
76
77 if (j != n) /* finally, divide by the pivot element. */
78 {
79 dum = 1.0 / a[j][j];
80#ifdef _OPENMP
81#pragma omp parallel for private(i) // OMPCheck: may not help much
82#endif
83 for (i = j + 1; i <= n; i++) a[i][j] = a[i][j] * dum;
84 }
85 }
86
87 free_dvector(vv, 1, n);
88}
nrerror
NRGLOBAL void nrerror()
TINY
#define TINY
Definition: luc.c:6
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
free_dvector
void free_dvector(double *v, long nl, long)
Definition: nrutil.c:470
dvector
double * dvector(long nl, long nh)
Definition: nrutil.c:63

Referenced by InvertMatrix().

◆ nrerror()

NRGLOBAL void nrerror ( )

Referenced by convert_matrix(), cvector(), d3tensor(), d4tensor(), dmatrix(), dvector(), f3tensor(), fmatrix(), fvector(), imatrix(), ivector(), ludcmp(), lvector(), submatrix(), and svdcmp().

◆ svbksb()

NRGLOBAL void svbksb ( double **  a,
double *  w,
double **  v,
int  matrow,
int  matcol,
double *  V,
double *  ChDen 
)

◆ svdcmp()

NRGLOBAL void svdcmp ( double **  a,
int  matrow,
int  matcol,
double *  w,
double **  v 
)

Definition at line 27 of file svdcmp.c.

27 {
28 int flag, i, its, j, jj, k, l, nm;
29 double c, f, h, s, x, y, z;
30 double anorm = 0.0, g = 0.0, scale = 0.0;
31 double *rv1;
32
33 if (m < n) nrerror("SVDCMP: A must be augmented with extra zeros.");
34
35 rv1 = dvector(1, n);
36
37 for (i = 1; i <= n; i++) /* householder reduction to bidiagonal form */
38 {
39 l = i + 1;
40 rv1[i] = scale * g;
41 g = s = scale = 0.0;
42
43 if (i <= m) {
44#ifdef _OPENMP
45#pragma omp parallel for private(k) reduction(+ : scale)
46#endif
47 for (k = i; k <= m; k++) scale += fabs(a[k][i]);
48 if (scale) {
49#ifdef _OPENMP
50#pragma omp parallel for private(k) reduction(+ : s)
51#endif
52 for (k = i; k <= m; k++) {
53 a[k][i] /= scale;
54 s += a[k][i] * a[k][i];
55 }
56
57 f = a[i][i];
58 g = -SIGNS(sqrt(s), f);
59 h = f * g - s;
60 a[i][i] = f - g;
61
62 if (i != n) {
63 for (j = l; j <= n; j++) {
64 s = 0.0;
65#ifdef _OPENMP
66#pragma omp parallel for private(k) reduction(+ : s)
67#endif
68 for (k = i; k <= m; k++) s += a[k][i] * a[k][j];
69 f = s / h;
70#ifdef _OPENMP
71#pragma omp parallel for private(k) // OMPCheck
72#endif
73 for (k = i; k <= m; k++) a[k][j] += f * a[k][i];
74 }
75 }
76#ifdef _OPENMP
77#pragma omp parallel for private(k)
78#endif
79 for (k = i; k <= m; k++) a[k][i] *= scale;
80 }
81 }
82 w[i] = scale * g;
83 g = s = scale = 0.0;
84 if (i <= m && i != n) {
85#ifdef _OPENMP
86#pragma omp parallel for private(k) reduction(+ : scale)
87#endif
88 for (k = l; k <= n; k++) scale += fabs(a[i][k]);
89 if (scale) {
90#ifdef _OPENMP
91#pragma omp parallel for private(k) reduction(+ : s)
92#endif
93 for (k = l; k <= n; k++) {
94 a[i][k] /= scale;
95 s += a[i][k] * a[i][k];
96 }
97 f = a[i][l];
98 g = -SIGNS(sqrt(s), f);
99 h = f * g - s;
100 a[i][l] = f - g;
101#ifdef _OPENMP
102#pragma omp parallel for private(k)
103#endif
104 for (k = l; k <= n; k++) rv1[k] = a[i][k] / h;
105 if (i != m) {
106 for (j = l; j <= m; j++) {
107 s = 0.0;
108#ifdef _OPENMP
109#pragma omp parallel for private(k) reduction(+ : s)
110#endif
111 for (k = l; k <= n; k++) s += a[j][k] * a[i][k];
112#ifdef _OPENMP
113#pragma omp parallel for private(k) // OMPCheck
114#endif
115 for (k = l; k <= n; k++) a[j][k] += s * rv1[k];
116 }
117 }
118#ifdef _OPENMP
119#pragma omp parallel for private(k)
120#endif
121 for (k = l; k <= n; k++) a[i][k] *= scale;
122 }
123 }
124 anorm = MAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
125 }
126
127 for (i = n; i >= 1; i--) /* acumulation of right-hand transformations */
128 {
129 if (i < n) {
130 if (g) {
131#ifdef _OPENMP
132#pragma omp parallel for private(j)
133#endif
134 for (j = l; j <= n; j++) /* double division to avoid possible */
135 v[j][i] = (a[i][j] / a[i][l]) / g; /* underflow */
136 for (j = l; j <= n; j++) {
137 s = 0.0;
138#ifdef _OPENMP
139#pragma omp parallel for private(k) reduction(+ : s)
140#endif
141 for (k = l; k <= n; k++) s += a[i][k] * v[k][j];
142#ifdef _OPENMP
143#pragma omp parallel for private(k) // OMPCheck
144#endif
145 for (k = l; k <= n; k++) v[k][j] += s * v[k][i];
146 }
147 }
148#ifdef _OPENMP
149#pragma omp parallel for private(j)
150#endif
151 for (j = l; j <= n; j++) v[i][j] = v[j][i] = 0.0;
152 }
153 v[i][i] = 1.0;
154 g = rv1[i];
155 l = i;
156 }
157
158 for (i = n; i >= 1; i--) /* accumulation of the left-hand transformations */
159 {
160 l = i + 1;
161 g = w[i];
162 if (i < n) {
163#ifdef _OPENMP
164#pragma omp parallel for private(j)
165#endif
166 for (j = l; j <= n; j++) a[i][j] = 0.0;
167 }
168 if (g) {
169 g = 1.0 / g;
170 if (i != n) {
171 for (j = l; j <= n; j++) {
172 s = 0.0;
173#ifdef _OPENMP
174#pragma omp parallel for private(k) reduction(+ : s)
175#endif
176 for (k = l; k <= m; k++) s += a[k][i] * a[k][j];
177 f = (s / a[i][i]) * g;
178#ifdef _OPENMP
179#pragma omp parallel for private(k) // OMPCheck
180#endif
181 for (k = i; k <= m; k++) a[k][j] += f * a[k][i];
182 }
183 }
184#ifdef _OPENMP
185#pragma omp parallel for private(j)
186#endif
187 for (j = i; j <= m; j++) a[j][i] *= g;
188 } else {
189#ifdef _OPENMP
190#pragma omp parallel for private(j)
191#endif
192 for (j = i; j <= m; j++) a[j][i] = 0.0;
193 }
194
195 ++a[i][i];
196 }
197
198 for (k = n; k >= 1; k--) /* diagonalization of the bidiagonal form */
199 for (its = 1; its <= MAX_ITS; its++) /* loop over alllowed values */
200 {
201 flag = 1;
202 for (l = k; l >= 1; l--) /* test for splitting */
203 {
204 nm = l - 1; /* note that rv1 is always zero */
205 if ((double)(fabs(rv1[l]) + anorm) == anorm) {
206 flag = 0;
207 break;
208 }
209 if ((double)(fabs(w[nm]) + anorm) == anorm) break;
210 }
211 if (flag) {
212 c = 0.0; /* cancellation of rv1[l], if l > 1 */
213 s = 1.0;
214 for (i = l; i <= k; i++) {
215 f = s * rv1[i];
216 rv1[i] = c * rv1[i];
217 if ((double)(fabs(f) + anorm) == anorm) break;
218 g = w[i];
219 h = PYTHAG(f, g);
220 w[i] = h;
221 h = 1.0 / h;
222 c = g * h;
223 s = (-f * h);
224#ifdef _OPENMP
225#pragma omp parallel for private(j, y, z) // OMPCheck
226#endif
227 for (j = 1; j <= m; j++) {
228 y = a[j][nm];
229 z = a[j][i];
230 a[j][nm] = y * c + z * s;
231 a[j][i] = z * c - y * s;
232 }
233 }
234 }
235 z = w[k];
236 if (l == k) /* convergence */
237 {
238 if (z < 0.0) /* singular value is made nonnegative */
239 {
240 w[k] = -z;
241#ifdef _OPENMP
242#pragma omp parallel for private(j)
243#endif
244 for (j = 1; j <= n; j++) v[j][k] = (-v[j][k]);
245 }
246 break;
247 }
248
249 if (its == MAX_ITS)
250 nrerror("no convergence in 1000 SVDCMP iterations.\n");
251
252 x = w[l]; /* shift from bottom 2-by-2 minor */
253 nm = k - 1;
254 y = w[nm];
255 g = rv1[nm];
256 h = rv1[k];
257 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
258 g = PYTHAG(f, 1.0);
259 f = ((x - z) * (x + z) + h * ((y / (f + SIGNS(g, f))) - h)) / x;
260
261 c = s = 1.0; /* next QR transformation */
262 for (j = l; j <= nm; j++) {
263 i = j + 1;
264 g = rv1[i];
265 y = w[i];
266 h = s * g;
267 g *= c;
268 z = PYTHAG(f, h);
269 rv1[j] = z;
270 c = f / z;
271 s = h / z;
272 f = x * c + g * s;
273 g = g * c - x * s;
274 h = y * s;
275 y *= c;
276#ifdef _OPENMP
277#pragma omp parallel for private(jj, x, z) // OMPCheck
278#endif
279 for (jj = 1; jj <= n; jj++) {
280 x = v[jj][j];
281 z = v[jj][i];
282 v[jj][j] = x * c + z * s;
283 v[jj][i] = z * c - x * s;
284 }
285
286 z = PYTHAG(f, h);
287 w[j] = z; /* rotation can be arbitrary if z = 0 */
288 if (z) {
289 z = 1.0 / z;
290 c = f * z;
291 s = h * z;
292 }
293 f = c * g + s * y;
294 x = c * y - s * g;
295#ifdef _OPENMP
296#pragma omp parallel for private(jj, y, z) // OMPCheck
297#endif
298 for (jj = 1; jj <= m; jj++) {
299 y = a[jj][j];
300 z = a[jj][i];
301 a[jj][j] = y * c + z * s;
302 a[jj][i] = z * c - y * s;
303 }
304 }
305 rv1[l] = 0.0;
306 rv1[k] = f;
307 w[k] = x;
308 }
309
310 free_dvector(rv1, 1, n);
311}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
SIGNS
#define SIGNS(a, b)
Definition: svdcmp.c:24
MAX_ITS
#define MAX_ITS
Definition: svdcmp.c:25
PYTHAG
#define PYTHAG(a, b)
Definition: svdcmp.c:18
MAX
#define MAX(a, b)
Definition: svdcmp.c:22

Referenced by DecomposeMatrixSVD().

◆ vector()

NRGLOBAL float * vector ( )