Debugging

Copy the function fdtest() from Appendix B into the file ``fdtest.m''. This program takes an initial $W$ containing Ns wave functions and computes the energy and gradient for that $W$ using your functions getE() and getgrad(). The function then forms a random direction and steps different distances along that direction, printing at each step the ratio of the actual change in energy to the change in energy expected from your gradient according to the formula Eq. (5). The number printed immediately below this ratio is a rough estimate of the amount of rounding error you can expect in this quantity.

To run this test, add the code block

%# Finite difference test
Ns=4; %# Number of states

randn('seed',0.2004);
W=(randn(prod(S),Ns)+i*randn(prod(S),Ns));

more off; %# View output as it is computed
fdtest(W);

to the bottom of your file ``sch.m''. Then run ``setup.m'' and ``sch.m'' and verify that you observe the ratio approach unity, gaining one order of magnitude per iteration before rounding error dominates the ratio.

Notes: The above code block sets the number of states for our problem to Ns=4, ``seeds'' the random number generator so that we all will get the same results, provides an initial random complex $W$ for ``fdtest.m'', does ``more off'' so that you can immediately see outputs as they are computed, and then calls fdtest.

Hints if your code fails the above test:

If your code fails the test, then it is helpful to debug the kinetic and potential energy parts separately. To test the kinetic energy part alone just delete (or comment out) the Vdual parts in both getE() and H(), and rerun. Then, repeat for the potential energy by putting back the Vdual parts and commenting out the L parts in both getE() and H(). If one of these works, but not the other, then you have isolated the problem.

If neither the potential nor the kinetic parts work, then the problem is likely in the algebra in your getgrad(). One way to test for this and to be able to debug H() and getE() independently of this extra algebra is to start with an initially orthonormal $W$. You can create such by including the statement ``W=W*inv(sqrtm(W'*O(W))));'' immediately before the call to fdtest. Be certain, however, to remove this statement and test again once you have identified your bug(s). For the rest of this problem set to function, it is critical that getgrad() and getE() work with non-orthonormal functions.