[ruby-list:39324] Re: 固有値、固有ベクトルの計算

皆さん今晩は、田鎖です。

DL を使った dgeev です。DL は、この頃使い出したので、本当にこれで良い
のかと言う部分が、多々ありますが、手軽で便利なのが気に入っています。自
分の環境は、linux(i386)ですが、他の環境(特にWindows)でもこのようなプロ
グラムが動くのかどうか非常に興味のあるところです。

						'04.03.05 tag


----------------------------------------------------------- ここから

#!/usr/local/bin/ruby 

#
# LAPACK support using DL module
#

require "dl/import"
require "dl/struct"

module LAPACK
  extend DL::Importable
  dlload "liblapack.so"
  dlload "libblas.so"

  extern "void dgeev_(char*,char*,int*,double*,int*,double*,double*,double*,int*,double*,int*,double*,int*,int *)"

  module_function

  def dgeev(a)
    jobvl = [?N].pack("C").to_ptr
    jobvr = [?V].pack("C").to_ptr
    _n = a.size
    n = [_n].pack("L").to_ptr
    aa = Marshal.load(Marshal.dump(a))
    lda = [_n].pack("L").to_ptr
    wr = ([0] * _n).pack("D*").to_ptr
    wi = ([0] * _n).pack("D*").to_ptr
    vl = ([0] * _n).pack("D*").to_ptr
    ldvl = [1].pack("L").to_ptr
    vr = ([0] * _n * _n).pack("D*").to_ptr
    ldvr = [_n].pack("L").to_ptr
    work = ([0] * 8 * _n).pack("D*").to_ptr
    lwork = [8 * _n].pack("L").to_ptr
    info = [0].pack("L").to_ptr
    dgeev_(jobvl, jobvr, n, aa.flatten.to_ptr, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
    err = info.to_a("L")[0]
    raise err if err != 0
    vr.to_a("D").each_with_index {|d, i| aa[i/_n][i%_n] = d}
    [aa, wr.to_a("D"), wi.to_a("D")]
  end

end

if __FILE__ == $0
  require "pp"

  a = [
    [1.0, 1.0, 2.0],
    [0.0, 2.0, 2.0],
    [-1.0, 1.0, 3.0],
  ]

  vr, wr, wi = LAPACK.dgeev(a)
  pp a, vr, wr, wi
end

----------------------------------------------------------- ここまで