How about calling your FORTRAN routines directly from ruby?
I recently started on a Fortran-oriented linear algebra package in
which I do this. However I put off working on it until I was
convinced a show-stopper crash in ruby/dl could be fixed. As it so
happens, it was fixed today [ruby-core:2095].
I have yet to form a project on rubyforge, but probably all you need
is the following example:
require ‘dl/import’
require ‘dl/struct’
require ‘evil’ # http://evil.rubyforge.org
module Lapack
extend DL::Importable
dlload “liblapack.so”
dlload “libblas.so”
dlload “libf2c.so”
typealias “doublereal”, “double”
typealias “integer”, “int”
typealias “ftnlen”, “short”
def self.prototype s
func, args = s.strip.gsub(%r!\s+!, " “).split(%r![()]!)
args = args.gsub(%r!(\w+[\s])\w*!){$1}.gsub(%r!\s+!,”")
extern “#{func}(#{args})”
end
prototype %{
void s_copy(char *a, char *b, ftnlen la, ftnlen lb);
}
prototype %{
int dgemm_(char *transa,
char *transb,
integer *m,
integer *n,
integer *k,
doublereal *alpha,
doublereal *a,
integer *lda,
doublereal *b,
integer *ldb,
doublereal *beta,
doublereal *c,
integer *ldc)
}
end
class DMatrix
ELEMTYPE = “d”
ELEMTYPE_GLOB = “d*”
attr_reader :vsize, :hsize
def initialize(columns, copy = true)
@vsize = columns[0].size
@hsize = columns.size
if copy
@columns = columns.map { |col| col.map { |e| e.to_f } }
else
@columns = columns
end
end
def self.columns(columns, copy = true)
self.new(columns, copy)
end
def (i, j) ; @columns[j][i] ; end
def =(i, j, e) ; @columns[j][i] = e.to_f ; end
def to_s
res = “”
(0…vsize).each { |i|
(0…hsize).each { |j|
res << sprintf(“% 10-.6f”, self[i,j])
}
res << “\n”
}
res
end
def data
@columns.flatten.pack(ELEMTYPE_GLOB)
end
def self.data(vsize, hsize, data)
flat = data.unpack(ELEMTYPE_GLOB)
cols = (0…hsize).map { flat.slice!(0, vsize) }
self.new(cols, false)
end
def *(other)
m = self.vsize
n = other.hsize
k = other.vsize
raise "matrix size mismatch" unless self.hsize == k
a = self
b = other
c_data = [0.0].pack(ELEMTYPE)*(m*n)
begin
# prevent dangling pointers
GC.disable
Lapack.dgemm_("N", # transa
"N", # transb
[m].to_ptr, # m
[n].to_ptr, # n
[k].to_ptr, # k
[1.0].to_ptr, # alpha
a.data.internal.ptr, # a
[m].to_ptr, # lda
b.data.internal.ptr, # b
[k].to_ptr, # ldb
[0.0].to_ptr, # beta
c_data.internal.ptr, # c
[n].to_ptr) # ldc
ensure
GC.enable
end
DMatrix.data(m, n, c_data)
end
end
a = DMatrix.columns [ [1,2], [3,4], [5,6], [7,8] ]
b = DMatrix.columns [ [9,10,11,12], [13,14,15,16] ]
[a, b, a*b].each { |m|
puts “-”*50
puts m
}
···
— Qubert qubert@orbit.gotdns.com wrote:
OK, I have an algorithm that I created to format a series of numbers
for a column centric FORTRAN routine to use. I don’t like what I have
created, although it works, so I would like to ask if anyone here has
a more cleaver way.
Do you Yahoo!?
Yahoo! Domains Claim yours for only $14.70/year
http://smallbusiness.promotions.yahoo.com/offer
