Mathematica test code

(* All of this is already defined before hand *)
L = 1200;
t = -1.0;
w = -1.0;
e = Table[-2 t Cos[(n \[Pi])/(L + 1.)], {n, 1, L}];
f = Table[Sqrt[2/(L + 1.)] Sin[(n \[Pi])/(L + 1)], {n, 1, L}];

(* Naive do loop *)
AbsoluteTiming[h = Table[0, {2 L}, {2 L}];]
AbsoluteTiming[Do[h[[i, i]] = e[[i]], {i, 1, L}];]
AbsoluteTiming[Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];]
AbsoluteTiming[
 Do[h[[i, j]] = w f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];]
AbsoluteTiming[Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];]

(* A more idiomatic way using built in functions *)
AbsoluteTiming[h1 = w ArrayPad[Outer[Times, f, f], {{0, L}, {L, 0}}];]
AbsoluteTiming[h1 += Transpose[h1];]
AbsoluteTiming[h1 += DiagonalMatrix[Join[e, e]];]

(* Compiled version of the do loop *)
cf = Compile[{{L, _Integer}, {e, _Real, 1}, {f, _Real, 1}},
   Module[{h},
    h = Table[0.0, {2 L}, {2 L}];
    Do[h[[i, i]] = e[[i]], {i, 1, L}];
    Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];
    Do[h[[i, j]] = -f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];
    Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];
    h], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

AbsoluteTiming[h2 = cf[L, e, f];]

(* Sanity checks. This may take a long time, so you should probably \
reduce L to something more sane *)
MatrixPlot[Chop[h - h1]]
MatrixPlot[Chop[h - h2]]
MatrixPlot[h]
MatrixPlot[h1]
MatrixPlot[h2]

k = 3 L/2;

(* Naive version *)
AbsoluteTiming[
 p = Table[If[i == j && j <= k, 1, 0], {i, 2 L}, {j, 2 L}];]

(* Create identity matrix, then pad it *)
AbsoluteTiming[
 p2 = ArrayPad[IdentityMatrix[k], {{0, 2 L - k}, {0, 2 L - k}}];]

(* Use a sparse array representation *)
AbsoluteTiming[p3 = SparseArray[{i_, i_} /; i <= k -> 1, {2 L, 2 L}];]

(* Sanity checks *)
MatrixPlot[Chop[p - p2]]
MatrixPlot[Chop[p - p3]]
MatrixPlot[p]
MatrixPlot[p2]
MatrixPlot[p3]