svdcmp[A,m,n] does a singular value decomposition of the m by n matrix A. It returns a list {U,W,V} where U is a m by n matrix which is orthogonal if m=n, W is a diagonal n by n matrix containing the singular values of A V is a n by n orthogonal matrix These matrices satisfy A = U.W.V^t. **************************************************************************** svdcmp:=proc(aa, m, n) local flag,i,its,j,jj,k,l,nm,anorm,c,f,g,h,s,scale,x,y,z,rv1,a,w,v,ww; a:= Array(1..m, 1..n); for i from 1 to m do for j from 1 to n do a[i,j]:=aa[i,j]; end do; end do; v:= Array(1..n, 1..n); w:=Array(1..n); rv1:=Array(1..n); g:=0.0; scale:=0.0; anorm:=0.0; for i from 1 to n do l:=i+1; rv1[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0; if (i <= m) then for k from i to m do scale := scale+abs(a[k,i]); end do; if (scale<>0.0) then for k from i to m do a[k,i] := a[k,i]/scale; s := s+a[k,i]*a[k,i]; end do; f:=a[i,i]; g := -mysign(sqrt(s),f); h:=f*g-s; a[i,i]:=f-g; for j from l to n do s:=0.0; for k from i to m do s := s+a[k,i]*a[k,j]; end do; f:=s/h; for k from i to m do a[k,j] := a[k,j]+f*a[k,i]; end do; end do; for k from i to m do a[k,i] := a[k,i]*scale; end do; end if; end if; w[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0; if (i <= m and i <> n) then for k from l to n do scale := scale+abs(a[i,k]); end do; if (scale<>0.0) then for k from l to n do a[i,k] := a[i,k]/scale; s := s+a[i,k]*a[i,k]; end do; f:=a[i,l]; g := -mysign(sqrt(s),f); h:=f*g-s; a[i,l]:=f-g; for k from l to n do rv1[k]:=a[i,k]/h; end do; for j from l to m do s:=0.0; for k from l to n do s := s+a[j,k]*a[i,k]; end do; for k from l to n do a[j,k] := a[j,k]+s*rv1[k]; end do; end do; for k from l to n do a[i,k] := a[i,k]*scale; end do; end if; end if; anorm:=fmax(anorm,(abs(w[i])+abs(rv1[i]))); end do; for i from n by -1 to 1 do if (i < n) then if (g<>0.0) then for j from l to n do v[j,i]:=(a[i,j]/a[i,l])/g; end do; for j from l to n do s:=0.0; for k from l to n do s := s+a[i,k]*v[k,j]; end do; for k from l to n do v[k,j] := v[k,j]+s*v[k,i]; end do; end do; end if; for j from l to n do v[i,j]:=0.0; v[j,i]:=0.0; end do; end if; v[i,i]:=1.0; g:=rv1[i]; l:=i; end do; for i from imin(m,n) by -1 to 1 do l:=i+1; g:=w[i]; for j from l to n do a[i,j]:=0.0; end do; if (g<>0.0) then g:=1.0/g; for j from l to n do s:=0.0; for k from l to m do s := s+a[k,i]*a[k,j]; end do; f:=(s/a[i,i])*g; for k from i to m do a[k,j] := a[k,j]+f*a[k,i]; end do; end do; for j from i to m do a[j,i] := a[j,i]*g; end do; else for j from i to m do a[j,i]=0.0; end do; end if; a[i,i]:=a[i,i]+1; end do; for k from n by -1 to 1 do for its from 1 to 30 do flag:=1; for l from k by -1 to 1 do nm:=l-1; if ((abs(rv1[l])+anorm) = anorm) then flag:=0; break; end if; if ((abs(w[nm])+anorm) = anorm) then break; end if; end do; if (flag<>0) then c:=0.0; s:=1.0; for i from l to k do f:=s*rv1[i]; rv1[i]:=c*rv1[i]; if ((abs(f)+anorm) = anorm) then break; end if; g:=w[i]; h:=pythag(f,g); w[i]:=h; h:=1.0/h; c:=g*h; s := -f*h; for j from 1 to m do y:=a[j,nm]; z:=a[j,i]; a[j,nm]:=y*c+z*s; a[j,i]:=z*c-y*s; end do; end do; end if; z:=w[k]; if (l = k) then if (z < 0.0) then w[k] := -z; for j from 1 to n do v[j,k] := -v[j,k]; end do; end if; break; end if; if (its = 30) then print("no convergence in 30 svdcmp iterations"); return(0); end if; 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+z)+h*((y/(f+mysign(g,f)))-h))/x; c:=1.0; s:=1.0; for j from l to nm do 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 := y*c; for jj from 1 to n do x:=v[jj,j]; z:=v[jj,i]; v[jj,j]:=x*c+z*s; v[jj,i]:=z*c-x*s; end do; z:=pythag(f,h); w[j]:=z; if (z<>0.0) then z:=1.0/z; c:=f*z; s:=h*z; end if; f:=c*g+s*y; x:=c*y-s*g; for jj from 1 to m do y:=a[jj,j]; z:=a[jj,i]; a[jj,j]:=y*c+z*s; a[jj,i]:=z*c-y*s; end do; end do; rv1[l]:=0.0; rv1[k]:=f; w[k]:=x; end do; end do; ww:=Array(1..n,1..n); for i from 1 to n do ww[i,i]:=w[i]; end do; return([a,ww,v]); end: pythag:=proc(a, b) local absa, absb; absa := abs(a); absb := abs(b); if (absa > absb) then return(absa*sqrt(1.0 + (absb/absa)*(absb/absa))); end if; if (absb = 0.0) then return(0.0); else return(absb*sqrt(1.0 + (absa/absb)*(absa/absb))); end if; end: mysign:=proc(a, b) if (b >= 0.0) then return(abs(a)); else return(-abs(a)); end if; end: fmax:=proc(a, b) if (a > b) then return(a); else return(b); end if; end: imin:=proc(a, b) if (a > b) then return(b); else return(a); end if; end: