struct SVD { Int m,n; MatDoub u,v; VecDoub w; Doub eps, tsh; SVD(MatDoub& a) : m(a.nrows()), n(a.ncols()), u(a), v(n,n), w(n) { eps = numeric_limits::epsilon(); decompose(); reorder(); tsh = 0.5*sqrt(m+n+1.)*w[0]*eps; } void solve(VecDoub_I &b, VecDoub_O &x, Doub thresh); void solve(MatDoub_I &b, MatDoub_O &x, Doub thresh); Int rank(Doub thresh); Int nullity(Doub thresh); MatDoub range(Doub thresh); MatDoub nullspace(Doub thresh); Doub inv_condition() { return (w[0] <= 0. || w[n-1] <= 0.) ? 0. : w[n-1]/w[0]; } void decompose(); void reorder(); Doub pythag(const Doub a, const Doub b); }; void SVD::solve(VecDoub_I &b, VecDoub_O &x, Doub thresh = -1.) { Int i,j,jj; Doub s; if (b.size() != m || x.size() != n) throw("SVD::solve bad sizes"); VecDoub tmp(n); tsh = (thresh >= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps); for (j=0;j tsh) { for (i=0;i= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps); for (j=0;j tsh) nr++; return nr; } Int SVD::nullity(Doub thresh = -1.) { Int j,nn=0; tsh = (thresh >= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps); for (j=0;j tsh) { for (i=0;i=0;i--) { if (i < n-1) { if (g != 0.0) { for (j=l;j=0;i--) { l=i+1; g=w[i]; for (j=l;j=0;k--) { for (its=0;its<30;its++) { flag=true; for (l=k;l>=0;l--) { nm=l-1; if (l == 0 || abs(rv1[l]) <= eps*anorm) { flag=false; break; } if (abs(w[nm]) <= eps*anorm) break; } if (flag) { c=0.0; s=1.0; for (i=l;i 1); for (k=0;k (m+n)/2) { for (i=0;i absb ? absa*sqrt(1.0+SQR(absb/absa)) : (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)))); }