C类型

Ctypes是一个非常有趣的python包，它允许您将共享对象库导入python并直接调用它们。我应该说，即使这叫做ctypes，它也可以用来调用用fortran编写的库中的函数。唯一复杂的是你需要知道一个fortran函数对C来说是什么样子。这很简单，所有的东西都是一个指针，所以如果你的fortran函数被称为foo（a，N），其中a是一个数组，N是它的长度，那么从C调用它需要一个指向双精度数组的指针和一个指向int的指针。另一个需要注意的是在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应该打电话总和动态库那你就可以了

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) 返回一个ctypes对象，该对象是指向a的基础数组的空指针。请注意，即使sum接受double * ，只要我们有一个指向正确数据的指针，它的类型是什么都无关紧要，因为它将被自动转换。

请注意，实际上还有其他方法传递所需的double数组。例如

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

此示例仅使用ctypes。Ctypes有C数据类型的包装器，例如

a=c_double(10.4)

将创建一个可以传递给C函数的ctypes double对象。请注意，有一个byref函数可以让您通过引用传递参数。这将在下一个示例中使用。c双倍 * 10，是一个python对象，它表示10个double和

a=(c_double*10)()

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

下面是一个使用ctypes直接调用lapack的示例。请注意，只有在您的系统上有一个lapack共享对象库时，这才有效。同样在linux上，该文件将是利伯拉帕克您可能会使用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.restype=c加倍。默认情况下，ctypes假定返回值是int。如果不是，则需要通过将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级别完成整个解决方案，而且我们的速度也会相同。

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