[project @ 2001-04-04 11:12:09 by sewardj]
[nofib.git] / real / linear / Matlib.lhs
1 \section{Matlib}
2
3 =======================================================================
4         FULL FACTORIZATION OF WELL TERMS
5
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  =======================================================================
13
14
15 Notes:
16
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)
21
22   The import of AbsDensematrix is to keep hbc happy.
23
24 \begin{code}
25 module Matlib (scale, precond, uncondition) where 
26
27 import List (genericLength)
28 import Matrix
29 import AbsDensematrix 
30
31
32 type Scale_function   = Matrix -> Vector -> (Matrix,Vector)
33 type Precond_function = Matrix -> Vector -> Vector
34
35 type Block_tuple = (Int,Int,Block)
36
37
38 scale       :: Scale_function
39 precond     :: Int -> Precond_function
40 uncondition :: Precond_function
41 \end{code}
42
43
44  =================================================
45  ============= SPARSE MATRIX SCALING =============
46  =================================================
47
48 Notes:
49
50     bdscl is Block diagonal scaling
51     BDSCL (A = L + D + U)  , V   ->    (L_new, I, U_new) , V_new 
52
53
54 \begin{code}
55 scale = bdscl
56
57
58 uncondition a v = mvmult (diag_matrix a) v
59
60
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]]
68
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"
73
74         diag_inverses = map inverse (diag_blocks m)
75         inverse (r,c,b) = binverse b
76 \end{code}
77
78
79      =========================================================
80      ================>>   PRECONDITIONING   <<================
81      =========================================================
82
83
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.
87
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. 
97
98
99   Given a sparse matrix, perform incomplete lu decomposition
100   through the maxrow'th row/col, inclusive.
101
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.
105
106       D_u =  inverse [ B_u_u - sum  (B_u_k * D_k * B_k_u)  ]
107                                k<u
108
109
110
111   Block-based Complete LU decomposition.
112
113   Gauss elimination starting with a given row/column through the
114   rest of the matrix.
115
116  ======== Notes on the functions to follow ==================
117
118
119    precond 
120
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)
124       
125         ------------------
126    pseudofactor
127
128    The partial LU factorization as described above.
129
130    The subfunctions bilu and bclu are also described
131    in more detail above.
132         ------------------
133    colsbefore
134   
135    Given a row of block tuples, return only the columns before the one given.
136         ------------------
137    backsubs
138
139       (L + D) * D_inv * (U + D) * p = r
140
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
147
148     returns the vector y which satisfies:
149  
150       (L+D) y = r
151
152         ==>
153
154       Dy = r - Ly
155
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)]
158
159   This is calculated in block matrix format in this manner:
160
161       y_i = Dinv_i * [ r_i - sum ( B_i_j * y_j ) ]
162                              j<k
163         -----------------------------
164     backward m y
165
166     returns the vector p which satisfies:
167
168       (Dinv * U + I) p = y
169
170   ==>
171       p = y - (Dinv U p)
172
173     where U are the upper blocks and Dinv are the diagonal blocks
174     of the matrix m.  [So... m = (L+Dinv+U)]
175
176     This is calculated in block matrix format in this manner:
177  
178         p_i = y_i - Dinv_i * sum [ B_i_j * p_j ]
179                              j>i
180
181 =============== End of Notes =======================
182
183
184
185          
186
187
188 \begin{code}
189 precond nwells a v
190     = backsubs (pseudofactor firstwell a) v
191     where
192         firstwell = (numrows a) - nwells
193   
194
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'
215
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
226
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))
233              else (i,j,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 ]
241
242
243
244 colsbefore n row = [(r,c,b) | (r,c,b) <- row, c < n ]
245
246
247   
248 backsubs m v
249    = v'
250       where
251          v' = ((backward m) . (forward m))  v
252
253
254
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
267
268
269  
270
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}
285
286
287
288   ============================================================
289   |  Matrix Operations which can deal with missing operands  |
290   |      (missing operands are considered zero matrices)     |
291   ============================================================
292
293 \begin{code}
294 mult [ ] [ ] = []
295 mult [x] [ ] = []
296 mult [ ] [y] = []
297 mult [x] [y] = [bmult x y]
298
299 sub  [ ] [ ] = []
300 sub  [x] [ ] = [x]
301 sub  [ ] [y] = [bneg y]
302 sub  [x] [y] = [bsub x y]
303
304 sumblocks xs = if blocks /= [] then [bsum blocks]
305                else []           
306                where blocks = concat xs
307 \end{code}
308
309
310
311  ==========================================================================
312  ================ Some general purpose matrix/vector stuff ================
313  ==========================================================================
314
315 Notes:
316
317      bsum sums up a bunch of blocks.
318
319      diag_matrix returns a matrix which consists of the diagonal blocks of 
320      a sparse matrix.
321
322      is_diag determines whether or not a block is on the diagonal.
323
324      vecsum sums up a bunch of vectors.
325
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.
328
329
330 \begin{code}
331 bsum :: [Block] -> Block
332 bsum [] = error "bsum: no blocks to sum"
333 bsum ms = foldl1 badd ms
334
335 diag_matrix :: Matrix -> Matrix
336 diag_matrix m = mkmatrix [filter is_diag (getrow i m) | i <- blocknums m]
337
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]
344
345 is_diag :: Block_tuple -> Bool
346 is_diag (r,c,b) = r==c
347
348
349 blocknums :: Matrix -> [Int]
350 blocknums m = [0..(numrows m)-1]
351
352
353 vecsum :: [Vec] -> Vec
354 vecsum [] = error "vecsum: no vecs to sum"
355 vecsum ms = foldl1 vecadd ms
356
357 okindex :: Int -> [a] -> Bool
358 okindex n m = (0<=n)&&(n<=((length m)-1))
359
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)
366
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")
375
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
386
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"
394
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
402      
403
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}