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

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

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