From: Jeff Mitchell <quixoticsycophant@...>
Date: 2004-05-20T12:35:56+09:00
Subject: Re: Fortran Format?

--- 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.

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
}



	
		
__________________________________
Do you Yahoo!?
Yahoo! Domains � Claim yours for only $14.70/year
http://smallbusiness.promotions.yahoo.com/offer