<![CDATA[
#define M 3
#define N 2
/* pseudoinverse */
int main() {
array double a[M][N] = {7, 8,
1, 3,
6, 4}; /* m-by-n matrix */
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.053595 -0.209150 0.264052
0.139869 0.228758 -0.201307
a*pinverse(a)*a =
7.000000 8.000000
1.000000 3.000000
6.000000 4.000000
pinverse(a)*a*pinverse(a) =
-0.053595 -0.209150 0.264052
0.139869 0.228758 -0.201307
]]>