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

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

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