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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathMatv_svd.c

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathMatv_svd.c,v $
 27  * Date    :  $Date: 2005/01/09 17:49:25 $
 28  * Version :  $Revision: 1.6 $
 29  * CVS Id  :  $Id: mathMatv_svd.c,v 1.6 2005/01/09 17:49:25 paul Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes : Matrix singular value decomposition.
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief  Matrix Singular Value Decomposition, using Mat.     
 40  * 
 41  *  NB, this version is for Mat and Vec, see mathNum_svd.c for an alternative version.
 42  * 
 43  *  See Numerical Recipes (3rd ed) Chapter 2.6 for more info on the method.
 44 */
 45 
 46 #include "mathMatv_svd.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <tina/sys/sysDef.h>
 54 #include <tina/sys/sysPro.h>
 55 #include <tina/math/math_MatvDef.h>
 56 #include <tina/math/mathMatv_mat.h>
 57 
 58 #define SIGN(u,v)               ( (v)>=0.0 ? fabs(u) : -fabs(u) )
 59 
 60 static double   radius(double u, double v)
 61 {
 62         double          w;
 63 
 64         u = fabs(u);
 65         v = fabs(v);
 66         if (u > v)
 67         {
 68                 w = v / u;
 69                 return (u * sqrt(1. + w * w));
 70         } else
 71         {
 72                 if (v)
 73                 {
 74                         w = u / v;
 75                         return (v * sqrt(1. + w * w));
 76                 } else
 77                         return 0.0;
 78         }
 79 }
 80 
 81 /*
 82  Given matrix a[m][n], m>=n, using svd decomposition a = p d q' to get
 83  p[m][n], diag d[n] and q[n][n].
 84 */
 85 void            svd(int m, int n, double **a, double **p, double *d, double **q)
 86 {
 87         int             flag, i, its, j, jj, k, l=0, nm=0, nm1 = n - 1, mm1 = m - 1;
 88         double          c, f, h, s, x, y, z;
 89         double          anorm = 0, g = 0, scale = 0;
 90         double         *r = tvector_alloc(0, n, double);
 91 
 92         for (i = 0; i < m; i++)
 93                 for (j = 0; j < n; j++)
 94                         p[i][j] = a[i][j];
 95         for (i = m; i < n; i++)
 96                 for (j = 0; j < n; j++)
 97                         p[i][j] = 0;
 98 
 99         /* Householder reduction to bidigonal form */
100         for (i = 0; i < n; i++)
101         {
102                 l = i + 1;
103                 r[i] = scale * g;
104                 g = s = scale = 0.0;
105                 if (i < m)
106                 {
107                         for (k = i; k < m; k++)
108                                 scale += fabs(p[k][i]);
109                         if (scale)
110                         {
111                                 for (k = i; k < m; k++)
112                                 {
113                                         p[k][i] /= scale;
114                                         s += p[k][i] * p[k][i];
115                                 }
116                                 f = p[i][i];
117                                 g = -SIGN(sqrt(s), f);
118                                 h = f * g - s;
119                                 p[i][i] = f - g;
120                                 if (i != nm1)
121                                 {
122                                         for (j = l; j < n; j++)
123                                         {
124                                                 for (s = 0.0, k = i; k < m; k++)
125                                                         s += p[k][i] * p[k][j];
126                                                 f = s / h;
127                                                 for (k = i; k < m; k++)
128                                                         p[k][j] += f * p[k][i];
129                                         }
130                                 }
131                                 for (k = i; k < m; k++)
132                                         p[k][i] *= scale;
133                         }
134                 }
135                 d[i] = scale * g;
136                 g = s = scale = 0.0;
137                 if (i < m && i != nm1)
138                 {
139                         for (k = l; k < n; k++)
140                                 scale += fabs(p[i][k]);
141                         if (scale)
142                         {
143                                 for (k = l; k < n; k++)
144                                 {
145                                         p[i][k] /= scale;
146                                         s += p[i][k] * p[i][k];
147                                 }
148                                 f = p[i][l];
149                                 g = -SIGN(sqrt(s), f);
150                                 h = f * g - s;
151                                 p[i][l] = f - g;
152                                 for (k = l; k < n; k++)
153                                         r[k] = p[i][k] / h;
154                                 if (i != mm1)
155                                 {
156                                         for (j = l; j < m; j++)
157                                         {
158                                                 for (s = 0.0, k = l; k < n; k++)
159                                                         s += p[j][k] * p[i][k];
160                                                 for (k = l; k < n; k++)
161                                                         p[j][k] += s * r[k];
162                                         }
163                                 }
164                                 for (k = l; k < n; k++)
165                                         p[i][k] *= scale;
166                         }
167                 }
168                 anorm = MAX(anorm, fabs(d[i]) + fabs(r[i]));
169         }
170 
171         /* Accumulation of right-hand transformations */
172         for (i = n - 1; i >= 0; i--)
173         {
174                 if (i < nm1)
175                 {
176                         if (g)
177                         {
178                                 for (j = l; j < n; j++)
179                                         q[j][i] = (p[i][j] / p[i][l]) / g;
180                                 for (j = l; j < n; j++)
181                                 {
182                                         for (s = 0.0, k = l; k < n; k++)
183                                                 s += p[i][k] * q[k][j];
184                                         for (k = l; k < n; k++)
185                                                 q[k][j] += s * q[k][i];
186                                 }
187                         }
188                         for (j = l; j < n; j++)
189                                 q[i][j] = q[j][i] = 0.0;
190                 }
191                 q[i][i] = 1.0;
192                 g = r[i];
193                 l = i;
194         }
195         /* Accumulation of left-hand transformations */
196         for (i = n - 1; i >= 0; i--)
197         {
198                 l = i + 1;
199                 g = d[i];
200                 if (i < nm1)
201                         for (j = l; j < n; j++)
202                                 p[i][j] = 0.0;
203                 if (g)
204                 {
205                         g = 1.0 / g;
206                         if (i != nm1)
207                         {
208                                 for (j = l; j < n; j++)
209                                 {
210                                         for (s = 0.0, k = l; k < m; k++)
211                                                 s += p[k][i] * p[k][j];
212                                         f = (s / p[i][i]) * g;
213                                         for (k = i; k < m; k++)
214                                                 p[k][j] += f * p[k][i];
215                                 }
216                         }
217                         for (j = i; j < m; j++)
218                                 p[j][i] *= g;
219                 } else
220                         for (j = i; j < m; j++)
221                                 p[j][i] = 0.0;
222                 ++p[i][i];
223         }
224         /* diagonalization of the bidigonal form */
225         for (k = n - 1; k >= 0; k--)
226         {                       /* loop over singlar values */
227                 for (its = 0; its < 30; its++)
228                 {               /* loop over allowed iterations */
229                         flag = 1;
230                         for (l = k; l >= 0; l--)
231                         {       /* test for splitting */
232                                 nm = l - 1;     /* note that r[l] is always
233                                                  * zero */
234                                 if (fabs(r[l]) + anorm == anorm)
235                                 {
236                                         flag = 0;
237                                         break;
238                                 }
239                                 if (fabs(d[nm]) + anorm == anorm)
240                                         break;
241                         }
242                         if (flag)
243                         {
244                                 c = 0.0;        /* cancellation of r[l], if
245                                                  * l>1 */
246                                 s = 1.0;
247                                 for (i = l; i <= k; i++)
248                                 {
249                                         f = s * r[i];
250                                         if (fabs(f) + anorm != anorm)
251                                         {
252                                                 g = d[i];
253                                                 h = radius(f, g);
254                                                 d[i] = h;
255                                                 h = 1.0 / h;
256                                                 c = g * h;
257                                                 s = (-f * h);
258                                                 for (j = 0; j < m; j++)
259                                                 {
260                                                         y = p[j][nm];
261                                                         z = p[j][i];
262                                                         p[j][nm] = y * c + z * s;
263                                                         p[j][i] = z * c - y * s;
264                                                 }
265                                         }
266                                 }
267                         }
268                         z = d[k];
269                         if (l == k)
270                         {       /* convergence */
271                                 if (z < 0.0)
272                                 {
273                                         d[k] = -z;
274                                         for (j = 0; j < n; j++)
275                                                 q[j][k] = (-q[j][k]);
276                                 }
277                                 break;
278                         }
279                         if (its == 30)
280                         {
281                                 error("svd: No convergence in 30 svd iterations", non_fatal);
282                                 return;
283                         }
284                         x = d[l];       /* shift from bottom 2-by-2 minor */
285                         nm = k - 1;
286                         y = d[nm];
287                         g = r[nm];
288                         h = r[k];
289                         f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
290                         g = radius(f, 1.0);
291                         /* next QR transformation */
292                         f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
293                         c = s = 1.0;
294                         for (j = l; j <= nm; j++)
295                         {
296                                 i = j + 1;
297                                 g = r[i];
298                                 y = d[i];
299                                 h = s * g;
300                                 g = c * g;
301                                 z = radius(f, h);
302                                 r[j] = z;
303                                 c = f / z;
304                                 s = h / z;
305                                 f = x * c + g * s;
306                                 g = g * c - x * s;
307                                 h = y * s;
308                                 y = y * c;
309                                 for (jj = 0; jj < n; jj++)
310                                 {
311                                         x = q[jj][j];
312                                         z = q[jj][i];
313                                         q[jj][j] = x * c + z * s;
314                                         q[jj][i] = z * c - x * s;
315                                 }
316                                 z = radius(f, h);
317                                 d[j] = z;       /* rotation can be arbitrary
318                                                  * id z=0 */
319                                 if (z)
320                                 {
321                                         z = 1.0 / z;
322                                         c = f * z;
323                                         s = h * z;
324                                 }
325                                 f = (c * g) + (s * y);
326                                 x = (c * y) - (s * g);
327                                 for (jj = 0; jj < m; jj++)
328                                 {
329                                         y = p[jj][j];
330                                         z = p[jj][i];
331                                         p[jj][j] = y * c + z * s;
332                                         p[jj][i] = z * c - y * s;
333                                 }
334                         }
335                         r[l] = 0.0;
336                         r[k] = f;
337                         d[k] = x;
338                 }
339         }
340         tvector_free(r, 0, double);
341 }
342 /*
343 void mat_svd(Mat *a, Mat *p, Vec *d, Mat *q)
344 {
345     Mat *at, *pt, *qt;
346     if(a->m >= a->n)
347     {
348         svd(a->m, a->n, a->el, p->el, d->el, q->el);
349         return;
350     }
351     at = mat_make(a->n, a->m);
352     pt = mat_make(a->m, a->m);
353     qt = mat_make(a->n, a->m);
354     mat_svd(at, qt, d, pt);
355     mat_transp(p, pt);
356     mat_transp(q, qt);
357     mat_free(at);
358     mat_free(pt);
359     mat_free(qt);
360 }
361 */
362 
363 void            mat_svd(Mat * a, Mat * p, Vec * d, Mat * q)
364 {
365         Mat            *at, *pt, *qt;
366         int             i, j;
367         if (a->m >= a->n)
368         {
369                 svd(a->m, a->n, a->el, p->el, d->el, q->el);
370                 return;
371         }
372         at = mat_make(a->n, a->m);
373         pt = mat_make(a->m, a->m);
374         qt = mat_make(a->n, a->m);
375         mat_transp(at, a);
376         mat_svd(at, qt, d, pt);
377         for (i = 0; i < a->m; i++)
378         {
379                 for (j = 0; j < a->m; j++)
380                         p->el[i][j] = pt->el[i][j];
381                 for (j = a->m; j < a->n; j++)
382                         p->el[i][j] = 0.0;
383         }
384         for (i = 0; i < a->n; i++)
385         {
386                 for (j = 0; j < a->m; j++)
387                         q->el[i][j] = qt->el[i][j];
388                 for (j = a->m; j < a->n; j++)
389                         q->el[i][j] = 0.0;
390         }
391         mat_free(at);
392         mat_free(pt);
393         mat_free(qt);
394 
395 }
396 

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