basic cythonized version of Jed's code

jed.pyx

cimport cython

cimport numpy as np

import numpy as np

DTYPE = np.float
ctypedef np.float DTYPE_t


from libc.math cimport exp, sqrt, pow, cos

cdef float pi = 3.14159265


@cython.cdivision(True)
@cython.boundscheck(False)
def logg(int m):
    D,xb = Cheb(m+1); 
    D = 2*D; xb = 0.5 + 0.5*xb; 
    # D2 = -(D*D)
    # D2 = D2(2:end-1,2:end-1); % Negative 1D Laplacian
    D2 = -(np.dot(D,D))[1:-1][:,1:-1]; 
    
    # xd = xb(2:end-1);
    xd = xb[1:-1]

    # sinx = sin(pi*xd);
    sinx = np.sin(pi*xd)

    # xexact_1 = (kron(sinx', sinx))
    # xexact = xexact_1(:)
    xexact = np.kron(sinx,sinx)

    #b = 2*pi^2*xexact;
    b = 2*np.pi**2*xexact

    # Compute action of inverse of L using sum factorization
    # See Lynch, Rice, and Thomas (1964)
    
    # [R, Lambda] = eig(D2);
    Lambda,R = np.linalg.eig(D2);
    
    #Lambda = diag(Lambda);
    Lambda = Lambda.reshape((m,1))

    #e = ones(m,1);
    e = np.ones((m,1));

    #Diag = kron(Lambda,e) + kron(e,Lambda);
    Diag = np.kron(Lambda,e) + np.kron(e,Lambda);

    #iDiag = reshape(1./Diag,m,m);
    iDiag = np.reshape(1./Diag,(m,m))
    # Rinv = inv(R);
    Rinv = np.linalg.inv(R)

    # y = iDiag .* (Rinv*reshape(b,m,m)*Rinv');
    y = iDiag * (np.dot(Rinv,np.dot(np.reshape(b,(m,m)),Rinv.T)))

    # x = (R*y*R')(:);                          % x = Linv * b
    x = np.dot(R,np.dot(y,R.T))

    error = np.linalg.norm(x.flatten() - xexact.flatten());

    return error

@cython.boundscheck(False)
def Cheb(int N):
  """Constructs Chebyshev differentiation matrix and nodes"""
  if (N==0):
   D=0
   x=1; 
   return D,x 

#  x = cos(pi*(0:N)/N)'; 
  x = np.cos(pi*np.arange(N+1)/N)[:, np.newaxis]
  #c  = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  c_1 = np.ones((1,N+1)); c_1[0,0] = 2; c_1[0,N] = 2
  c = (c_1*(-1)**np.arange(N+1)).T

  #X = repmat(x,1,N+1);
  X = np.tile(x,(1,N+1))

  #dX = X-X';                  
  dX = X-X.T

  #D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
  D  = (c*(1./c).T)/(dX+(np.eye(N+1)))
  # D  = D - diag(sum(D'));                 % diagonal entries
  D  = D - np.diag(np.sum(D,1)); 

  return D,x

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [Extension("jed", ["jed.pyx"])]

import numpy as np

setup(
  name = 'Jed',
  cmdclass = {'build_ext': build_ext},
  ext_modules = ext_modules,
  include_dirs=[np.get_include()]
)

setup.py file if needed