1 \section{Matlib}
3 =======================================================================
4         FULL FACTORIZATION OF WELL TERMS
6  This file contains all the functions required by the full
7  factorization preconditioning strategy employed in the different
8  iterative linear system solver algorithms
9  Several general purpose matrix/vector operations are also included.
10  -- written by Brian D. Moe (Summer 1990)
11  Converted to Haskell by Ian R McLelland and Cordelia V Hall (Nov 1992)
12  =======================================================================
15 Notes:
17   Scale_function and Precond_function were originally exported
18   but hbc will not allow this, and as they are not used by the
19   only importing module(Test) they have been left out of the
20   export list here(irm)
22   The import of AbsDensematrix is to keep hbc happy.
24 \begin{code}
25 module Matlib (scale, precond, uncondition) where
27 import List (genericLength)
28 import Matrix
29 import AbsDensematrix
32 type Scale_function   = Matrix -> Vector -> (Matrix,Vector)
33 type Precond_function = Matrix -> Vector -> Vector
35 type Block_tuple = (Int,Int,Block)
38 scale       :: Scale_function
39 precond     :: Int -> Precond_function
40 uncondition :: Precond_function
41 \end{code}
44  =================================================
45  ============= SPARSE MATRIX SCALING =============
46  =================================================
48 Notes:
50     bdscl is Block diagonal scaling
51     BDSCL (A = L + D + U)  , V   ->    (L_new, I, U_new) , V_new
54 \begin{code}
55 scale = bdscl
58 uncondition a v = mvmult (diag_matrix a) v
61 bdscl :: Matrix -> Vector -> (Matrix,Vector)
62 bdscl m vs
63    = (m', vs')
64      where
65         vs' = mkvector [ bvecmult (row_factor r) (vsubscript r vs)
66                           | r <- [0..(vsize vs)-1] ]
67         m' = mkmatrix [map scale1 (getrow i m) | i <- [0..(numrows m)-1]]
69         scale1 (r,c,b) = if r==c then (r,c, bident (fst (bsize b)))
70                          else (r,c, bmult (row_factor r) b)
71         row_factor n = if (okindex n diag_inverses) then  diag_inverses !! n
72                        else error "bdscl:dinverse in matrix"
74         diag_inverses = map inverse (diag_blocks m)
75         inverse (r,c,b) = binverse b
76 \end{code}
79      =========================================================
80      ================>>   PRECONDITIONING   <<================
81      =========================================================
84    The preconditioner given here is as specified by Julio Diaz.
85    It has been designed to help condition systems of equations
86    of the type found in oil reservoir models.
88    The basic idea is that the matrix is put through
89    an incomplete LU factorization, where only the diagonal blocks
90    are altered.  The rows and columns which represent wells in the
91    reservoir model (and occupy the last rows & columns) are then
92    treated to full Gaussian elimination.  Once this has been
93    accomplished (by the function "pseudofactor"), back substitution
94    is performed on the vector to be conditioned.  This vector
95    supplies a better guess than might otherwise be used for the
96    direction the solution vector must move.
99   Given a sparse matrix, perform incomplete lu decomposition
100   through the maxrow'th row/col, inclusive.
102   The resultant matrix is the sum of the old lower and upper
103   matrices and the matrix made up of the diagonal blocks which
104   are the D_u's as below.
106       D_u =  inverse [ B_u_u - sum  (B_u_k * D_k * B_k_u)  ]
107                                k<u
111   Block-based Complete LU decomposition.
113   Gauss elimination starting with a given row/column through the
114   rest of the matrix.
116  ======== Notes on the functions to follow ==================
119    precond
121    Given the number of rows/columns which represent well terms,
122    the system matrix and some vector, do the preconditioning on
123    the vector (as described above)
125         ------------------
126    pseudofactor
128    The partial LU factorization as described above.
130    The subfunctions bilu and bclu are also described
131    in more detail above.
132         ------------------
133    colsbefore
135    Given a row of block tuples, return only the columns before the one given.
136         ------------------
137    backsubs
139       (L + D) * D_inv * (U + D) * p = r
141   [NOTE: D_inv is the diagonal block matrix where the diagonals are inverses
142        of those in matrix D.  The diagonal blocks come already inverted, since
143        those inverses are already calculated in pseudofactor and the original
144        blocks are not needed here. ]
145        ---------------------------
146     forward m r
148     returns the vector y which satisfies:
150       (L+D) y = r
152         ==>
154       Dy = r - Ly
156   where (L+D) are the lower blocks and the *inverses* of the diagonal blocks
157   of the matrix m.  [So... m = (L+Dinv+U)]
159   This is calculated in block matrix format in this manner:
161       y_i = Dinv_i * [ r_i - sum ( B_i_j * y_j ) ]
162                              j<k
163         -----------------------------
164     backward m y
166     returns the vector p which satisfies:
168       (Dinv * U + I) p = y
170   ==>
171       p = y - (Dinv U p)
173     where U are the upper blocks and Dinv are the diagonal blocks
174     of the matrix m.  [So... m = (L+Dinv+U)]
176     This is calculated in block matrix format in this manner:
178         p_i = y_i - Dinv_i * sum [ B_i_j * p_j ]
179                              j>i
181 =============== End of Notes =======================
188 \begin{code}
189 precond nwells a v
190     = backsubs (pseudofactor firstwell a) v
191     where
192         firstwell = (numrows a) - nwells
195 pseudofactor :: Int -> Matrix -> Matrix
196 pseudofactor n mat
197    = mat'
198      where
199         mat' = mkmatrix [ newrow k | k <- blocknums mat ]
200         newrow k = (bilu k) ++ (bclu k)
201         maxblock = last (blocknums mat)
202         bilu k = if (k>=n) then []
203                  else  map bilu' (colsbefore n (getrow k mat))
204         bilu' (r,c,oldb)
205            = if r/=c then (r,c,oldb)
206              else  (r,c,binverse newb)
207               where
208                  newb = if sumpart == [] then oldb
209                        else  oldb bsub (bsum sumpart)
210                  k=r
211                  sumpart = [ b'uk bmult (b' k k) bmult (b' k u)
212                            | (u,k,b'uk) <- colsbefore k (getrow k mat') ]
213         b  i j = msubscript i j mat
214         b' i j = msubscript i j mat'
216         first_col_in_row r
217            = if getrow r mat == [] then  maxblock
218              else colno (head (getrow r mat))
219         colno (r,c,b) = c
220         firstcols = [first_col_in_row r | r<-[n..maxblock]]
221         firstcol r = firstcols !! (r-n)
222         exist r c
223            = if c>r then exist c r
224              else if c>=(firstcol r) then True
225                    else False
227         bclu r
228            = if (r<n) then [ bclu' r c  | c <- [n..maxblock], exist r c ]
229              else [ bclu' r c  | c <- blocknums mat, exist r c ]
230         bclu' i j
231            = if result == [] then error "pseudofactor"
232              else if i==j then (i,j,binverse (head result))
234              where
235              result =  if i<j then [b' i i] mult ([b i j] sub sumpartU)
236                       else [b i j] sub sumpartL
237              sumpartU = sumblocks [ [bik] mult [b k j]
238                                       | (x,k,bik) <- (getrow i mat), k<i ]
239              sumpartL = sumblocks [ [bik] mult [b' k k] mult [b' k j]
240                                       | (x,k,bik) <- (getrow i mat), k<j ]
244 colsbefore n row = [(r,c,b) | (r,c,b) <- row, c < n ]
248 backsubs m v
249    = v'
250       where
251          v' = ((backward m) . (forward m))  v
255 forward m rs
256    = mkvector [(y k) | k <- (blocknums m) ]
257      where
258         y k = bvecmult (dinv k) (terms k)
259         dinv k = b k k
260         r k = vsubscript k rs
261         terms k = if (sumpart k)==[] then (r k)
262                   else  (r k) vecsub (vecsum (sumpart k))
263         sumpart k = sumparts!!k
264         sumparts = [ [ b_ij bvecmult (y j) | (i,j, b_ij) <- (getrow k m) , j<i ]
265                       | k <- (blocknums m) ]
266         b i j = msubscript i j m
271 backward :: Matrix -> Vector -> [Vec]
272 backward m ys
273    = mkvector ps
274      where
275        ps = [ (p' k) | k <- (blocknums m) ]
276        p' k =  if (terms k) == [] then (y k)
277               else (y k) vecsub ((b k k) bvecmult (vecsum (terms k)))
278        p k = ps !! k
279        b i j = msubscript i j m
280        y k = vsubscript k ys
281        termss  = [ [b_i_j bvecmult (p j) | (i,j, b_i_j) <- (getrow k m), j>i ]
282                    | k <- (blocknums m) ]
283        terms k = termss!!k
284 \end{code}
288   ============================================================
289   |  Matrix Operations which can deal with missing operands  |
290   |      (missing operands are considered zero matrices)     |
291   ============================================================
293 \begin{code}
294 mult [ ] [ ] = []
295 mult [x] [ ] = []
296 mult [ ] [y] = []
297 mult [x] [y] = [bmult x y]
299 sub  [ ] [ ] = []
300 sub  [x] [ ] = [x]
301 sub  [ ] [y] = [bneg y]
302 sub  [x] [y] = [bsub x y]
304 sumblocks xs = if blocks /= [] then [bsum blocks]
305                else []
306                where blocks = concat xs
307 \end{code}
311  ==========================================================================
312  ================ Some general purpose matrix/vector stuff ================
313  ==========================================================================
315 Notes:
317      bsum sums up a bunch of blocks.
319      diag_matrix returns a matrix which consists of the diagonal blocks of
320      a sparse matrix.
322      is_diag determines whether or not a block is on the diagonal.
324      vecsum sums up a bunch of vectors.
326      okindex checks that the number given won't be out of range for the given
327      list. WARNING : WILL HANG IF LIST IS INFINITE.
330 \begin{code}
331 bsum :: [Block] -> Block
332 bsum [] = error "bsum: no blocks to sum"
333 bsum ms = foldl1 badd ms
335 diag_matrix :: Matrix -> Matrix
336 diag_matrix m = mkmatrix [filter is_diag (getrow i m) | i <- blocknums m]
338 diag_blocks :: Matrix -> [Block_tuple]
339 diag_blocks m
340    = if (length(diags) == numrows m) then concat diags
341      else  error "matlib:diag_blocks:  given matrix with missing diagonal block"
342      where
343         diags = [filter is_diag (getrow i m) | i <- blocknums m]
345 is_diag :: Block_tuple -> Bool
346 is_diag (r,c,b) = r==c
349 blocknums :: Matrix -> [Int]
350 blocknums m = [0..(numrows m)-1]
353 vecsum :: [Vec] -> Vec
354 vecsum [] = error "vecsum: no vecs to sum"
355 vecsum ms = foldl1 vecadd ms
357 okindex :: Int -> [a] -> Bool
358 okindex n m = (0<=n)&&(n<=((length m)-1))
360 testmat :: Matrix -> [Char]
361 testmat m
362    = if result == [] then  "Matrix is probably ok"
363      else result
364      where
365         result = (rowsnumbered m) ++ (symmetric m) ++ (sorted m)
367 rowsnumbered :: Matrix -> [Char]
368 rowsnumbered m
369    = concat [ goodrow i (getrow i m) | i <- blocknums m ]
370      where
371         goodrow i row = concat  (map (isrow i) row)
372         isrow i (r,c,b)
373            =if i==r then  []
374             else  ("Row " ++ (show i) ++ " is misnumbered\n")
376 symmetric :: Matrix -> [Char]
377 symmetric m
378    = concat [ symmetric_row (getrow i m) | i <- blocknums m ]
379      where
380         symmetric_row row = concat [ has_corresponding elem | elem <- row ]
381         has_corresponding (r,c,b)
382            = if exists c r then []
383              else "Cannot find corresponding block for " ++ (show (r,c))++"\n"
384         exists r c = (filter (iscol c) (getrow r m)) /= []
385         iscol c (i,j,b) = c==j
387 sorted :: Matrix -> [Char]
388 sorted m
389    = concat [ sortedrow i (getrow i m) | i <- blocknums m ]
390      where
391         sortedrow i row
392            = if row == (sort row) then []
393              else "Row number " ++ (show i) ++ " is not properly sorted.\n"
395 sort :: (Ord a) => [a] -> [a]
396 sort xs = if (n == 1) then xs
397           else merge (sort us) (sort vs)
398           where
399            n = genericLength xs
400            us = take (n div 2) xs
401            vs = drop (n div 2) xs
404 merge :: (Ord a) => [a] -> [a] -> [a]
405 merge [] ys = ys
406 merge (x:xs) [] = (x:xs)
407 merge (x:xs) (y:ys) = if (x <= y) then (x:merge xs (y:ys))
408                     else (y:(merge (x:xs) ys))
409 \end{code}