Generating the partners in a multi-dimensional irreducible representation. I am trying to block diagonalize a Hermitian matrix using the irreducible representations of its symmetry group.  
Using the group's character table, it is straightforward to generate a set of projection matrices that will project an arbitrary vector onto the irreducible sub-space, as follows
$$\mathbf{P}^{(\Gamma_n)} = \frac{l_n}{h} \sum_R \chi^{(\Gamma_n)}(R) P_R$$
where $\Gamma_n$ is the representation, $l_n$ is its dimensionality, $h$ is the size of the group, $\chi^{(\Gamma_n)}(R)$ is the character for the $R^\text{th}$ operation in the representation, and $P_R$ is the $R^\text{th}$ matrix operator. For example, looking at the 2-dimensional $E$ representation of the symmetry of the square in $\mathbb{R}^3$ gives 
$$
\mathbf{P}^E = \left(
\begin{array}{ccc}
 1 & 0 & 0 \\
 0 & 1 & 0 \\
 0 & 0 & 0 \\
\end{array}
\right)
$$
For this simple case, it is apparent that given $(\begin{array}{ccc}1&0&0\end{array})$ it's partner in the representation is $(\begin{array}{ccc}0&1&0\end{array})$. But, for cases where you have multiple representations of the same type, it is not clear how to generate the basis for each one as $\mathbf{P}^{(\Gamma_n)}$ projects onto all of the representations of $\Gamma_n$ present. For instance, consider the case of the symmetry of the triangle where d-orbitals (quadratic functions) form two $E$ representations.
Both Tinkham and Dresselhaus, et al., give the formula
$$\mathrm{P}^{(\Gamma_n)}_{kl} = \frac{l_n}{h} \sum_R D^{(\Gamma_n)}(R)^\ast_{kl}P_R$$
where 
$$D^{(\Gamma_n)}(R)^\ast_{kl} = \mathbf{e}^{(\Gamma_n)}_k P_R \mathbf{e}^{(\Gamma_n)\ast}_l$$ 
given that $\mathbf{e}^{(\Gamma_n)}_k$ is the $k^\text{th}$ vector in the irreducible representation $\Gamma_n$ which is related to $\mathbf{P}^{(\Gamma_n)}$ via
$$\mathbf{P}^{(\Gamma_n)} = \sum_k \mathrm{P}^{(\Gamma_n)}_{kk}. $$
When $\mathrm{P}^{(\Gamma_n)}_{kl}$ is applied to $k^\text{th}$ vector gives the $l^\text{th}$ vector in the representation, in other words, given a vector in the representation you can generate all its partners. 
My question is how do I generate $D^{(\Gamma_n)}(R)^\ast_{kl}$ when I only know $\mathbf{P}^{(\Gamma_n)}$? This seems like a chicken and egg problem. 
 A: Better late than never!
This sounds more like an algorithm problem than an actual math problem, so I'll try to answer it as such.
In addition to P(Γn) you also have the projection operators of any subgroups.
If you lower the symmetry the degenerate irreducible representations will decompose into 1-dimensional ones and the partner functions will belong to different symmetry species in the lower symmetry group.
E.g. in your triangle example functions in the E representations of D3h will decompose into eigenspaces in A' and A'' in one of the Cs subgroups containing a σv. Project your E functions into A' and A'', choose a function in A' and apply the symmetry operations of D3h. This will result in functions which have some (can be zero) component in a 1-dimensional subspace of the A'' eigenspace. This is your partner function.
If you're a programmer I've written some code on the dev branch of libmsym which does this, in addition to using direct product information for molecular orbitals (not finished, so it's ugly). I'll be posting a link to my thesis, with a more detailed description, there too when it's ready.
Note: If you're just looking at the spherical harmonics atomic orbitals the partner functions for E representations will be ±m
