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