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
svdcmp.c File Reference
#include <math.h>
#include "NR.h"

Go to the source code of this file.

Macros

#define PYTHAG(a, b)
 
#define MAX(a, b)    (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
 
#define SIGNS(a, b)   ((b) >= 0.0 ? fabs(a) : -fabs(a))
 
#define MAX_ITS   1000
 

Functions

void svdcmp (double **a, int m, int n, double *w, double **v)
 

Macro Definition Documentation

◆ MAX

#define MAX (   a,
  b 
)     (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))

Definition at line 22 of file svdcmp.c.

◆ MAX_ITS

#define MAX_ITS   1000

Definition at line 25 of file svdcmp.c.

◆ PYTHAG

#define PYTHAG (   a,
  b 
)
Value:
((at = fabs(a)) > (bt = fabs(b)) \
? (ct = bt / at, at * sqrt(1.0 + ct * ct)) \
: (bt ? (ct = at / bt, bt * sqrt(1.0 + ct * ct)) : 0.0))

Definition at line 18 of file svdcmp.c.

◆ SIGNS

#define SIGNS (   a,
  b 
)    ((b) >= 0.0 ? fabs(a) : -fabs(a))

Definition at line 24 of file svdcmp.c.

Function Documentation

◆ svdcmp()

void svdcmp ( double **  a,
int  m,
int  n,
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}
nrerror
NRGLOBAL void nrerror()
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
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
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().