読者です 読者をやめる 読者になる 読者になる

Lapack C interface, lapacke

http://www.netlib.org/lapack/lapacke.html

lapackeはlapackのパッケージに入っている。たとえば、
lapack-3.4.2/lapacke
に入って、make lapackeすればコンパイルできる。
make lapackeする前に、lapack-3.4.2/make.incを書く必要がある。
Ubuntu <= 12.04の場合、lapack-devをインストールしておいて、

CFLAGS = -O3 -fPIC
# LAPACKLIB    = liblapack.a
# BLASLIB      = ../../librefblas.a

liblapacke.aができるので、lapack-3.4.2/lapacke/includeとともに、必要なところにコピーすれば使える。使うときは、-llapack -lblasとしてリンクを忘れないように。

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

int main (int argc, const char * argv[])
{
  lapack_int info, n, lda;
  double w[6];
  int i,j;
  lapack_complex_double a[36];
  double a_d[] = {
    0.459759,  0.000000,  -0.003394,  0.000000,  0.008884,  -0.000000,  -0.411255,  -0.005905,  -0.006039,  -0.136799,  0.020047,  0.041309,  -0.003394,  -0.000000,  0.459759,  0.000000,  -0.008884,  -0.000000,  -0.006039,  -0.136799,  -0.411255,  -0.005905,  -0.020047,  -0.041309,  0.008884,  0.000000,  -0.008884,  0.000000,  0.491964,  0.000000,  0.020047,  0.041309,  -0.020047,  -0.041309,  -0.411255,  -0.005905,  -0.411255,  0.005905,  -0.006039,  0.136799,  0.020047,  -0.041309,  0.459759,  0.000000,  -0.003394,  -0.000000,  0.008884,  0.000000,  -0.006039,  0.136799,  -0.411255,  0.005905,  -0.020047,  0.041309,  -0.003394,  0.000000,  0.459759,  0.000000,  -0.008884,  0.000000,  0.020047,  -0.041309,  -0.020047,  0.041309,  -0.411255,  0.005905,  0.008884,  0.000000,  -0.008884,  -0.000000,  0.491964,  0.000000};

  n = 6;
  lda = 6;

  for (i = 0; i < 36; i++) {
    a[i] = lapack_make_complex_double(a_d[i * 2], a_d[i * 2 + 1]);
  }
  info = LAPACKE_zheev(LAPACK_ROW_MAJOR,'V', 'U', n, a, lda, w);

    
  for(i = 0; i < n; i++) {
    if (w[i] < 0) {
      printf("%f\n", -sqrt(-w[i]) * 15.633302);
    } else {
      printf("%f\n", sqrt(w[i]) * 15.633302);
    }
    for (j = 0; j < n; j++) {
      printf("%f %f\n",
	     lapack_complex_double_real(a[j * n + i]),
	     lapack_complex_double_imag(a[j * n + i]));
    }
  } 
  printf("\n");
  return(info);
}