RutherfordBoeing.jl

module RutherfordBoeing

export RutherfordBoeingData


type RBMeta
  # Metadata attached to Rutherford-Boeing data.
  title :: String
  key   :: String

  mxtype :: String
  nrow :: Int
  ncol :: Int
  nnzero :: Int
  neltvl :: Int

  hermitian :: Bool
  assembled :: Bool
  pattern_only :: Bool

  ptrfmt :: String
  indfmt :: String
  valfmt :: String

  dattyp :: String
  positn :: String
  orgniz :: Char
  caseid :: String
  numerf :: Char
  auxfm1 :: String
  auxfm2 :: String
  auxfm3 :: String
end


RBDataType = Union(Array{Float64,2}, Array{Complex64,2}, SparseMatrixCSC)


type RutherfordBoeingData

  meta :: RBMeta
  data :: RBDataType

  function RutherfordBoeingData(file_name :: String)
    data_types = ['r' => Float64, 'c' => Complex64, 'i' => Int64]
    rb = open(file_name)

    # Read header.
    seekstart(rb)
    line  = readline(rb)
    title = strip(line[1:72])      # A72
    key   = strip(line[72:end-1])  # A8
    buffer1 = readline(rb)         # A80
    buffer2 = readline(rb)         # A80

    println(buffer1)
    println(length(buffer1))
    println(buffer2)
    println(length(buffer2))

    if buffer2[3] in ['a', 'e']
      println("matrix")
      # Read a matrix.
      mxtype = buffer2[1:3]
      @printf("mxtype = %s\n", mxtype)

      nrow, ncol, nnzero, neltvl = map(int, split(chomp(buffer2[4:end])))
      @printf("nrow=%d ncol=%d nnzero=%d neltvl=%d\n", nrow, ncol, nnzero, neltvl)

      pattern_only = (mxtype[1] in ['p', 'q'])
      if pattern_only
        ptrfmt, indfmt = split(strip(readline(rb)))
        @printf("ptrfmt=%s indfmt=%s\n", ptrfmt, indfmt)
      else
        data_type = data_types[mxtype[1]]
        ptrfmt, indfmt, valfmt = split(strip(readline(rb)))
        @printf("ptrfmt=%s indfmt=%s valfmt=%s\n", ptrfmt, indfmt, valfmt)
      end

      # Read integer data.
      np1 = (mxtype[2:3] == "re") ? (2*ncol + 1) : (ncol + 1)
      @printf("np1=%d\n", np1)

      ip  = read_array(rb, np1, ptrfmt)
      display(ip'); println()

      ind = read_array(rb, nnzero, indfmt)
      display(ind'); println()

      nreal = neltvl > 0 ? neltvl : nnzero
      @printf("nreal=%d\n", nreal)

      if pattern_only
        vals = ones(Float64, nreal)  # To be able to build a SparseMatrixCSC
      else
        vals = read_array(rb, nreal, valfmt)
      end
      display(vals'); println()
      data = SparseMatrixCSC(nrow, ncol, ip, ind, vals)

      meta = RBMeta(title, key, mxtype, nrow, ncol, nnzero, neltvl,
                    mxtype[2] == 's', mxtype[3] == 'a', pattern_only,
                    ptrfmt, indfmt, pattern_only ? "" : valfmt,
                    "", "", '\0', "", '\0', "", "", "")

    else

      println("data")
      # Read supplementary data.
      dattyp = buffer1[1:3]
      positn = buffer1[4]
      orgniz = buffer1[5]
      caseid = buffer1[7:14]
      numerf = buffer1[16]
      nrow, ncol, nnzero = map(int, buffer1[17:end])
      auxfm1, auxfm2, auxfm3 = split(strip(buffer2))

      if (dattyp == "rhs" & orgniz == 's') |
         (dattyp in ["ipt", "icv", "ord"])
        if dattyp == "ord"
          nnzero = nrow * ncol
          auxfm = auxfm1
        else
          ip = read_array(rb, ncol+1, auxfm1)
          auxfm = auxfm2
        end

        ind = read_array(rb, nnzero, auxfm)
      end

      data_type = data_types[numerf]
      if dattyp != "rhs"
        nnzero = nrow * ncol
      end
      if dattyp == "rhs" & orgniz == 's'
        auxfm = auxfm3
      else
        auxfmt = auxfm1
      end

      val = read_array(rb, nnzero, auxfm)

      if (dattyp == "rhs" & orgniz == 's') |
         (dattyp in ["ipt", "icv", "ord"])
        data = SparseMatrixCSC(nrow, ncol, ip, ind, vals)

      else
        if dattyp != "rhs"
          data = reshape(vals, (nrow, ncol))
        else
          data = vals
        end
      end

      meta = RBMeta(title, key, mxtype, nrow, ncol, nnzero, 0,
                    false, orgniz != 'e', false,
                    "", "", "",
                    dattyp, positn, orgniz, caseid, numerf,
                    auxfm1, auxfm2, auxfm3)
    end

    close(rb)

    new(meta, data)
  end
end


# Helper functions.

function decode_int_fmt(fmt :: String)
  if fmt[1] == '('
    fmt = uppercase(fmt[2:end-1])
  end
  return map(int, split(fmt, 'I'))
end


function decode_real_fmt(fmt :: String)
  fmt = join(split(fmt))  # Remove all white spaces.
  if fmt[1] == '('
    fmt = uppercase(fmt[2:end-1])
  end
  scale = 0
  if (',' in fmt)  # Process scale factor, e.g., 1P,5D16.9
    scale, fmt = split(fmt, ',')
    scale, _ = split(scale, 'P')
  elseif ('P' in fmt)
    scale, fmt = split(fmt, 'P')
  end
  scale = int(scale)

  fmt1 = split(fmt, '.')[1]
  if search(fmt1, 'E') > 0
    (npl, len) = map(int, split(fmt1, 'E'))
  elseif search(fmt1, 'D') > 0
    (npl, len) = map(int, split(fmt1, 'D'))
  elseif search(fmt1, 'F') > 0
    (npl, len) = map(int, split(fmt1, 'F'))
  else
    error("Malformed real format")
  end

  return (npl, len, scale)
end


function standardize_real(number_as_str :: String)
  s = join(split(number_as_str))  # for numbers in the form "0.24555165E 00".
  # change ".16000000+006" to ".16000000e+006". The first char could be +/-.
  if (search(s, 'E') + search(s, 'D') + search(s, 'e') + search(s, 'd')) == 0
    if search(s[2:end], '+') > 0
      s = s[1:1] * join(split(s[2:end], '+'), "e+")
    elseif search(s[2:end], '-') > 0
      s = s[1:1] * join(split(s[2:end], '-'), "e-")
    end
  end
  return s
end


function read_array(io :: IO, n :: Int, fmt :: String; is_complex :: Bool = false)
  if 'I' in fmt
    scale = 0
    (npl, len) = decode_int_fmt(fmt)
    conv = int
    typ = Int64
  else
    (npl, len, scale) = decode_real_fmt(fmt)
    conv = float
    typ = Float64
    if is_complex
      n *= 2
    end
  end

  x = zeros(typ, n)
  for j = 1 : div(n, npl)
      line = join(split(uppercase(readline(io)), 'D'), 'e')
      chunk = [line[len*(i-1)+1:len*i] for i = 1 : npl]
      if typ == Float64
        chunk = map(standardize_real, chunk)
      end
      x[npl * (j-1) + 1 : npl * j] = map(conv, chunk)
  end
  rem = mod(n, npl)
  if rem > 0
    line = join(split(uppercase(readline(io)), 'D'), 'e')
    chunk = [line[len*(i-1)+1:len*i] for i = 1 : rem]
    if typ == Float64
      chunk = map(standardize_real, chunk)
    end
    x[end-rem+1 : end] = map(conv, chunk)
  end
  if scale != 0
    x /= 10^scale
  end
  return is_complex ? [Complex64(x[i], x[i+1]) for i = 1 : 2 : n-1] : x
end


# Displaying Rutherford-Boeing matrix instances.

# import Base.show, Base.print
# function show(io :: IO, rb :: RutherfordBoeingData)
#   s  = @sprintf("Rutherford-Boeing data %s of type %s\n", rb.meta.key, rb.meta.mxtype)
#   s *= @sprintf("%d rows, %d cols, %d nonzeros\n", rb.meta.nrow, rb.meta.ncol, rb.meta.nnzero)
#   s *= @sprintf("%d right-hand sides, %d guesses, %d solutions\n",
#                 size(rb.rhs,2), size(rb.guess,2), size(rb.sol,2))
#   print(io, s)
# end

# function print(io :: IO, rb :: RutherfordBoeingData)
#   @printf("Rutherford-Boeing data %s\n", rb.meta.key)
#   @printf("%s\n", rb.meta.title)
#   @printf("%d rows, %d cols, %d nonzeros\n", rb.meta.nrow, rb.meta.ncol, rb.meta.nnzero)
#   if rb.meta.mxtype[1] == 'P'
#     dtype = "pattern only"
#   elseif rb.meta.mxtype[1] == 'R'
#     dtype = "real"
#   else
#     dtype = "complex"
#   end
#   herm = rb.meta.hermitian ? "hermitian" : "non-hermitian"
#   assm = rb.meta.assembled ? "assembled" : "elemental"
#   @printf("(%s, %s, %s)\n", dtype, herm, assm)
#   @printf("%d right-hand sides, %d guesses, %d solutions\n",
#            size(rb.rhs,2), size(rb.guess,2), size(rb.sol,2))
#   if size(rb.rhs,2) > 0
#     display(rb.rhs'); @printf("\n")
#   end
#   if size(rb.guess,2) > 0
#     display(rb.guess'); @printf("\n")
#   end
#   if size(rb.sol,2) > 0
#     display(rb.sol'); @printf("\n")
#   end
# end

end