<![CDATA[
#define M 3
#define N 3
/* pseudoinverse */
int main() {
/* singular matrix rank(a) == 1 */
array double a[M][N] = {1, 2, 3,
1, 2, 3,
2, 4, 6};
int r;
array double pinv[N][M];
r = rank(a);
pinv = pinverse(a);
printf("rank(a) = %d\n", r);
printf("pinverse(a) = \n%f\n", pinv);
printf("a*pinverse(a)*a = \n%f\n", a*pinv*a);
printf("pinverse(a)*a*pinverse(a) = \n%f\n", pinv*a*pinv);
}
]]>
<![CDATA[
rank(a) = 1
pinverse(a) =
0.011905 0.011905 0.023810
0.023810 0.023810 0.047619
0.035714 0.035714 0.071429
a*pinverse(a)*a =
1.000000 2.000000 3.000000
1.000000 2.000000 3.000000
2.000000 4.000000 6.000000
pinverse(a)*a*pinverse(a) =
0.011905 0.011905 0.023810
0.023810 0.023810 0.047619
0.035714 0.035714 0.071429
]]>