clear Wk*

N = 500;
Nx = 10000;

[x w] = hermpts(Nx);
tic;
B = hermiteS(0:N-1,x);
fprintf('Basis auswerten:         %s\n' ,szeit(toc));

Vks = diag(sparse(w' .* (V(x) + VH(x))));
tic;
Wk = B * Vks * B.';
fprintf('Referenz: %s\n', szeit(toc));

dx = gpuArray(x);
dw = gpuArray(w);
dVk1 = (dw' .* (V(dx) + VH(dx))) * gpuArray(ones(1,N));

tic;
dB = cuhermite(0:N-1,dx);
fprintf('Basis auf GPU auswerten: %s\n' ,szeit(toc));

tic;
dWk1 = dB * (dVk1 .* dB.');
fprintf('Rechnung auf GPU:    %s      ', szeit(toc));
serr(Wk, gather(dWk1));

%dVks = diag(sparse(dw' .* (V(dx) + VH(dx))));

dVk2 = dw' .* (V(dx) + VH(dx));
kern = parallel.gpu.CUDAKernel('diagMult.ptx','diagMult.cu');
kern.ThreadBlockSize = 128;

tic;
dWk2 = dB * feval(kern, dB.', dVk2, dB.', Nx, N);
fprintf('Rechnung mit Kernel: %s      ', szeit(toc));
serr(Wk, gather(dWk2));

whos dVk1 dVk2;