For debugging purposes, you should now run your software for DFT calculations in a simple-harmonic oscillator potential (``dft.m'') with the parameters S=[32; 32; 32] and R=diag([8 8 8]). For convenience, you may wish to comment out the parts of ``dft.m'' generating the plots, as we will not be looking at them. (Do keep the call to getPsi(), however.) Take note of your final converged values, as we will be checking that your reduced basis calculation gives essentially the same result. (I get E=46.594, epsilon=[6.320;7.742;7.742;7.742].)
Of the operators, it turns out that Linv(), cJ() and cJdag() are never
called to act on the wave functions, and so need no modification.
(Ultimately, the fact that your software runs unmodified without octave
generating any ``nonconformant'' size errors confirms this.) Also, as
shown in lecture, the overlap operator is still just multiplication by
the same constant,
, and so requires no
modification. This, then, leaves the following operators to be
modified.