CTYPE

Ctype是一个非常有趣的Python包，它允许您将共享对象库导入到Python中并直接调用它们。我应该说，即使这被称为ctype，它也可以用来调用用Fortran编写的库中的函数。唯一的复杂之处是你需要知道Fortran函数对于C是什么样子的。这很简单，一切都是一个指针，所以如果你的Fortran函数被调用为foo(A，N)，其中A是一个数组，N是它的长度，那么从C调用它，它接受一个指向双精度数组的指针和一个指向整型数组的指针。另一件需要注意的事情是，在C语言中，Fortran函数通常带有下划线。也就是说，Fortran函数foo将显示为C中的foo(通常是这种情况，但取决于编译器)。话虽如此，下面的例子是用C编写的。

例如，假设您编写了以下简单的C程序

#include <stdio.h>

int sum(double *x,int n)
{
  int i;
  double counter;
  counter = 0;
  for(i=0;i<n;i++)
    {
      counter=counter+x[i];

    }
  return counter;
}

您想要从Python调用的。首先通过执行以下操作创建共享对象库(在命令行中)

$ gcc -c sum.c
$ gcc -shared -o sum.so sum.o

请注意，在OSX上，-Shared应替换为-Dynamiclib，而Sum.so应称为Sum.dylib，然后您可以

from ctypes import *
my_sum=CDLL('sum.so')
a=numpy.array(range(10),dtype=float)
my_sum.sum(a.ctypes.data_as(c_void_p),int(10))

请注意，这里 a.ctypes.data_as(c_void_p) 返回一个ctype对象，它是指向基础数组的空指针。请注意，即使sum接受一个双精度*，只要我们有指向正确数据的指针，它的类型并不重要，因为它将被自动强制转换。

请注意，实际上还有其他方法可以传入所需的Double数组。例如

a=(c_double*10)()
for i in range(10):
   a[i]=i
my_sum.sum(a,int(10))

此示例仅使用ctype。Ctype具有用于C数据类型的包装器，例如

a=c_double(10.4)

将创建一个ctype双重对象，该对象可以传递给C函数。请注意，有一个byref函数允许您通过引用传递参数。这将在下一个示例中使用。C Double*10，是一个表示10个Double和

a=(c_double*10)()

将设置为等于10个双精度数组。我发现当数据是数学数据时，这种方法通常不如使用NumPy数组有用，因为NumPy数组更好地集成到了Python和Sage中。

下面是一个使用ctype直接调用Lapack的示例。请注意，只有当您的系统上有Lapack共享对象库时，这才能起作用。同样在Linux上，文件将是liblapack.so，您可能会使用dgesv(OSX使用CLAPACK，因此没有下划线)。

from ctypes import *
def ctypes_solve(m,b,n):
    a=CDLL('/usr/lib/liblapack.dylib')
    import numpy
    p=(c_int*n)()
    size=c_int(n)
    ones=c_int(1)
    ok=c_int(0)
    a.dgesv(byref(size),byref(ones),m.ctypes.data_as(c_void_p),
            byref(size),p,b.ctypes.data_as(c_void_p),byref(size),byref(ok))

为完整起见，让我们考虑一种使用C型求解拉普拉斯方程的方法。假设您已经用C语言编写了一个简单的求解器，并且您希望从Python中调用它，这样您就可以轻松地测试不同的边界条件。您的C程序可能如下所示。

#include <math.h>
#include <stdio.h>

double timestep(double *u,int nx,int ny,double dx,double dy)
{
  double tmp, err, diff,dx2,dy2,dnr_inv;
  dx2=dx*dx;
  dy2=dy*dy;
  dnr_inv=0.5/(dx2+dy2);
  err = 0.0;
  int i,j;

for (i=1; i<nx-1; ++i) {
  for (j=1; j<ny-1; ++j) {
    tmp = u[i*nx+j];
    u[i*nx+j] = ((u[(i-1)*nx+j] + u[(i+1)*nx+j])*dy2 +
          (u[i*nx+j-1] + u[i*nx+j+1])*dx2)*dnr_inv;
    diff = u[i*nx+j] - tmp;
    err += diff*diff;
  }
}

 return sqrt(err);
}

double solve_in_C(double *u,int nx,int ny,double dx,double dy)
{
  double err;
  int iter;
  iter = 0;
  err = 1;
    while(iter <10000 && err > 1e-6)
      {
    err=timestep(u,nx,ny,dx,dy);
    iter++;
      }

  return err;
}

我们可以通过在命令行运行来编译它

$ gcc -c laplace.c
$ gcc -shared -o laplace.so laplace.o

现在在SAGE(笔记本或命令行)中执行

from ctypes import *
laplace=CDLL('/home/jkantor/laplace.so')
laplace.timestep.restype=c_double
laplace.solve_in_C.restype=c_double
import numpy
u=numpy.zeros((51,51),dtype=float)
pi_c=float(pi)
x=numpy.arange(0,pi_c+pi_c/50,pi_c/50,dtype=float)
u[0,:]=numpy.sin(x)
u[50,:]=numpy.sin(x)

def solve(u):
  iter =0
  err = 2
  n=c_int(int(51))
  pi_c=float(pi/50)
  dx=c_double(pi_c)
  while(iter <5000 and err>1e-6):
     err=laplace.timestep(u.ctypes.data_as(c_void_p),n,n,dx,dx)
     iter+=1
     if(iter %50==0):
        print((err,iter))
  return (u,err,iter)

请注意laplace.timestep.rest ype=c Double行。默认情况下，ctype假定返回值为整型。如果它们不是，您需要通过将RESTYPE设置为正确的返回类型来告诉它。如果执行上述代码，则solve(U)将求解系统。它可以与Fortran解决方案相媲美，只需大约0.2秒。或者，您也可以这样做

n=c_int(int(51))
dx=c_double(float(pi/50))
laplace.solve_in_C(n.ctypes.data_as(c_void_p),n,n,dx,dx)

完全用C语言计算，这是非常快的。诚然，我们可以让我们的Fortran例程在Fortran级别上完成整个解决方案，我们将拥有相同的速度。

正如我在前面所说的，您可以使用ctype轻松地调用用Fortran编写的共享对象库。关键是它必须是一个共享对象库，并且所有Fortran参数都是通过引用传递的，即作为指针或使用byref。而且，即使我们使用非常简单的数据类型，也可以处理更复杂的C结构。有关ctype的详细信息，请参阅http://python.net/crew/theller/ctypes/