Let E = √D.
ECE is symmetric.
We use this to find x and λ such that ECEx = λx.
Multiply on left by E.
EECEx = λEx
DCEx = λEx
(DC)(Ex) = λEx
Thus Ex is an eigenvector of DC and λ is its eigenvalue!
Algol 68 implementation
Here is some C code where
void lev(int n, R C[n][n], R d[n], R eval[n], R g[n][n]);
is the prototype where d[j] is Djj. The n eigenvectors are placed in g and the eigenvalues go to eval. The two blocks that follow “if(1)” are respectively optional checks for
Routine void evt(int n, R m[n][n], R val[n], R vec[n][n])
checks that eigenvectors have been found.
It has the same arguments as eigen except that all the arguments are inputs to evt.
Index j chooses which eigenvector to test.
Index y chooses which component of the eigenvector to compute.
Vector u is DCEx