#
# Program for singular value decomposition
#
#    A = U*S*transpose(V)
#
SVD:=proc(  A::matrix,  # input matrix
            U::evaln,   # output
            S::evaln,   # output
            V::evaln    # output
          )
   local m, n, s, i;

   m:=linalg[rowdim](A);
   n:=linalg[coldim](A);
   s:=evalf(Svd(A, U, V));

   S:=matrix(m,n, [seq(0,i=1..m*n)]);
   for i to min(m,n) do
      S[i,i]:=s[i]
   od;

   print(` `)

end;