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