» Using LAPACK in C Cheng Meng’s Small Talk

Using LAPACK in C

二月 22, 2008 at 1:00 下午 · Filed under Default, Technology

简单的说就是如何在C程序中调用Fortran子程序。我不想再学一门语言，所以直接调用库函数是最好的办法。google了一下，在VS2005里调试成功。简单来说你只需要：

编译好的LAPACK库

相应的头文件

知道Fortran数组内存布局的特点

知道如何在C程序中调用DLL

就没问题。dll和头文件可以在这里下载。建议动手前查看一下头文件中的f2c.h，里面定义了一些常用的数据类型，混个眼熟。

怎么在C程序中链接外部DLL库，这个和编译器有关，就不在这儿罗嗦了。

主要的问题是怎样在C的数据类型和Fortran中的类型间进行转换。基本的类型，比如REAL, DOUBLEREAL之类的定义，都在f2c.h里写好了。最要紧的一点：Fortran中二维数组是按列存储的，而C是按行存储。当然，我个人的意见是干脆不要用C提供的多维数组，直接用一维数组，元素按Fortran惯例存放就行了，顺便定义一套宏什么的，用起来省事，看起来顺眼。此外，Fortran中函数传参是call by reference，和C的call by value截然相反，所以所有参数都必须是指针。记住这两条就差不多了。

我要用LAPACK对角化Hamiltonian，感觉应该用zheev, compute all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix

附精美范例一则，在VS 2005下编译通过：

#include <stdio.h>
#include "f2c.h"
#include "zgeev.h"

#define SIZE 3

int main(int argc, char *argv[]) {
  doublecomplex A[SIZE][SIZE], b[SIZE];
  doublecomplex DUMMY[1][1],WORK[2*SIZE];
  doublecomplex AT[SIZE*SIZE];
  int i, j, ok, c1, c2, c3;
  char c4;

  /* initializing A */

  for (i = 0; i < SIZE; i++) {
    for (j = 0; j < SIZE; j++) {
      AT[j + SIZE*i].r = A[j][i].r;
      AT[j + SIZE*i].i = A[j][i].i;
    }
  }

  c1 = SIZE;
  c2 = 2*SIZE;
  c3 = 1;
  c4 = ‘N’;

  zgeev_(&c4, &c4, &c1, AT, &c1, b,
        DUMMY, &c3, DUMMY, &c3, WORK, &c2, WORK, &ok);

  if (ok == 0)
    for (i = 0; i < SIZE; i++)
      printf("%f + %fi\n", b[i].r, b[i].i);
  else
    printf("An error occured");

  return 0;
}