<![CDATA[
#define M 3
#define N 3
/* pseudoinverse */
int main() {
/* non-singular matrix rank(a) == 3 */
array double a[M][N] = {1, 2, 3,
5, 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) = 3
pinverse(a) =
-0.500000 1.000000 -0.500000
-0.500000 -2.000000 1.500000
0.833333 1.000000 -0.833333
a*pinverse(a)*a =
1.000000 2.000000 3.000000
5.000000 5.000000 6.000000
7.000000 8.000000 9.000000
pinverse(a)*a*pinverse(a) =
-0.500000 1.000000 -0.500000
-0.500000 -2.000000 1.500000
0.833333 1.000000 -0.833333
]]>