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

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

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