ATLAS Offline Software
ALFA_SvdCalc.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /*
6  * svdcomp - SVD decomposition routine.
7  * Takes an mxn matrix a and decomposes it into udv, where u,v are
8  * left and right orthogonal transformation matrices, and d is a
9  * diagonal matrix of singular values.
10  *
11  * This routine is adapted from svdecomp.c in XLISP-STAT 2.1 which is
12  * code from Numerical Recipes adapted by Luke Tierney and David Betz.
13  *
14  * Input to dsvd is as follows:
15  * a = mxn matrix to be decomposed, gets overwritten with u
16  * m = row dimension of a
17  * n = column dimension of a
18  * w = returns the vector of singular values of a
19  * v = returns the right orthogonal transformation matrix
20 */
21 
22 #include "ALFA_Geometry/ALFA_SvdCalc.h"
23 #include <cmath>
24 #include <cstdio>
25 #include <cstdlib>
26 
27 static double PYTHAG(double a, double b)
28 {
29  double at = fabs(a), bt = fabs(b), ct, result;
30 
31  if (at > bt) { ct = bt / at; result = at * sqrt(1.0 + ct * ct); }
32  else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); }
33  else result = 0.0;
34  return(result);
35 }
36 
37 
38 int dsvd(double **a, int m, int n, double *w, double **v)
39 {
40  int flag, i, its, j, jj, k, l=0, nm=0;
41  double c, f, h, s, x, y, z;
42  double anorm = 0.0, g = 0.0, scale = 0.0;
43  double *rv1;
44 
45  if (m < n)
46  {
47  fprintf(stderr, "#rows must be > #cols \n");
48  return(0);
49  }
50 
51  rv1 = (double *)malloc((unsigned int) n*sizeof(double));
52 
53 /* Householder reduction to bidiagonal form */
54  for (i = 0; i < n; i++)
55  {
56  /* left-hand reduction */
57  l = i + 1;
58  rv1[i] = scale * g;
59  g = s = scale = 0.0;
60  if (i < m)
61  {
62  for (k = i; k < m; k++)
63  scale += fabs((double)a[k][i]);
64  if (scale)
65  {
66  for (k = i; k < m; k++)
67  {
68  a[k][i] = (double)((double)a[k][i]/scale);
69  s += ((double)a[k][i] * (double)a[k][i]);
70  }
71  f = (double)a[i][i];
72  g = -SIGN(sqrt(s), f);
73  h = f * g - s;
74  a[i][i] = (double)(f - g);
75  if (i != n - 1)
76  {
77  for (j = l; j < n; j++)
78  {
79  for (s = 0.0, k = i; k < m; k++)
80  s += ((double)a[k][i] * (double)a[k][j]);
81  f = s / h;
82  for (k = i; k < m; k++)
83  a[k][j] += (double)(f * (double)a[k][i]);
84  }
85  }
86  for (k = i; k < m; k++)
87  a[k][i] = (double)((double)a[k][i]*scale);
88  }
89  }
90  w[i] = (double)(scale * g);
91 
92  /* right-hand reduction */
93  g = s = scale = 0.0;
94  if (i < m && i != n - 1)
95  {
96  for (k = l; k < n; k++)
97  scale += fabs((double)a[i][k]);
98  if (scale)
99  {
100  for (k = l; k < n; k++)
101  {
102  a[i][k] = (double)((double)a[i][k]/scale);
103  s += ((double)a[i][k] * (double)a[i][k]);
104  }
105  f = (double)a[i][l];
106  g = -SIGN(sqrt(s), f);
107  h = f * g - s;
108  a[i][l] = (double)(f - g);
109  for (k = l; k < n; k++)
110  rv1[k] = (double)a[i][k] / h;
111  if (i != m - 1)
112  {
113  for (j = l; j < m; j++)
114  {
115  for (s = 0.0, k = l; k < n; k++)
116  s += ((double)a[j][k] * (double)a[i][k]);
117  for (k = l; k < n; k++)
118  a[j][k] += (double)(s * rv1[k]);
119  }
120  }
121  for (k = l; k < n; k++)
122  a[i][k] = (double)((double)a[i][k]*scale);
123  }
124  }
125  anorm = MAX(anorm, (fabs((double)w[i]) + fabs(rv1[i])));
126  }
127 
128  /* accumulate the right-hand transformation */
129  for (i = n - 1; i >= 0; i--)
130  {
131  if (i < n - 1)
132  {
133  if (g)
134  {
135  for (j = l; j < n; j++)
136  v[j][i] = (double)(((double)a[i][j] / (double)a[i][l]) / g);
137  /* double division to avoid underflow */
138  for (j = l; j < n; j++)
139  {
140  for (s = 0.0, k = l; k < n; k++)
141  s += ((double)a[i][k] * (double)v[k][j]);
142  for (k = l; k < n; k++)
143  v[k][j] += (double)(s * (double)v[k][i]);
144  }
145  }
146  for (j = l; j < n; j++)
147  v[i][j] = v[j][i] = 0.0;
148  }
149  v[i][i] = 1.0;
150  g = rv1[i];
151  l = i;
152  }
153 
154  /* accumulate the left-hand transformation */
155  for (i = n - 1; i >= 0; i--)
156  {
157  l = i + 1;
158  g = (double)w[i];
159  if (i < n - 1)
160  for (j = l; j < n; j++)
161  a[i][j] = 0.0;
162  if (g)
163  {
164  g = 1.0 / g;
165  if (i != n - 1)
166  {
167  for (j = l; j < n; j++)
168  {
169  for (s = 0.0, k = l; k < m; k++)
170  s += ((double)a[k][i] * (double)a[k][j]);
171  f = (s / (double)a[i][i]) * g;
172  for (k = i; k < m; k++)
173  a[k][j] += (double)(f * (double)a[k][i]);
174  }
175  }
176  for (j = i; j < m; j++)
177  a[j][i] = (double)((double)a[j][i]*g);
178  }
179  else
180  {
181  for (j = i; j < m; j++)
182  a[j][i] = 0.0;
183  }
184  ++a[i][i];
185  }
186 
187  /* diagonalize the bidiagonal form */
188  for (k = n - 1; k >= 0; k--)
189  { /* loop over singular values */
190  for (its = 0; its < 30; its++)
191  { /* loop over allowed iterations */
192  flag = 1;
193  for (l = k; l >= 0; l--)
194  { /* test for splitting */
195  nm = l - 1;
196  if (fabs(rv1[l]) + anorm == anorm)
197  {
198  flag = 0;
199  break;
200  }
201  if (fabs((double)w[nm]) + anorm == anorm)
202  break;
203  }
204  if (flag)
205  {
206  c = 0.0;
207  s = 1.0;
208  for (i = l; i <= k; i++)
209  {
210  f = s * rv1[i];
211  if (fabs(f) + anorm != anorm)
212  {
213  g = (double)w[i];
214  h = PYTHAG(f, g);
215  w[i] = (double)h;
216  h = 1.0 / h;
217  c = g * h;
218  s = (- f * h);
219  for (j = 0; j < m; j++)
220  {
221  y = (double)a[j][nm];
222  z = (double)a[j][i];
223  a[j][nm] = (double)(y * c + z * s);
224  a[j][i] = (double)(z * c - y * s);
225  }
226  }
227  }
228  }
229  z = (double)w[k];
230  if (l == k)
231  { /* convergence */
232  if (z < 0.0)
233  { /* make singular value nonnegative */
234  w[k] = (double)(-z);
235  for (j = 0; j < n; j++)
236  v[j][k] = (-v[j][k]);
237  }
238  break;
239  }
240  if (its >= 30) {
241  free((void*) rv1);
242  fprintf(stderr, "No convergence after 30,000! iterations \n");
243  return(0);
244  }
245 
246  /* shift from bottom 2 x 2 minor */
247  x = (double)w[l];
248  nm = k - 1;
249  y = (double)w[nm];
250  g = rv1[nm];
251  h = rv1[k];
252  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
253  g = PYTHAG(f, 1.0);
254  f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
255 
256  /* next QR transformation */
257  c = s = 1.0;
258  for (j = l; j <= nm; j++)
259  {
260  i = j + 1;
261  g = rv1[i];
262  y = (double)w[i];
263  h = s * g;
264  g = c * g;
265  z = PYTHAG(f, h);
266  rv1[j] = z;
267  c = f / z;
268  s = h / z;
269  f = x * c + g * s;
270  g = g * c - x * s;
271  h = y * s;
272  y = y * c;
273  for (jj = 0; jj < n; jj++)
274  {
275  x = (double)v[jj][j];
276  z = (double)v[jj][i];
277  v[jj][j] = (double)(x * c + z * s);
278  v[jj][i] = (double)(z * c - x * s);
279  }
280  z = PYTHAG(f, h);
281  w[j] = (double)z;
282  if (z)
283  {
284  z = 1.0 / z;
285  c = f * z;
286  s = h * z;
287  }
288  f = (c * g) + (s * y);
289  x = (c * y) - (s * g);
290  for (jj = 0; jj < m; jj++)
291  {
292  y = (double)a[jj][j];
293  z = (double)a[jj][i];
294  a[jj][j] = (double)(y * c + z * s);
295  a[jj][i] = (double)(z * c - y * s);
296  }
297  }
298  rv1[l] = 0.0;
299  rv1[k] = f;
300  w[k] = (double)x;
301  }
302  }
303  free((void*) rv1);
304  return(1);
305 }
ALFA_SvdCalc.h
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
get_generator_info.result
result
Definition: get_generator_info.py:21
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
MAX
#define MAX(x, y)
Definition: ALFA_SvdCalc.h:28
SIGN
#define SIGN(a, b)
Definition: ALFA_SvdCalc.h:29
get_generator_info.stderr
stderr
Definition: get_generator_info.py:40
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
x
#define x
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
beamspotman.n
n
Definition: beamspotman.py:731
extractSporadic.h
list h
Definition: extractSporadic.py:97
master.flag
bool flag
Definition: master.py:29
hist_file_dump.f
f
Definition: hist_file_dump.py:135
TrigInDetValidation_Base.malloc
malloc
Definition: TrigInDetValidation_Base.py:132
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
calibdata.ct
ct
Definition: calibdata.py:418
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
python.PyAthena.v
v
Definition: PyAthena.py:154
a
TList * a
Definition: liststreamerinfos.cxx:10
CalibCoolCompareRT.nm
nm
Definition: CalibCoolCompareRT.py:110
y
#define y
h
dsvd
int dsvd(double **a, int m, int n, double *w, double **v)
Definition: ALFA_SvdCalc.cxx:38
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
python.compressB64.c
def c
Definition: compressB64.py:93
fitman.k
k
Definition: fitman.py:528
Generated on Wed Jan 8 2025 21:06:37 for ATLAS Offline Software by   doxygen 1.8.18