Commit 9399bea2 authored by Peter Bruin's avatar Peter Bruin

simplify matsmall_direct_sum_projections

parent d3747a63
......@@ -368,38 +368,23 @@ matsmall_inverse_image (GEN A, GEN V, unsigned long p, GEN T) {
/*
Given a decomposition of a vector space as a direct sum of
subspaces, compute the projection maps onto all of the subspaces.
TODO: this assumes the subspaces are non-trivial (because of the
transpose).
*/
GEN
matsmall_direct_sum_projections (GEN spaces, unsigned long p, GEN T) {
pari_sp av = avma;
GEN A, A_inv_t, projections;
int i, j, dim = 0, count;
GEN A, A_inv, projections;
int i, j, dim;
/* Compute the dimension of the total space. */
for (i = 1; i < lg (spaces); i++)
dim += lg (gel (spaces, i)) - 1;
/*
Compute the change-of-bases matrix from the direct sum
of the subspaces to the ambient space.
*/
A = cgetg (dim + 1, t_MAT);
for (count = 0, i = 1; i < lg (spaces); i++)
for (j = 1; j < lg (gel (spaces, i)); j++)
gel (A, ++count) = gmael (spaces, i, j);
A_inv_t = matsmall_transpose (matsmall_inv (A, p, T), p, T);
projections = cgetg (lg (spaces), t_VEC);
for (count = 0, i = 1; i < lg (spaces); i++) {
/*
Select the correct range of rows from A^{-1}. We have
transposed that matrix because extracting columns is easier
than extracting rows.
*/
gel (projections, i) = cgetg (lg (gel (spaces, i)), t_MAT);
for (j = 1; j < lg (gel (spaces, i)); j++)
gmael (projections, i, j) = gel (A_inv_t, ++count);
gel (projections, i) = matsmall_transpose (gel (projections, i), p, T);
A = shallowconcat1(spaces);
A_inv = matsmall_inv(A, p, T);
projections = cgetg(lg(spaces), t_VEC);
for (i = j = 1; i < lg(spaces); i++, j += dim) {
dim = lg(gel(spaces, i)) - 1;
gel(projections, i) = rowslice(A_inv, j, j + dim - 1);
}
return gerepilecopy (av, projections);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment