\section{Matlib}

=======================================================================
        FULL FACTORIZATION OF WELL TERMS

 This file contains all the functions required by the full
 factorization preconditioning strategy employed in the different
 iterative linear system solver algorithms
 Several general purpose matrix/vector operations are also included.
 -- written by Brian D. Moe (Summer 1990)
 Converted to Haskell by Ian R McLelland and Cordelia V Hall (Nov 1992)
 =======================================================================


Notes:

  Scale_function and Precond_function were originally exported
  but hbc will not allow this, and as they are not used by the 
  only importing module(Test) they have been left out of the 
  export list here(irm)

  The import of AbsDensematrix is to keep hbc happy.

\begin{code}
module Matlib (scale, precond, uncondition) where 

import List (genericLength)
import Matrix
import AbsDensematrix 


type Scale_function   = Matrix -> Vector -> (Matrix,Vector)
type Precond_function = Matrix -> Vector -> Vector

type Block_tuple = (Int,Int,Block)


scale       :: Scale_function
precond     :: Int -> Precond_function
uncondition :: Precond_function
\end{code}


 =================================================
 ============= SPARSE MATRIX SCALING =============
 =================================================

Notes:

    bdscl is Block diagonal scaling
    BDSCL (A = L + D + U)  , V   ->    (L_new, I, U_new) , V_new 


\begin{code}
scale = bdscl


uncondition a v = mvmult (diag_matrix a) v


bdscl :: Matrix -> Vector -> (Matrix,Vector)
bdscl m vs
   = (m', vs')
     where
        vs' = mkvector [ bvecmult (row_factor r) (vsubscript r vs)
                          | r <- [0..(vsize vs)-1] ]
        m' = mkmatrix [map scale1 (getrow i m) | i <- [0..(numrows m)-1]]

        scale1 (r,c,b) = if r==c then (r,c, bident (fst (bsize b)))
                         else (r,c, bmult (row_factor r) b)
        row_factor n = if (okindex n diag_inverses) then  diag_inverses !! n
                       else error "bdscl:dinverse in matrix"

	diag_inverses = map inverse (diag_blocks m)
        inverse (r,c,b) = binverse b
\end{code}


     =========================================================
     ================>>   PRECONDITIONING   <<================
     =========================================================


   The preconditioner given here is as specified by Julio Diaz.
   It has been designed to help condition systems of equations
   of the type found in oil reservoir models.

   The basic idea is that the matrix is put through
   an incomplete LU factorization, where only the diagonal blocks
   are altered.  The rows and columns which represent wells in the
   reservoir model (and occupy the last rows & columns) are then
   treated to full Gaussian elimination.  Once this has been 
   accomplished (by the function "pseudofactor"), back substitution
   is performed on the vector to be conditioned.  This vector
   supplies a better guess than might otherwise be used for the
   direction the solution vector must move. 


  Given a sparse matrix, perform incomplete lu decomposition
  through the maxrow'th row/col, inclusive.

  The resultant matrix is the sum of the old lower and upper
  matrices and the matrix made up of the diagonal blocks which
  are the D_u's as below.

      D_u =  inverse [ B_u_u - sum  (B_u_k * D_k * B_k_u)  ]
                               k<u



  Block-based Complete LU decomposition.

  Gauss elimination starting with a given row/column through the
  rest of the matrix.

 ======== Notes on the functions to follow ==================


   precond 

   Given the number of rows/columns which represent well terms,
   the system matrix and some vector, do the preconditioning on
   the vector (as described above)
      
        ------------------
   pseudofactor

   The partial LU factorization as described above.

   The subfunctions bilu and bclu are also described
   in more detail above.
        ------------------
   colsbefore
  
   Given a row of block tuples, return only the columns before the one given.
        ------------------
   backsubs

      (L + D) * D_inv * (U + D) * p = r

  [NOTE: D_inv is the diagonal block matrix where the diagonals are inverses
       of those in matrix D.  The diagonal blocks come already inverted, since
       those inverses are already calculated in pseudofactor and the original
       blocks are not needed here. ]
       ---------------------------
    forward m r

    returns the vector y which satisfies:
 
      (L+D) y = r

        ==>

      Dy = r - Ly

  where (L+D) are the lower blocks and the *inverses* of the diagonal blocks
  of the matrix m.  [So... m = (L+Dinv+U)]

  This is calculated in block matrix format in this manner:

      y_i = Dinv_i * [ r_i - sum ( B_i_j * y_j ) ]
                             j<k
        -----------------------------
    backward m y

    returns the vector p which satisfies:

      (Dinv * U + I) p = y

  ==>
      p = y - (Dinv U p)

    where U are the upper blocks and Dinv are the diagonal blocks
    of the matrix m.  [So... m = (L+Dinv+U)]

    This is calculated in block matrix format in this manner:
 
        p_i = y_i - Dinv_i * sum [ B_i_j * p_j ]
                             j>i

=============== End of Notes =======================



         


\begin{code}
precond nwells a v
    = backsubs (pseudofactor firstwell a) v
    where
        firstwell = (numrows a) - nwells
  

pseudofactor :: Int -> Matrix -> Matrix
pseudofactor n mat
   = mat'
     where
        mat' = mkmatrix [ newrow k | k <- blocknums mat ]
        newrow k = (bilu k) ++ (bclu k)
        maxblock = last (blocknums mat)
        bilu k = if (k>=n) then []                 
                 else  map bilu' (colsbefore n (getrow k mat))
        bilu' (r,c,oldb)
           = if r/=c then (r,c,oldb)
             else  (r,c,binverse newb)
              where
                 newb = if sumpart == [] then oldb
                       else  oldb `bsub` (bsum sumpart)
                 k=r
                 sumpart = [ b'uk `bmult` (b' k k) `bmult` (b' k u)
                           | (u,k,b'uk) <- colsbefore k (getrow k mat') ]
        b  i j = msubscript i j mat
        b' i j = msubscript i j mat'

        first_col_in_row r
           = if getrow r mat == [] then  maxblock                   
             else colno (head (getrow r mat))
        colno (r,c,b) = c
        firstcols = [first_col_in_row r | r<-[n..maxblock]]
        firstcol r = firstcols !! (r-n)
        exist r c
           = if c>r then exist c r
             else if c>=(firstcol r) then True
                   else False

        bclu r
           = if (r<n) then [ bclu' r c  | c <- [n..maxblock], exist r c ]
             else [ bclu' r c  | c <- blocknums mat, exist r c ]
        bclu' i j
           = if result == [] then error "pseudofactor"
             else if i==j then (i,j,binverse (head result))
             else (i,j,head result)
             where
             result =  if i<j then [b' i i] `mult` ([b i j] `sub` sumpartU)  
                      else [b i j] `sub` sumpartL   
             sumpartU = sumblocks [ [bik] `mult` [b k j]
                                      | (x,k,bik) <- (getrow i mat), k<i ]
             sumpartL = sumblocks [ [bik] `mult` [b' k k] `mult` [b' k j]
                                      | (x,k,bik) <- (getrow i mat), k<j ]



colsbefore n row = [(r,c,b) | (r,c,b) <- row, c < n ]


  
backsubs m v
   = v'
      where
         v' = ((backward m) . (forward m))  v



forward m rs
   = mkvector [(y k) | k <- (blocknums m) ] 
     where
        y k = bvecmult (dinv k) (terms k)
        dinv k = b k k
        r k = vsubscript k rs
        terms k = if (sumpart k)==[] then (r k)                              
                  else  (r k) `vecsub` (vecsum (sumpart k))
        sumpart k = sumparts!!k
        sumparts = [ [ b_ij `bvecmult` (y j) | (i,j, b_ij) <- (getrow k m) , j<i ]
                      | k <- (blocknums m) ]
        b i j = msubscript i j m


 

backward :: Matrix -> Vector -> [Vec]
backward m ys
   = mkvector ps
     where
       ps = [ (p' k) | k <- (blocknums m) ]
       p' k =  if (terms k) == [] then (y k)
              else (y k) `vecsub` ((b k k) `bvecmult` (vecsum (terms k)))
       p k = ps !! k
       b i j = msubscript i j m
       y k = vsubscript k ys
       termss  = [ [b_i_j `bvecmult` (p j) | (i,j, b_i_j) <- (getrow k m), j>i ]
                   | k <- (blocknums m) ]
       terms k = termss!!k
\end{code}



  ============================================================
  |  Matrix Operations which can deal with missing operands  |
  |      (missing operands are considered zero matrices)     |
  ============================================================

\begin{code}
mult [ ] [ ] = []
mult [x] [ ] = []
mult [ ] [y] = []
mult [x] [y] = [bmult x y]

sub  [ ] [ ] = []
sub  [x] [ ] = [x]
sub  [ ] [y] = [bneg y]
sub  [x] [y] = [bsub x y]

sumblocks xs = if blocks /= [] then [bsum blocks]
               else []           
               where blocks = concat xs
\end{code}



 ==========================================================================
 ================ Some general purpose matrix/vector stuff ================
 ==========================================================================

Notes:

     bsum sums up a bunch of blocks.

     diag_matrix returns a matrix which consists of the diagonal blocks of 
     a sparse matrix.

     is_diag determines whether or not a block is on the diagonal.

     vecsum sums up a bunch of vectors.

     okindex checks that the number given won't be out of range for the given
     list. WARNING : WILL HANG IF LIST IS INFINITE.


\begin{code}
bsum :: [Block] -> Block
bsum [] = error "bsum: no blocks to sum"
bsum ms = foldl1 badd ms

diag_matrix :: Matrix -> Matrix
diag_matrix m = mkmatrix [filter is_diag (getrow i m) | i <- blocknums m]

diag_blocks :: Matrix -> [Block_tuple]
diag_blocks m
   = if (length(diags) == numrows m) then concat diags
     else  error "matlib:diag_blocks:  given matrix with missing diagonal block"
     where
	diags = [filter is_diag (getrow i m) | i <- blocknums m]

is_diag :: Block_tuple -> Bool
is_diag (r,c,b) = r==c


blocknums :: Matrix -> [Int]
blocknums m = [0..(numrows m)-1]


vecsum :: [Vec] -> Vec
vecsum [] = error "vecsum: no vecs to sum"
vecsum ms = foldl1 vecadd ms

okindex :: Int -> [a] -> Bool
okindex n m = (0<=n)&&(n<=((length m)-1))

testmat :: Matrix -> [Char]
testmat m
   = if result == [] then  "Matrix is probably ok"
     else result
     where
	result = (rowsnumbered m) ++ (symmetric m) ++ (sorted m)

rowsnumbered :: Matrix -> [Char]
rowsnumbered m
   = concat [ goodrow i (getrow i m) | i <- blocknums m ]
     where
	goodrow i row = concat  (map (isrow i) row)
	isrow i (r,c,b)
	   =if i==r then  []
	    else  ("Row " ++ (show i) ++ " is misnumbered\n")

symmetric :: Matrix -> [Char]
symmetric m
   = concat [ symmetric_row (getrow i m) | i <- blocknums m ]
     where
	symmetric_row row = concat [ has_corresponding elem | elem <- row ]
	has_corresponding (r,c,b)
	   = if exists c r then []     
	     else "Cannot find corresponding block for " ++ (show (r,c))++"\n"
	exists r c = (filter (iscol c) (getrow r m)) /= []
	iscol c (i,j,b) = c==j

sorted :: Matrix -> [Char]
sorted m
   = concat [ sortedrow i (getrow i m) | i <- blocknums m ]
     where
	sortedrow i row
	   = if row == (sort row) then []
	     else "Row number " ++ (show i) ++ " is not properly sorted.\n"

sort :: (Ord a) => [a] -> [a]
sort xs = if (n == 1) then xs
          else merge (sort us) (sort vs)
          where 
           n = genericLength xs
           us = take (n `div` 2) xs
           vs = drop (n `div` 2) xs
     

merge :: (Ord a) => [a] -> [a] -> [a]
merge [] ys = ys
merge (x:xs) [] = (x:xs)
merge (x:xs) (y:ys) = if (x <= y) then (x:merge xs (y:ys))
                    else (y:(merge (x:xs) ys))
\end{code}