Skip to content

Instantly share code, notes, and snippets.

@ahmadia
Created October 15, 2012 16:17
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save ahmadia/3893361 to your computer and use it in GitHub Desktop.
basic cythonized version of Jed's code
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
@ahmadia
Copy link
Author

ahmadia commented Oct 15, 2012

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment