~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/numeric/svd.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)SVD
  2 **/
  3 
  4 #include <stdio.h>
  5 #include <math.h>
  6 #include <tina/sys.h>
  7 #include <tina/sysfuncs.h>
  8 #include <tina/math.h>
  9 
 10 /* smc 2/2/93 */
 11 #ifdef Print
 12 #undef Print
 13 #endif
 14 /* smc 2/2/93 */
 15 
 16 #define Print(x)  { if(out_text!=NULL) out_text(x); }
 17 
 18 /* This macro is now redundant as function exists but kept incase
 19  * function is bad re-implementation of macro */
 20 #define Radius(u,v)          ( (Au=Abs((u))) > (Av=Abs((v))) ? \
 21                                (Aw=Av/Au, Au*sqrt(1.+Aw*Aw)) : \
 22                                (Av ? (Aw=Au/Av, Av*sqrt(1.+Aw*Aw)) : 0.0) )
 23 
 24 #define Sign(u,v)               ( (v)>=0.0 ? Abs(u) : -Abs(u) )
 25 
 26 #define Max(u,v)                ( (u)>(v)?  (u) : (v) )
 27 
 28 #define Abs(x)                  ( (x)>0.0?  (x) : (-(x)) )
 29 
 30 static void (*out_text) () = NULL;
 31 
 32 static double radius(double u, double v)
 33 {
 34     double  Au, Av, Aw;
 35 
 36     Au = Abs(u);
 37     Av = Abs(v);
 38 
 39     if (Au > Av)
 40     {
 41         Aw = Av / Au;
 42         return (Au * sqrt(1. + Aw * Aw));
 43     } else
 44     {
 45         if (Av)
 46         {
 47             Aw = Au / Av;
 48             return (Av * sqrt(1. + Aw * Aw));
 49         } else
 50             return 0.0;
 51     }
 52 }
 53 
 54 
 55 
 56 void    svd_err_func(void (*text) ( /* ??? */ ))
 57 {
 58     out_text = text;
 59 }
 60 
 61 /*************************** SVDcmp *****************************************
 62 * Given matrix A[m][n], m>=n, using svd decomposition A = U W V' to get     *
 63 * U[m][n], W[n][n] and V[n][n], where U occupies the position of A.         *
 64 * NOTE: if m<n, A should be filled up to square with zero rows.             *
 65 *       A[m][n] has been destroyed by U[m][n] after the decomposition.      *
 66 ****************************************************************************/
 67 void    SVDcmp(double **a, int m, int n, double *w, double **v)
 68 {
 69     /* BUG `nm' may be used uninitialized in this function */
 70     int     flag, i, its, j, jj, k, l, nm, nm1 = n - 1, mm1 = m - 1;
 71     double  c, f, h, s, x, y, z;
 72     double  anorm = 0.0, g = 0.0, scale = 0.0;
 73     double *rv1;
 74 
 75     if (m < n)
 76         Print("SVDCMP: You must augment A with extra zero rows");
 77     rv1 = (double *) ralloc((unsigned) n * sizeof(double));
 78 
 79     /* Householder reduction to bidigonal form */
 80     for (i = 0; i < n; i++)
 81     {
 82         l = i + 1;
 83         rv1[i] = scale * g;
 84         g = s = scale = 0.0;
 85         if (i < m)
 86         {
 87             for (k = i; k < m; k++)
 88                 scale += Abs(a[k][i]);
 89             if (scale)
 90             {
 91                 for (k = i; k < m; k++)
 92                 {
 93                     a[k][i] /= scale;
 94                     s += a[k][i] * a[k][i];
 95                 }
 96                 f = a[i][i];
 97                 g = -Sign(sqrt(s), f);
 98                 h = f * g - s;
 99                 a[i][i] = f - g;
100                 if (i != nm1)
101                 {
102                     for (j = l; j < n; j++)
103                     {
104                         for (s = 0.0, k = i; k < m; k++)
105                             s += a[k][i] * a[k][j];
106                         f = s / h;
107                         for (k = i; k < m; k++)
108                             a[k][j] += f * a[k][i];
109                     }
110                 }
111                 for (k = i; k < m; k++)
112                     a[k][i] *= scale;
113             }
114         }
115         w[i] = scale * g;
116         g = s = scale = 0.0;
117         if (i < m && i != nm1)
118         {
119             for (k = l; k < n; k++)
120                 scale += Abs(a[i][k]);
121             if (scale)
122             {
123                 for (k = l; k < n; k++)
124                 {
125                     a[i][k] /= scale;
126                     s += a[i][k] * a[i][k];
127                 }
128                 f = a[i][l];
129                 g = -Sign(sqrt(s), f);
130                 h = f * g - s;
131                 a[i][l] = f - g;
132                 for (k = l; k < n; k++)
133                     rv1[k] = a[i][k] / h;
134                 if (i != mm1)
135                 {
136                     for (j = l; j < m; j++)
137                     {
138                         for (s = 0.0, k = l; k < n; k++)
139                             s += a[j][k] * a[i][k];
140                         for (k = l; k < n; k++)
141                             a[j][k] += s * rv1[k];
142                     }
143                 }
144                 for (k = l; k < n; k++)
145                     a[i][k] *= scale;
146             }
147         }
148         anorm = Max(anorm, (Abs(w[i]) + Abs(rv1[i])));
149     }
150 
151     /* Accumulation of right-hand transformations */
152     for (i = n - 1; i >= 0; i--)
153     {
154         if (i < nm1)
155         {
156             if (g)
157             {
158                 /* double division to avoid possible underflow */
159                 for (j = l; j < n; j++)
160                     v[j][i] = (a[i][j] / a[i][l]) / g;
161                 for (j = l; j < n; j++)
162                 {
163                     for (s = 0.0, k = l; k < n; k++)
164                         s += a[i][k] * v[k][j];
165                     for (k = l; k < n; k++)
166                         v[k][j] += s * v[k][i];
167                 }
168             }
169             for (j = l; j < n; j++)
170                 v[i][j] = v[j][i] = 0.0;
171         }
172         v[i][i] = 1.0;
173         g = rv1[i];
174         l = i;
175     }
176     /* Accumulation of left-hand transformations */
177     for (i = n - 1; i >= 0; i--)
178     {
179         l = i + 1;
180         g = w[i];
181         if (i < nm1)
182             for (j = l; j < n; j++)
183                 a[i][j] = 0.0;
184         if (g)
185         {
186             g = 1.0 / g;
187             if (i != nm1)
188             {
189                 for (j = l; j < n; j++)
190                 {
191                     for (s = 0.0, k = l; k < m; k++)
192                         s += a[k][i] * a[k][j];
193                     f = (s / a[i][i]) * g;
194                     for (k = i; k < m; k++)
195                         a[k][j] += f * a[k][i];
196                 }
197             }
198             for (j = i; j < m; j++)
199                 a[j][i] *= g;
200         } else
201             for (j = i; j < m; j++)
202                 a[j][i] = 0.0;
203         ++a[i][i];
204     }
205     /* diagonalization of the bidigonal form */
206     for (k = n - 1; k >= 0; k--)
207     {                           /* loop over singlar values */
208         for (its = 0; its < 30; its++)
209         {                       /* loop over allowed iterations */
210             flag = 1;
211             for (l = k; l >= 0; l--)
212             {                   /* test for splitting */
213                 nm = l - 1;     /* note that rv1[l] is always zero */
214                 if (Abs(rv1[l]) + anorm == anorm)
215                 {
216                     flag = 0;
217                     break;
218                 }
219                 if (Abs(w[nm]) + anorm == anorm)
220                     break;
221             }
222             if (flag)
223             {
224                 c = 0.0;        /* cancellation of rv1[l], if l>1 */
225                 s = 1.0;
226                 for (i = l; i <= k; i++)
227                 {
228                     f = s * rv1[i];
229                     if (Abs(f) + anorm != anorm)
230                     {
231                         g = w[i];
232                         h = radius(f, g);
233                         w[i] = h;
234                         h = 1.0 / h;
235                         c = g * h;
236                         s = (-f * h);
237                         for (j = 0; j < m; j++)
238                         {
239                             y = a[j][nm];
240                             z = a[j][i];
241                             a[j][nm] = y * c + z * s;
242                             a[j][i] = z * c - y * s;
243                         }
244                     }
245                 }
246             }
247             z = w[k];
248             if (l == k)
249             {                   /* convergence */
250                 if (z < 0.0)
251                 {
252                     w[k] = -z;
253                     for (j = 0; j < n; j++)
254                         v[j][k] = (-v[j][k]);
255                 }
256                 break;
257             }
258             if (its == 30)
259                 Print("No convergence in 30 SVDCMP iterations");
260             x = w[l];           /* shift from bottom 2-by-2 minor */
261             nm = k - 1;
262             y = w[nm];
263             g = rv1[nm];
264             h = rv1[k];
265             f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
266             g = radius(f, 1.0);
267             /* next QR transformation */
268             f = ((x - z) * (x + z) + h * ((y / (f + Sign(g, f))) - h)) / x;
269             c = s = 1.0;
270             for (j = l; j <= nm; j++)
271             {
272                 i = j + 1;
273                 g = rv1[i];
274                 y = w[i];
275                 h = s * g;
276                 g = c * g;
277                 z = radius(f, h);
278                 rv1[j] = z;
279                 c = f / z;
280                 s = h / z;
281                 f = x * c + g * s;
282                 g = g * c - x * s;
283                 h = y * s;
284                 y = y * c;
285                 for (jj = 0; jj < n; jj++)
286                 {
287                     x = v[jj][j];
288                     z = v[jj][i];
289                     v[jj][j] = x * c + z * s;
290                     v[jj][i] = z * c - x * s;
291                 }
292                 z = radius(f, h);
293                 w[j] = z;       /* rotation can be arbitrary id z=0 */
294                 if (z)
295                 {
296                     z = 1.0 / z;
297                     c = f * z;
298                     s = h * z;
299                 }
300                 f = (c * g) + (s * y);
301                 x = (c * y) - (s * g);
302                 for (jj = 0; jj < m; jj++)
303                 {
304                     y = a[jj][j];
305                     z = a[jj][i];
306                     a[jj][j] = y * c + z * s;
307                     a[jj][i] = z * c - y * s;
308                 }
309             }
310             rv1[l] = 0.0;
311             rv1[k] = f;
312             w[k] = x;
313         }
314     }
315     rfree((void *) rv1);
316 }
317 
318 /******************************* SVDbksb ****************************
319 * Given A[m][n], b[m], solves A x = b in the svd form U W V'x = b   *
320 * so x = V U'b/W                                                    *
321 * No input quantities are destroyed, so sequential calling is OK    *
322 ********************************************************************/
323 void    SVDbksb(double **u, double *w, double **v, int m, int n, double *b, double *x)
324 {
325     int     jj, j, i;
326     double  s, *tmp;
327 
328     tmp = (double *) ralloc((unsigned) n * sizeof(double));
329     for (j = 0; j < n; j++)
330     {                           /* calculate <U,b> */
331         s = 0.0;
332         if (w[j])
333         {                       /* nonzero result only if nonzero w[j] */
334             for (i = 0; i < m; i++)
335                 s += u[i][j] * b[i];
336             s /= w[j];          /* premultiplied by inverse W */
337         }
338         tmp[j] = s;
339     }
340     for (j = 0; j < n; j++)
341     {                           /* V(U'b/W) */
342         s = 0.0;
343         for (jj = 0; jj < n; jj++)
344             s += v[j][jj] * tmp[jj];
345         x[j] = s;
346     }
347     rfree((void *) tmp);
348 }
349 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.