<![CDATA[
int main() {
array double a[2][3] = {7, 8, 1,
3, 6, 4}; /* m-by-n matrix */
int m = 2, n = 3, i;
int lmn = min(m,n);
array double s[lmn], sm[m][n];
array double u[m][m];
array double vt[n][n];
svd(a, s, u, vt);
printf("u =\n%f\n", u);
printf("s =\n%f\n", s);
printf("vt =\n%f\n", vt);
for(i=0; i<lmn; i++)
sm[i][i] = s[i];
printf("u*s*v^t =\n%f\n", u*sm*transpose(vt));
}
]]>
<![CDATA[
u =
-0.818910 -0.573922
-0.573922 0.818910
s =
12.851503 3.136698
vt =
-0.580020 -0.497570 0.644981
-0.777715 0.102681 -0.620174
-0.242353 0.861325 0.446525
u*s*v^t =
7.000000 8.000000 1.000000
3.000000 6.000000 4.000000
]]>