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

Linux Cross Reference
Tina4/src/math/matvec/mat_svd.c

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

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

~ [ 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.