int svd(int m,int n,int withu,int withv,double eps,double tol, double *a, double *q, double *u, double *v, double *vt) { 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<m;i++) { for (j=0;j<n;j++) u[i*n + j] = a[i*n + j]; } /* Householder's reduction to bidiagonal form. */ g = x = 0.0; for (i=0;i<n;i++) { e[i] = g; s = 0.0; l = i+1; for (j=i;j<m;j++) s += (u[j*n+i]*u[j*n+i]); if (s < tol) g = 0.0; else { f = u[i*n+i]; g = (f < 0) ? sqrt(s) : -sqrt(s); h = f * g - s; u[i*n+i] = f - g; for (j=l;j<n;j++) { s = 0.0; for (k=i;k<m;k++) s += (u[k*n+i] * u[k*n+j]); f = s / h; for (k=i;k<m;k++) u[k*n+j] += (f * u[k*n+i]); } /* end j */ } /* end s */ q[i] = g; s = 0.0; for (j=l;j<n;j++) s += (u[i*n+j] * u[i*n+j]); if (s < tol) g = 0.0; else { f = u[i*n+i+1]; g = (f < 0) ? sqrt(s) : -sqrt(s); h = f * g - s; u[i*n+i+1] = f - g; for (j=l;j<n;j++) e