\section{Input}
Represent a sparse matrix where each "element" of the
matrix is a dense block
Roger L. Wainwright 6-8-90
================================================================
TEST DATA INPUT FUNCTIONS
These functions create test data sets for the linear
system solver.
================================================================
Notes:
gmat n provides a large sparse matrix for data set n.
rhside n provides the b vector for data set n.
soln_vect n provides the x vector for data set n.
wells n provides the length of wells in (gmat n) split
a_easy n provides some "easy" a matrix.
a_hard n provides some "hard" a matrix.
x1 n provides a vector of all 1's for (a_easy n)
AbsDensematrix is imported to keep hbc happy.
mkbigvec and mksparse contain some yucky kludges.
(r,c,b) = (r-1,c-1,..) to make indices zero origin.
\begin{code}
module Input
(gmat,rhside,soln_vect,wells,split,blm',mkbigvec,
mksparse,a_easy,a_hard,x1) where
import List (transpose)
import Matrix
import AbsDensematrix
import Utils
mkbigvec :: [[Float]] -> Vector
mkbigvec v = mkvector (map mkvec v)
mksparse :: [[(Int,Int,[[Float]])]] -> Matrix
mksparse m
= mkmatrix (map (map f) m)
where
f (r,c,b) = (r-1,c-1,mkblock b)
\end{code}
#######################################################################
################ GCOMP DATA INPUT ROUTINES ############################
Note:
readmtx, split, splitmtx was
written by Dan Vasicek and copied from his directory 6/22/90
with minor modifications.
In the code that follows "abc" replaces force(ReadFile x)
as a temporary measure(irm) except in myread in which it
replaces:
f (system ("cat "++fname))
where
f (a,b,c) = a
\begin{code}
directory :: String
directory = "./Data/"
file1 :: Int -> String
file1 n = directory ++ "gcomp.mat" ++ (show n) ++ ".one"
file2a :: Int -> String
file2a n = directory ++ "gcomp.mat" ++ (show n) ++ ".twoa"
file2b :: Int -> String
file2b n = directory ++ "gcomp.mat" ++ (show n) ++ ".twob"
file2c :: Int -> String
file2c n = directory ++ "gcomp.mat" ++ (show n) ++ ".twoc"
file3 :: Int -> String
file3 n = directory ++ "gcomp.mat" ++ (show n) ++ ".three"
file4 :: Int -> String
file4 n = directory ++ "gcomp.mat" ++ (show n) ++ ".four"
file5 :: Int -> String
file5 n = directory ++ "gcomp.mat" ++ (show 2) ++ ".soln"
file6 :: Int -> String
file6 n = directory ++ "gcomp.mat" ++ (show n) ++ ".well"
file :: String
file = "gconem.bmtx3" -- original file
resn = 11
resm = 21
resnm = resn * resm
readmtx rtype x
= map (map numval) (map f y)
where
y = if (rtype == 1) then drop 1 (lines ("abc"))
else
if (rtype == 2) then init (lines ("abc"))
else lines ("abc")
f y = [take 15 (drop 1 y ) ,
take 15 (drop 16 y ) ,
take 15 (drop 31 y ) ,
take 15 (drop 46 y ) ]
readsoln x
= map (map numval) (map f y)
where y = lines ("abc")
f y = [take 15 (drop 8 y ) ,
take 15 (drop 24 y ) ,
take 15 (drop 40 y ) ,
take 15 (drop 56 y ) ]
myreadmtx rtype x
= map (map numval) (map f y)
where y = if (rtype == 1) then drop 1 (lines (myread x))
else
if (rtype == 2) then init (lines (myread x))
else lines (myread x)
f y = [take 15 (drop 1 y ) ,
take 15 (drop 16 y ) ,
take 15 (drop 31 y ) ,
take 15 (drop 46 y ) ]
myread fname
= "abc"
readsolnvect x
= collect_by4 (concat z)
where z = map (map numval) (map f y)
y = lines ("abc")
f y = [take 21 (drop 11 y )]
collect_by4 [] = []
collect_by4 (a:b:c:d:xs)
= [a,b,c,d]: (collect_by4 xs)
\end{code}
******* TEST VARIOUS OPERATIONS ************************
file1 contains a header line
file4 contains an END line
input1' looks like: [ [4 values] [4 values] ... [4 values] ]
input1 looks like: [blocks] as shown below
[ [ [4 values] [4 values] [4values] [4values] ]
[ [ ] [ ] [ ] [ ] ]
. .
. .
[ [ ] [ ] [ ] [ ] ] ]
\begin{code}
split m [] = []
split m xs = (take m xs):(split m (drop m xs))
splitmtx n m xs = split n (split m xs)
\end{code}
************************* READ FILES ****************************
\begin{code}
input1' n = readmtx 1 (file1 n)
input1 n = split 4 (input1' n)
input2a' n = readmtx 0 (file2a n)
input2a n = split 4 (input2a' n)
input2b' n = readmtx 0 (file2b n)
input2b n = split 4 (input2b' n)
input2c' n = readmtx 0 (file2c n)
input2c n = split 4 (input2c' n)
input3' n = readmtx 0 (file3 n)
input3 n = split 4 (input3' n)
input4' n = readmtx 2 (file4 n)
input5' n = readsolnvect (file5 n)
\end{code}
********************** FILE 1 *******************************
\begin{code}
lowerband' n = take (resnm-resn) (drop resn (input1 n))
lowerband n = map transpose (lowerband' n)
block_lowerband n = convert_lowerband (lowerband' n)
\end{code}
********************** FILE 2a *******************************
\begin{code}
ldiagband' n = drop 1 (input2a n)
ldiagband n = map transpose (ldiagband n)
block_ldiagband n = convert_ldiagband (ldiagband' n)
\end{code}
********************** FILE 2b *******************************
\begin{code}
mdiagband' = input2b
mdiagband n = map transpose (mdiagband' n)
block_mdiagband n = convert_mdiagband (mdiagband' n)
\end{code}
********************** FILE 2c *******************************
\begin{code}
udiagband' n = init(input2c n)
udiagband n = map transpose (udiagband' n)
block_udiagband n = convert_udiagband (udiagband' n)
\end{code}
********************** FILE 3 *******************************
\begin{code}
upperband' n = take (resnm-resn) (input3 n)
upperband n = map transpose (upperband' n)
block_upperband n = convert_upperband (upperband' n)
\end{code}
********************** FILE 4 *******************************
\begin{code}
rhside n = mkbigvec (input4' n)
\end{code}
********************** FILE 5 *******************************
\begin{code}
soln_vect n = mkbigvec (input5' n)
\end{code}
******* END TEST VARIOUS OPERATIONS *********************
\begin{code}
convert_lowerband :: [[[Float]]] -> Block_list
convert_lowerband b
= [(i+resn,i,head (drop(i-1) b) ) | i <-[1..length b]]
convert_ldiagband :: [[[Float]]] -> Block_list
convert_ldiagband b
= [(i+1,i,head (drop(i-1)b)) | i <-[1..length b]]
convert_mdiagband :: [[[Float]]] -> Block_list
convert_mdiagband b
= [(i,i,head (drop(i-1)b)) | i <-[1..length b]]
convert_udiagband :: [[[Float]]] -> Block_list
convert_udiagband b
= [(i,i+1,head (drop(i-1)b))| i <-[1..length b]]
convert_upperband :: [[[Float]]] -> Block_list
convert_upperband b
= [(i,i+resn, head (drop(i-1) b)) | i <-[1..length b]]
\end{code}
------------------ more testing ------------------------------------------
Note:
In make_lmat(below) that
Block_list :: [(Int,Int,[[Float]])]
\begin{code}
wells :: Int -> Int
wells 6 = 1
wells 7 = 3
wells 8 = 7
wells n = 0
gvect :: Vector
gvect = mkbigvec (rep resnm [1,1,1,1])
gmat n = make_lmat (block_lowerband n)
(block_ldiagband n) (block_mdiagband n) (block_udiagband n)
(block_upperband n)
gvect0 :: Vector
gvect0 = mkbigvec (rep resnm [0,0,0,0])
make_lmat :: Block_list -> Block_list -> Block_list ->
Block_list -> Block_list -> Matrix
make_lmat lower ldiag mdiag udiag upper
= mksparse [row i | i<- [1..resnm]]
where
row i = combine i lower ldiag mdiag udiag upper
combine i low ld md ud up =
if (i==1) then (take 1 md) ++ (take 1 ud) ++ (take 1 up)
else if (i <=resn) then diag3 ++ upper
else if (i >(resnm - resn)) then lower ++ diag3
else if (i == resnm) then
[last low] ++ [last ld] ++ [last md]
else lower ++ diag3 ++ upper
where
diag3 = (take 1 (drop (i-2) ld )) ++
(take 1 (drop (i-1) md )) ++
(take 1 (drop (i-1) ud ))
upper = take 1(drop(i-1) up)
lower = take 1(drop (i-resn-1) low)
\end{code}
===========================================================
=================== DEBUGGING CONSTANTS ===================
===========================================================
\begin{code}
x0 size = mkbigvec [ [0,0,0,0] | i<-[1..size*size]]
x1 size = mkbigvec [ [1,1,1,1] | i<-[1..size*size]]
a_easy size
= if (size > 1) then blm size
else mksparse [[(1,1,[[11,-1, 0, 0],
[-1,11,-1, 0],
[ 0,-1,11,-1],
[ 0, 0,-1,11]])]]
a_hard size
= if size > 1 then blm_hard size
else error "ummm. system size to small for this model."
off_block = [[-1,0,0,0],
[0,-1,0,0],
[0,0,-1,0],
[0,0,0,-1]]
d_block = [[11,-1, 0, 0],
[-1,11,-1, 0],
[ 0,-1,11,-1],
[ 0, 0,-1,11]]
hard_off_block = (map (map (/90))
[[-1,-2,-3,-4],
[-4,-1,-2,-3],
[-3,-4,-1,-2],
[-2,-3,-4,-1]])
hard_d_block = (map (map (/90))
[[90,-2,-3,-4],
[-4,90,-2,-3],
[-3,-4,90,-2],
[-2,-3,-4,40]])
\end{code}
-----------------------------------------------------------
-----------------------------------------------------------
CONSTRUCT A LIST MATRIX
for a specific pattern using off_block and d_block for data.
Note we assume 4x4 blocks in this pattern.
br (build row) will build a block-row
blm (build list matrix) will build a list matrix
\begin{code}
blm :: Int -> Matrix
blm n = mksparse (blm' n)
blm' :: Int -> [[(Int,Int,[[Float]])]]
blm' n = [row i | i<- [1..(n*n)]]
where row i = br i n
br :: Int -> Int -> [(Int,Int,[[Float]])]
br rw n =
if ((rw == 1) || (rw == (n*n-n+1))) then
[block i| i<- [1,2,n+1,n+2]]
else if ((rw == n) || (rw == n*n)) then [block i | i<-[n-1,n,2*n-1,2*n]]
else if ((rw `mod` n) == 1) then [block i | i<-[1,2,n+1,n+2,2*n+1,2*n+2]]
else if ((rw `mod` n) == 0) then [block i | i<-[n-1,n,2*n-1,2*n,3*n-1,3*n]]
else if ((rw >1) && (rw <n)) || ((rw >(n*n-n+1)) && (rw < n*n)) then [block i | i<-[m,1+m,2+m,n+m,n+1+m,n+2+m]]
else [block i | i<-[m,1+m,2+m,n+m,n+1+m,n+2+m,
2*n+m,2*n+1+m,2*n+2+m]]
where
block i = if (rw == (i+offset)) then (rw,i+offset,d_block)
else (rw,i+offset,off_block)
m = (rw `mod` n) -1
offset = if (rw <= n) then 0
else (((rw-1) `div` n) -1) * n
blm_hard :: Int -> Matrix
blm_hard n
= mksparse hard
where
hard = map (map f) (blm' n)
f (r,c,b) = if (r/=c) then (r,c,hard_off_block)
else (r,c,hard_d_block)
\end{code}
======================================================
Include wells in a given matrix
(Copied from Roger Wainwright and modified)
=======================================================
Note:
include_wells :: matrix -> [[num]] -> [[num]] -> [[num]] -> [[num]] -> matrix
include_well :: matrix -> [num]-> [num]-> [num]-> [num] -> matrix
include_row_well :: matrix -> [num] -> [num] -> matrix
include_col_well :: matrix -> [num] -> [num] -> matrix
include_dw :: matrix -> [num] -> matrix
append_row :: matrix -> (num, [num]) -> matri
In append_row the line:
colnum = length m
is required to append row wells before
col wells.
\begin{code}
include_wells m (a:as) (b:bs) (c:cs) (d:ds)
= if (as == []) then include_well m a b c d
else include_wells m' as bs cs ds
where m' = include_well m a b c d
include_well m well_pos wellg wellh welldw
= include_dw m'' welldw
where m'' = include_col_well m' well_pos wellg
m' = include_row_well m well_pos wellh
include_row_well m well_pos wellh
= m'
where m' = m ++ [new_row]
new_row = map2 f well_pos (split 4 wellh)
f x y = (row_num, x,[y])
row_num = length m +1
include_col_well m well_pos wellg
= foldl append_row m x
where x = zip2 well_pos (split 4 wellg)
include_dw m welldw = append_row m (length m,welldw)
append_row m (rownum, dat)
= m'
where m' = (take (rownum-1) m) ++
[appendrow] ++
drop rownum m
appendrow = row ++ [(rownum,colnum,data')]
data' = split 1 dat
row = m !! (rownum-1)
colnum = length m
numval :: String -> Float
numval x = (read x)::Float
\end{code}
----------------------------------------------------------------
---------- THAT'S ALL FOLKS -------------------------------------
-----------------------------------------------------------------
|
