getE()

function E=getE(W)

Input:

• W: Expansion coefficients for $N_s$ unconstrained wave functions, stored as an $\prod S_k\,\times\,N_s$ matrix

Global variables:

• gbl_Vdual: Dual potential coefficients stored as a $S_k\,\times\,1$ column vector.

Output:

• E: Energies summed over $N_s$ states

Using the expressions in Eq. (8) and your diagouter() operator above, produce a function of the above prototype in the file ``getE.m''.

As a quick test, although you will generally find machine-precision sized imaginary parts due to rounding, your output should always be real. You may verify this behavior with

setup; sch; %# Make sure your global variables are all set
W=randn(prod(S),4)+i*randn(prod(S),4); %# Put 4 random wavefunctions in W
getE(W)

After verifying that your output is indeed real to machine precision for a few different random input $W$'s, it is best to modify your code to take the real part of E using real() before returning, so as to avoid dealing with complex numbers in inappropriate places later.