<![CDATA[
#define M 3
#define N 3
int main() {
/* singular matrix rank(a) == 3 */
array double a[M][N] = {1, 2, 3,
5, 5, 6,
7, 8, 9};
int r = rank(a);
array double orth[M][r];
orthonormalbase(a, orth);
printf("rank(a) = %d\n", r);
printf("orth from orthonormalbase(a, orth) = \n%f\n", orth);
printf("transpose(orth)*orth = \n%f\n", transpose(orth)*orth);
}
]]>
<![CDATA[
rank(a) = 3
orth from orthonormalbase(a, orth) =
-0.210138 0.955081 -0.208954
-0.541507 -0.291649 -0.788486
-0.814010 -0.052541 0.578470
transpose(orth)*orth =
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
]]>