function M = mlaplacian(g,m)
% size of g = 3^(m+1)

M = zeros(3^(m+1)) ;

for k = 1:3^(m+1)
  v = address(k,m) ;
  vnbhd = vneighbors(v,m) ;
  gnbhd = gneighbors(v,m) ;
  
  for j = 1:4
    l = indexsg(gnbhd(j,:),m) ;

    n1 = indexsg(vnbhd(j,:),m) ;
    n2 = alternateindex(n1,m) ;
    
    M(k,k) = M(k,k) + g(l) ;
    M(k,n1) = M(k,n1) - g(l) ;
    if(n2 ~= n1)
      M(k,n2) = M(k,n2) - g(l) ;
    end 
  end
end
