hello,
i'm writing a program in c which i wish to use the svdcmp method, but i'm not using the nr2 vector/matrix utilities (where they are created with indexing starting from 1). so i copy the c++ decompose() method which seems almost the same however i cannot get the correct results. from my output the svdcmp gives:

Decomposition matrices:
Matrix u
0.945570 -0.035867 -0.283998
0.202496 -0.125480 1.065605
0.347207 1.182848 -0.239236
Diagonal of matrix w
42.164364 19.190076 0.945570
Matrix v-transpose
-0.526086 0.730554 -0.435343
0.832413 0.547171 -0.087707
-0.174133 0.408527 0.895982

Check product against original matrix:
Original matrix:
23.588258 -6.362339 6.410960
-6.589769 28.441677 -12.996656
11.415051 -15.084915 11.724426
Product u*w*(v-transpose):
-21.500914 28.640419 -17.537102
-6.671670 5.331614 -2.603019
11.232515 23.022928 -8.566856

and from the original the svdcmp gives:
Matrix u
-0.465114 0.827668 -0.314059
0.713037 0.560516 0.421188
-0.524639 0.028036 0.850863
Diagonal of matrix w
42.168671 19.200581 2.778243
Matrix v-transpose
-0.513623 0.738779 -0.436344
0.841102 0.534006 -0.085935
-0.169523 0.411148 0.895667

Check product against original matrix:
Original matrix:
23.588301 -6.362340 6.410960
-6.589770 28.441700 -12.996700
11.415100 -15.084900 11.724400
Product u*w*(v-transpose):
23.588295 -6.362344 6.410962
-6.589768 28.441696 -12.996696
11.415096 -15.084901 11.724397

so as you can see they are different but i've checked my version with the decompose() method and it's word for word exact. could anyone offer any insight into what i'm doin wrong?
sincerely,
juan

void svd(double **a, int m, int n, double* w, double **v)
{
int flag, i, its, j, jj, k, l, nm;
double anorm, c, f, g, h, s, scale, x, y, z, *rv1;
rv1 = (double*) malloc(sizeof(double)*n);
if(rv1==NULL) {
printf("Could not malloc rv1 in svd()");
return;
}

g = scale = anorm = 0.0;
for (i = 0; i < n; ++i) {
l = i+2;
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 != 0.0) {
for (k = i; k < m; k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f = a[i][i];
g = -sign2(sqrt(s),f);
h = f*g-s;
a[i][i] = f - g;

for(j = l-1; 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+1 <= m) && (i+1 != n)) {
for(k = l-1; k < n; ++k) {
scale += fabs(a[i][k]);
}
if (scale != 0) {
for (k = l-1; k < n; ++k) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f = a[i][l-1];
g = -sign2(sqrt(s), f);
h = f*g - s;
a[i][l-1] = f - g;
for(k = l-1; k < n; ++k) {
rv1[k] = a[i][k]/h;
}
for(j = l-1; j < m; ++j) {
for(s = 0.0, k=l-1; k < n; ++k) {
s += a[j][k]*a[i][k];
}
for(k=l-1; k < n; ++k) {
a[j][k] += s*rv1[k];
}
}
for (k = l-1; k < n; ++k) {
a[i][k] *= scale;
}
}
}
anorm = max(anorm, (fabs(w[i])+fabs(rv1[i])));
}

for (i = n-1; i >= 0; i--) {
if(i < n-1) {
if (g != 0) {
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;
}

for(i=min(m, n)-1; i >= 0; --i) {
l = i + 1;
g = w[i];
for(j = l; j < n; ++j) {
a[i][j] = 0.0;
}
if(g != 0) {
g = 1.0/g;
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];
}
for(k = n-1; k >= 0; --k) {
for(its = 0; its < 30; ++its) {
flag = 1;
for(l = k; l >=0; --l) {
nm = l - 1;
if(l==0||(abs(rv1[l]) <= tol*anorm)) {
flag = 0;
break;
}
if(fabs(w[nm]) <= tol*anorm) {
break;
}
}
if(flag) {
c = 0.0;
s = 1.0;
for(i=l; i < k+1; ++i) {
f = s*rv1[i];
rv1[i] = c*rv1[i];
if(fabs(f) <= tol*anorm) {
break;
}
g = w[i];
h = pythag(f,g);
w[i] = h;
h = 1.0/h;
c = g*h;
s = -f*h;
for(j = 0; 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) {
if(z < 0.0) {
w[k] = -z;
for(j = 0; j < n; ++j) {
v[j][k] = -v[j][k];
}
}
break;
}
if(its==29) {
printf("no convergence in 30 iterations.\n");
}
x = w[l];
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+h)+h*((y/(f+sign2(g,f)))-h))/x;
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 *= c;
for(jj=0; 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=0; 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(rv1);
}

p.s. sign2==SIGN
max==MAX / min == MIN

hi,
so i cut/paste the SVD::decompose() method to my own routine and change it so it works on it's own in C, and i still get the different results. i just changed the NR types to the C type double/int etc, and set rv1 to double* rv1 = (double*)malloc(n*sizeof(double));
and the results are below. so they're different from above but still not the same as svdcmp(), which should have equivalent functionality as SVD::decompose, right?

Decomposition matrices:
Matrix u
0.920641 -0.216529 -0.285165
0.201668 -0.126781 1.065608
0.355188 1.180460 -0.239316
Diagonal of matrix w
42.168048 19.188613 0.920641
Matrix v-transpose
-0.513760 0.738837 -0.436086
0.840339 0.535781 -0.082272
-0.172861 0.408728 0.896136

Check product against original matrix:
Original matrix:
23.588258 -6.362339 6.410960
-6.589769 28.441677 -12.996656
11.415051 -15.084915 11.724426
Product u*w*(v-transpose):
-23.391141 26.349440 -16.823010
-6.582905 5.380581 -2.629153
11.378051 23.112126 -8.592543