<![CDATA[
#define M 3
#define N 3
/* pseudoinverse */
int main() {
/* singular matrix rank(a) == 2 */
array double a[M][N] = {1, 2, 3,
4, 5, 6,
7, 8, 9};
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) = 2
pinverse(a) =
-0.638889 -0.166667 0.305556
-0.055556 0.000000 0.055556
0.527778 0.166667 -0.194444
a*pinverse(a)*a =
1.000000 2.000000 3.000000
4.000000 5.000000 6.000000
7.000000 8.000000 9.000000
pinverse(a)*a*pinverse(a) =
-0.638889 -0.166667 0.305556
-0.055556 0.000000 0.055556
0.527778 0.166667 -0.194444
]]>