#define M 2
#define N 3
/* 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);
}
rank(a) = 2
pinverse(a) =
0.128000 -0.104000
0.030769 0.061538
-0.142154 0.235692
a*pinverse(a)*a =
7.000000 8.000000 1.000000
3.000000 6.000000 4.000000
pinverse(a)*a*pinverse(a) =
0.128000 -0.104000
0.030769 0.061538
-0.142154 0.235692