Hi,

I tried to write svdcmp in php but it doesnt work properly. The numerical results are OK, but the program swaps the "coordinates" of the entries. For example, for the matrix
[1,2,-1]
[2.3,4,4]
[-2,5.1,1]
the resulting w vector is [2.006, 3.779, 7.483] instead of [7.483, 3.779, 2.006] and the other resulting matrices are confused as well. Can anyone figure out what might be the problem?

Here is "my" function:

function svdcmp ($A) {
$m = sizeof($A)-1;
$n = sizeof($A[1])-1;
// The 0th element of the arrays is always 0
$g=$scale=$anorm=0;
$l=$nm=0;
for ($i=1;$i<=$n;$i++) {
$l = $i + 1;
$rv1[$i] = $scale*$g;
$g=$s=$scale=0;
if ($i<=$m) {
for($k=$i;$k<=$m;$k++) $scale += abs($A[$k][$i]);
if ($scale!=0) {
for ($k=$i;$k<=$m;$k++) {
$A[$k][$i] /= $scale;
$s += $A[$k][$i]*$A[$k][$i];
}
$f = $A[$i][$i];
$g = -1*($f>=0 ? abs(sqrt($s)) : -1*abs(sqrt($s)));
$h = $f*$g-$s;
$A[$i][$i] = $f - $g;
for ($j=$l;$j<=$n;$j++) {
$s = 0;
for ($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;
if ($i<=$m and $i!=$n) {
for ($k=$l;$k<=$n;$k++) $scale += abs($A[$i][$k]);
if ($scale!=0) {
for($k=$l;$k<=$n;$k++) {
$A[$i][$k] /= $scale;
$s += $A[$i][$k]*$A[$i][$k];
}
$f = $A[$i][$l];
$g = -1*($f>=0 ? abs(sqrt($s)) : -1*abs(sqrt($s)));
$h = $f*$g-$s;
$A[$i][$l] = $f - $g;
for($k=$l;$k<=$n;$k++) $rv1[$k] = $A[$i][$k] / $h;
for($j=$l;$j<=$m;$j++) {
$s = 0;
for($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 = max($anorm, (abs($w[$i]) + abs($rv1[$i])));
}
for($i=$n;$i>=1;$i--) {
if ($i<$n) {
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++) {
$s = 0;
for($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;
}
$v[$i][$i] = 1;
$g = $rv1[$i];
$l = $i;
}
for($i=min($m,$n);$i>=1;$i--) {
$l = $i + 1;
$g = $w[$i];
for($j=$l;$j<=$n;$j++) $A[$i][$j] = 0;
if ($g!=0) {
$g = 1 / $g;
for($j=$l;$j<=$n;$j++) {
$s = 0;
for($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;
}
$A[$i][$i] = $A[$i][$i] + 1;
}
for($k=$n;$k>=1;$k--) {
for($its=1;$its<=30;$its++) {
$flag = 1;
for($l=$k;$l>=1;$l--) {
$nm = $l - 1;
if ((abs($rv1[$l])+$anorm)==$anorm) {
$flag = 0;
break;
}
if ((abs($w[$nm])+$anorm)==$anorm) break;
}
if ($flag==1) {
$c = 0;
$s = 1;
for($i=$l;$i<=$k;$i++) {
$f = $s*$rv1[$i];
$rv1[$i] = $c*$rv1[$i];
if ((abs($f)+$anorm)==$anorm) break;
$g = $w[$i];
$h = pythag($f,$g);
$w[$i] = $h;
$h = 1 / $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) {
if ($z<0) {
$w[$k] = -$z;
for($j=1;$j<=$n;$j++) $v[$j][$k] = -$v[$j][$k];
}
break;
}
if ($its==30) die ('No convergence in 30 svdcmp iterations');
$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*$h*$y);
$g = pythag($f,1);
$signhez = ($f>=0 ? abs($g) : -1*abs($g));
$f = (($x-$z)*($x+$z)+$h*(($y/($f+$signhez))-$h))/$x;
$c = $s = 1;
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=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!=0) {
$z = 1 / $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;
$rv1[$k] = $f;
$w[$k] = $x;
}
}
$back['U'] = $A;
$back['w'] = $w;
$back['V'] = $v;
return $back;
}