<![CDATA[
int main() {
array double a[3][3] = {7, 8, 1,
3, 6, 4,
3, 5, 7}; /* m-by-n matrix */
array float b[3][3] = {7, 8, 1,
3, 6, 4,
3, 5, 7}; /* m-by-n matrix */
int m = 3, n = 3;
int lmn = min(m,n);
array double s[lmn];
array double u[m][m];
array double vt[n][n];
array float fs[lmn];
array float fu[m][m];
array float fvt[n][n];
int status;
svd(a, s, NULL, NULL);
printf("s =\n%f\n", s);
svd(a, s, u, vt);
printf("u =\n%f\n", u);
printf("s =\n%f\n", s);
printf("vt =\n%f\n", vt);
printf("u*s*v^t =\n%f\n", u*diagonalmatrix(s)*transpose(vt));
status = svd(b, fs, fu, fvt);
if(status ==0) {
printf("fu =\n%f\n", fu);
printf("fs =\n%f\n", fs);
printf("fvt =\n%f\n", fvt);
printf("u*s*v^t =\n%f\n", fu*diagonalmatrix(fs)*transpose(fvt));
}
else
printf("error: numerical error in svd()\n");
}
]]>
<![CDATA[
s =
15.071167 5.471953 0.957939
u =
-0.660443 0.704551 0.259659
-0.511773 -0.169318 -0.842271
-0.549458 -0.689158 0.472395
s =
15.071167 5.471953 0.957939
vt =
-0.517995 0.430638 0.739075
-0.736603 0.214678 -0.641349
-0.434853 -0.876621 0.206007
u*s*v^t =
7.000000 8.000000 1.000000
3.000000 6.000000 4.000000
3.000000 5.000000 7.000000
fu =
-0.660443 0.704551 0.259659
-0.511773 -0.169318 -0.842271
-0.549458 -0.689158 0.472395
fs =
15.071167 5.471953 0.957939
fvt =
-0.517995 0.430638 0.739075
-0.736603 0.214678 -0.641349
-0.434853 -0.876621 0.206007
u*s*v^t =
7.000000 8.000000 1.000000
3.000000 6.000000 4.000000
3.000000 5.000000 7.000000
]]>