/* svd.c -- Singular value decomposition. Translated to 'C' from the * original Algol code in "Handbook for Automatic Computation, * vol. II, Linear Algebra", Springer-Verlag. * * (C) 2000, C. Bond. All rights reserved. * * This is almost an exact translation from the original, except that * an iteration counter is added, to prevent stalls. This change is * described in the notes for the algorithm and corresponds to * similar changes in other translations. * * Returns an error code = 0, if no errors are encountered, and 'k' if * a failure to converge at the 'kth' singular value. * * NOTE: This version requires that the dimensions of matrix U (u) * be MxM rather than MxN. This is not required for the correct * operation of the algorithm, but permits the full orthogonal matrix U * to be returned to the caller. Otherwise, only the first N columns of * U are available. */ #include /* for array allocation of 'e' */ #include /* for 'fabs' */ int svd(int m,int n,int withu,int withv,double eps,double tol, double **a,double *q,double **u,double **v) { int i,j,k,l,l1,iter,retval; double c,f,g,h,s,x,y,z; double *e; e = (double *)calloc(n,sizeof(double)); retval = 0; /* Copy 'a' to 'u' */ for (i=0;i x) x = y; } /* end i */ /* accumulation of right-hand transformations */ if (withv) { for (i=n-1;i>=0;i--) { if (g != 0.0) { h = u[i][i+1] * g; for (j=l;j=0;i--) { l = i + 1; g = q[i]; for (j=l;j=0;k--) { iter = 0; test_f_splitting: for (l=k;l>=0;l--) { if (fabs(e[l]) <= eps) goto test_f_convergence; if (fabs(q[l-1]) <= eps) goto cancellation; } /* end l */ /* cancellation of e[l] if l > 0 */ cancellation: c = 0.0; s = 1.0; l1 = l - 1; for (i=l;i<=k;i++) { f = s * e[i]; e[i] *= c; if (fabs(f) <= eps) goto test_f_convergence; g = q[i]; h = q[i] = sqrt(f*f + g*g); c = g / h; s = -f / h; if (withu) { for (j=0;j 30) { retval = k; break; } x = q[l]; y = q[k-1]; g = e[k-1]; h = e[k]; f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y); g = sqrt(f*f + 1.0); f = ((x-z)*(x+z) + h*(y/((f<0)?(f-g):(f+g))-h))/x; /* next QR transformation */ c = s = 1.0; for (i=l+1;i<=k;i++) { g = e[i]; y = q[i]; h = s * g; g *= c; e[i-1] = z = sqrt(f*f+h*h); c = f / z; s = h / z; f = x * c + g * s; g = -x * s + g * c; h = y * s; y *= c; if (withv) { for (j=0;j