home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
vis-ftp.cs.umass.edu
/
vis-ftp.cs.umass.edu.tar
/
vis-ftp.cs.umass.edu
/
pub
/
Software
/
ASCENDER
/
ascender.tar.Z
/
ascender.tar
/
Triangulate
/
svalue.c
< prev
next >
Wrap
C/C++ Source or Header
|
1995-04-13
|
12KB
|
547 lines
/* @(#)svalue.c 1.10 12/20/94 */
/*>>
Singular value decomposition routines taken from the book
" Numerical Recipes in C " (Cambridge University Press 1988)
<<*/
#include "cvar.h"
#include <math.h>
#include "rh_util.h"
#include "alloc.h"
#include "array_utils.s"
#include "error_message.s"
/*>>
Singular value decomposition routines.
<<*/
#define NUMITER 30
FUNCTION_DEF ( void svbksb, (
double **u,
double w[],
double **v,
int m,
int n,
double b[],
double x[]))
/* Solves A.X = B for a vector X, where A is specified by the arrays
u[1..m][1..n], w[1..n], v[1..n][1..n] as returned by svdcmp.
m and n are dimensions of A, and will be equal for square matrices.
b[1..m] is the input right-hand side.
x[1..n] is the output solution vector.
No input quantities are destroyed, so the routine may be called
sequentially with different b's.
*/
{
int jj, j, i;
double s, *tmp;
tmp = VECTOR(1,n,double);
for (j=1; j<=n; j++)
{
/* Calculate U^T.B */
s = 0.0;
if (w[j]) /* Nonzero result only if w[j] is non-zero. */
{
for (i=1; i<=m; i++)
s += u[i][j] * b[i];
s /= w[j];
}
tmp[j] = s;
}
/* Matrix multiply by V to get answer */
for (j=1; j<=n; j++)
{
s = 0.0;
for (jj=1; jj<=n; jj++)
s += v[j][jj] * tmp[jj];
x[j] = s;
}
/* Free temporaries */
FREE (tmp, 1);
}
static double at, bt, ct;
/* Pythag computes (sqrt(a*a + b*b)) without destructive over or under flow */
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
(ct=bt/at, at*sqrt(1.0+ct*ct)) : \
(bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)):0.0))
static double maxarg1, maxarg2;
#define FloatMAX(a,b) (maxarg1=(a), maxarg2=(b), (maxarg1) > (maxarg2) ? \
(maxarg1) : (maxarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
FUNCTION_DEF ( int svdcmp, (
double **a,
int m,
int n,
double *w,
double **v))
/* Given a matrix a[1..m][1..n], this routine computes its singular value
decomposition, A = U.W.V^T. The matrix U replaces a on output.
The diagonal matrix of singular values is output as a vector w[1..n].
The matrix V (not the transpose V^T) is output as v[1..n][1..n].
m must be greater or equal to n; If it is smaller, then a should
be filled up to square with zero rows.
The routine returns 1 on success, 0 on failure. In case of failure,
a good thing to do is to permute the rows and try again.
*/
{
int flag, i, its, j, jj, k, l, nm;
double c, f, h, s, x, y, z;
double anorm = 0.0, g = 0.0, scale = 0.0;
double *rv1;
int returnval = 1;
/* Check the bounds of the array */
if (m < n)
{
error_message ("SVDCMP : matrix has fewer rows (%d) than columns (%d)",
m, n);
bail_out (1);
}
/* Take a temporary array */
rv1 = VECTOR(1,n,double);
/* Householder reduction to bidiagonal form */
for (i=1; i<=n; i++)
{
l = i+1;
rv1[i] = scale*g;
g = s = scale = 0.0;
if (i <= m)
{
for (k=i; k<=m; k++) scale += fabs(a[k][i]);
if (scale)
{
for (k=i; k<=m; k++)
{
a[k][i] /= scale;
s += a[k][i] * a[k][i];
}
f = a[i][i];
g = - SIGN(sqrt(s), f);
h = f*g-s;
a[i][i] = f-g;
if (i != n)
{
for (j=l; j<=n; j++)
{
for (s=0.0, k=i; k<=m; k++) s += a[k][i] * a[k][j];
f = s/h;
for (k=i; k<=m; k++) a[k][j] += f*a[k][i];
}
}
for (k=i; k<=m; k++) a[k][i] *= scale;
}
}
w[i] = scale*g;
g = s = scale = 0.0;
if (i<=m && i != n)
{
for (k=l; k<=n;k++) scale += fabs (a[i][k]);
if (scale)
{
for (k=l; k<=n; k++)
{
a[i][k] /= scale;
s += a[i][k] * a[i][k];
}
f = a[i][l];
g = -SIGN(sqrt(s),f);
h = f*g-s;
a[i][l] = f-g;
for (k=l; k<=n;k++)
rv1[k] = a[i][k]/h;
if (i != m)
{
for (j=l; j<=m; j++)
{
for (s=0.0, k=l; k<=n; k++) s += a[j][k] * a[i][k];
for (k=l; k<=n;k++) a[j][k] += s*rv1[k];
}
}
for (k=l; k<=n; k++) a[i][k] *= scale;
}
}
anorm = FloatMAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
}
/* Accumulation of right-hand transformations */
for (i=n; i>=1; i--)
{
if (i<n)
{
if (g)
{
for (j=l; j<=n; j++)
v[j][i] = (a[i][j] / a[i][l]) / g;
for (j=l; j<=n; j++)
{
for (s=0.0, k=l; k<=n; k++) s += a[i][k] * v[k][j];
for (k=l; k<=n; k++) v[k][j] += s*v[k][i];
}
}
for (j=l; j<=n; j++) v[i][j] = v[j][i] = 0.0;
}
v[i][i] = 1.0;
g = rv1[i];
l = i;
}
/* Accumulation of left-hand transformations. */
for (i=n; i>=1; i--)
{
l = i+1;
g = w[i];
if (i < n )
for (j=l; j<=n; j++) a[i][j]= 0.0;
if (g)
{
g = 1.0/g;
if (i != n)
{
for (j=l; j<=n; j++)
{
for (s = 0.0, k=l; k<=m; k++) s += a[k][i]*a[k][j];
f = (s/a[i][i]) * g;
for (k=i; k<=m; k++) a[k][j] += f * a[k][i];
}
}
for (j=i; j<=m; j++) a[j][i] *= g;
}
else
{
for (j=i; j<=m; j++) a[j][i] = 0.0;
}
++a[i][i];
}
/* diagonalization of the bidiagonal form */
for (k=n; k>=1; k--) /* Loop over singular values */
{
for (its = 1; its <=NUMITER; its++) /* Loop over allowed iterations */
{
flag = 1;
for (l=k; l>=1; l--) /* Test for splitting */
{
double rv_plus_anorm, m_plus_anorm;
nm = l-1; /* Note that rv1[1] is always zero */
/* Bug fix, RIH : Change of the convergence criteria */
rv_plus_anorm = (double)(fabs(rv1[l]) + anorm);
if ((rv_plus_anorm == anorm)
|| ((its==NUMITER) && (rv_plus_anorm <= 1.00001 * anorm)))
{
flag = 0;
break;
}
m_plus_anorm = (double)(fabs(w[nm]) + anorm);
if ((m_plus_anorm == anorm) ||
((its==NUMITER) && (m_plus_anorm <= 1.00001 * anorm)))
break; /* Bug fix, RIH */
}
if (flag)
{
c = 0.0;
s = 1.0;
for (i=l; i<=k; i++)
{
f = s*rv1[i];
if ((double)(fabs(f) + anorm) != anorm) /* Bug fix, RIH */
{
g = w[i];
h = PYTHAG(f, g);
w[i] = h;
h = 1.0/h;
c = g*h;
s = (-f*h);
for (j=1; j <=m; j++)
{
y = a[j][nm];
z = a[j][i];
a[j][nm] = y*c + z*s;
a[j][i] = z*c - y*s;
}
}
}
}
z = w[k];
if ( l == k)
{
/* Convergence */
if (z < 0.0)
{
/* Singular value is made non-negative */
w[k] = -z;
for (j = 1; j<=n; j++) v[j][k] = (-v[j][k]);
}
break;
}
if (its == NUMITER)
{
returnval = 0;
goto cleanup;
}
x = w[l]; /* Shift from bottom 2-by-2 minor */
nm = k-1;
y = w[nm];
g = rv1[nm];
h = rv1[k];
f = ((y-z)*(y+z) + (g-h)*(g+h))/(2.0*h*y);
g = PYTHAG(f, 1.0);
f = ((x-z)*(x+z) + h*((y/(f+SIGN(g,f)))-h))/x;
/* Next QR transformation */
c = s = 1.0;
for (j=l; j<=nm; j++)
{
i = j+1;
g = rv1[i];
y = w[i];
h = s*g;
g = c*g;
z = PYTHAG(f,h);
rv1[j] = z;
c = f/z;
s = h/z;
f = x*c + g*s;
g = g*c - x*s;
h = y*s;
y = y*c;
for (jj = 1; jj <= n; jj++)
{
x = v[jj][j];
z = v[jj][i];
v[jj][j] = x*c + z*s;
v[jj][i] = z*c - x*s;
}
z = PYTHAG(f,h);
w[j] = z;
if (z)
{
z = 1.0 / z;
c = f*z;
s = h*z;
}
f = (c*g) + (s*y);
x = (c*y) - (s*g);
for (jj=1; jj<=m; jj++)
{
y = a[jj][j];
z = a[jj][i];
a[jj][j] = y*c + z*s;
a[jj][i] = z*c - y*s;
}
}
rv1[l] = 0.0;
rv1[k] = f;
w[k] = x;
}
}
/* Free the temporaries */
cleanup:
FREE (rv1, 1);
/* Return success */
return (returnval);
}
FUNCTION_DEF ( int solve, (
double **a,
double *x,
double *b,
int nrow,
int ncol))
{
/* Find the least-squares solution to an over-constrained linear system.
* The inputs are
* a[1..nrow][1..ncol]
* b[1..nrow]
* The output is
* x[1..ncol]
*
* The routine returns 1 on success, 0 on failure. It fails if the SVD
* fails.
*/
int i, j;
double *w, **v, **u, *w2;
double besterror = HUGE;
double bestthreshold = 0.0;
int returnval = 1;
/* Allocate matrices */
u = MATRIX (1, nrow, 1, ncol, double);
w = VECTOR (1, ncol, double);
w2 = VECTOR (1, ncol, double);
v = MATRIX (1, ncol, 1, ncol, double);
/* Copy a to u */
for (i=1; i<=nrow; i++)
for (j=1; j<=ncol; j++)
u[i][j] = a[i][j];
/* Now do the singular value decomposition */
if (! svdcmp(u, nrow, ncol, w, v))
{
returnval = 0;
goto cleanup;
}
/* Go into a loop of back-substitution, to find the best */
for (j=1; j<=ncol; j++)
{
/* Threshold on the i-th singular value */
double wmin = rhAbs(w[j]);
/* Set near-zero values to zero */
for (i=1; i<=ncol; i++)
{
if ((rhAbs(w[i]) < wmin)) w2[i] = 0.0;
else w2[i] = w[i];
}
/* Back substitute */
svbksb (u, w2, v, nrow, ncol, b, x);
/* Calculate the difference */
{
int i, k;
double error = 0.0;
for (i=1; i<=nrow; i++)
{
double bt = 0.0;
for (k=1; k<=ncol; k++) bt += a[i][k] * x[k];
error += (bt-b[i]) * (bt-b[i]);
}
/* If this is better than the best error, then record */
if (error <= besterror)
{
besterror = error;
bestthreshold = wmin;
}
}
}
/* Repeat the calculation with the best values */
{
double wmin = bestthreshold;
/* Set near-zero values to zero */
for (i=1; i<=ncol; i++)
{
if ((rhAbs(w[i]) < wmin)) w2[i] = 0.0;
else w2[i] = w[i];
}
/* Back substitute */
svbksb (u, w2, v, nrow, ncol, b, x);
}
/* Free the temporaries */
cleanup:
FREE (u, 1);
FREE (w, 1);
FREE (w2, 1);
FREE (v, 1);
/* return success */
return (returnval);
}
#define NEAR_ZERO 1.0e-12
FUNCTION_DEF ( int solve_simple, (
double **a,
double *x,
double *b,
int nrow,
int ncol))
{
/* Find the least-squares solution to an over-constrained linear system.
* The inputs are
* a[1..nrow][1..ncol]
* b[1..nrow]
* The output is
* x[1..ncol]
* It is allowable for x and b to be the same vector.
*
* The routine returns 1 on success, 0 on failure. It fails if the SVD
* fails.
*/
int i, j;
double *w, **v, **u, *w2;
double wmax;
int returnval = 1;
/* Allocate matrices */
u = MATRIX (1, nrow, 1, ncol, double);
w = VECTOR (1, ncol, double);
w2 = VECTOR (1, ncol, double);
v = MATRIX (1, ncol, 1, ncol, double);
/* Copy a to u */
for (i=1; i<=nrow; i++)
for (j=1; j<=ncol; j++)
u[i][j] = a[i][j];
/* Now do the singular value decomposition */
if (! svdcmp(u, nrow, ncol, w, v))
{
returnval = 0;
goto cleanup;
}
/* Find the maximum singular value */
wmax = -1.0;
for (j=1; j<=ncol; j++)
{
/* Threshold on the i-th singular value */
double wcurr = rhAbs(w[j]);
if (wcurr > wmax) wmax = wcurr;
}
/* Repeat the calculation with the best values */
{
double wmin = wmax * NEAR_ZERO;
/* Set near-zero values to zero */
for (i=1; i<=ncol; i++)
{
if ((rhAbs(w[i]) < wmin)) w2[i] = 0.0;
else w2[i] = w[i];
}
/* Back substitute */
svbksb (u, w2, v, nrow, ncol, b, x);
}
/* Free the temporaries */
cleanup:
FREE (u, 1);
FREE (w, 1);
FREE (w2, 1);
FREE (v, 1);
/* return success */
return (returnval);
}