function [D, Dt, DtD] = DiffOper_d(row, col)
N = row * col;
Dv = spdiags([reshape([-ones(row-1,col); zeros(1,col)],N,1) ...
            reshape([zeros(1,col); ones(row-1,col)],N,1)], [0 1], N, N);
Dh = spdiags([reshape([-ones(row,col-1) zeros(row,1)],N,1) ...
            reshape([zeros(row,1) ones(row,col-1)],N,1)], [0 row], N, N);
D = sparse([Dh;Dv]);   
Dt = D';
DtD = Dt*D;
end