Created
October 17, 2014 17:21
-
-
Save dpo/4fe9287d6c01d759064b to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment