Debugging

The ultimate test of a good preconditioner is whether it actually improves the rate of convergence. You may note, however, that if you replace your current call to lm() in ``dft.m'' with pclm(), the energy may begin to increase momentarily. This is not uncommon when starting a minimization far from the minimum. To make a more appropriate comparison replace your current call to lm() with

%# Converge
W=sd(W,20); %# 20 iterations of simple sd() to get nearer to the minimum
W=W*inv(sqrtm(W'*O(W))); %# Restart as orthonormal functions

[Wlm, Elm]=lm(W,50); %# 50 iterations of lm() from W, while recording results
[Wpclm, Epclm]=pclm(W,50); %# 50 iterations of pclm from same W

%# Plot results on log scale with nice labels
grid on;
xlabel("Iteration (#) ->"); ylabel("Error (H)");
semilogy([1:length(Elm)],Elm-43.3371147782040,'r-;lm;', ...
         [1:length(Epclm)],Epclm-43.3371147782040,'b-;pclm;');
pause; %# Wait for keyboard input before continuing

This code block first ``relaxes'' W with 20 iterations and orthonormalizes the result to give a good starting point relatively near the minimum. The code block then runs both lm() and pclm() from that this starting point, each for 50 iterations, and plots on a log scale the difference of each algorithm from the fully relaxed result, with line minimization in red, and preconditioned line minimization in blue. If the preconditioner is working, you should notice both some initial difference coming from the first iteration where the preconditioner eliminates nearly all of the high frequency errors and a difference in slope (exponential convergence rate) at the higher iteration numbers.
Note: After the plot, the program will halt at the ``pause;'' statement, giving the option of exiting by hitting Ctrl-C before generating the plots.