1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 module Cg(cgiters,Cg_state (..),showcg_state) 21 where 22 23 import Matrix (Vector , Matrix , 24 Block_list , Col_pos , Row_pos , 25 vsub,vdot,svmult,vadd,mvmult,norm) 26 27 import AbsDensematrix 28 29 import Utils 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 data Cg_state = Cg_stateC Vector Vector Vector Vector Int 46 type Precond_function = Matrix -> Vector -> Vector 47 type Scale_function = Matrix -> Vector -> (Matrix,Vector) 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 cgiters :: Scale_function -> Precond_function -> Matrix -> Vector -> [Cg_state] 74 cgiters scale precond a0 b0 75 = iterate cgiter (Cg_stateC x0 r0 p0 q0 0) 76 where 77 x0 = b0 `vsub` b0 78 (a,b) = scale a0 b0 79 r0 = b 80 p0 = precondition r0 81 q0 = a `mvmult` p0 82 precondition = precond a 83 cgiter (Cg_stateC x r p q cnt) 84 = (Cg_stateC x' r' p' q' cnt') 85 where 86 qq = q `vdot` q 87 cnt' = cnt + 1 88 x' = x `vadd` (alpha `svmult` p) 89 r' = r `vsub` (alpha `svmult` q) 90 alpha = rq / qq 91 rq = r `vdot` q 92 p'' = precondition r' 93 q'' = a `mvmult` p'' 94 beta = qp / qq 95 qp = q `vdot` q'' 96 p' = p'' `vsub` (beta `svmult` p) 97 q' = q'' `vsub` (beta `svmult` q) 98 99 100 101 showcg_state :: Vector -> Cg_state -> [Char] 102 showcg_state soln (Cg_stateC x r p q cnt) 103 = rjustify 5 (show cnt) ++ (spaces 4) ++ 104 rjustify 25 (show (norm err)) ++ (spaces 3) ++ 105 rjustify 25 (show (norm r)) ++ "\n" 106 where 107 err = vsub x soln 108 alpha = if qq /= 0 then rq / qq 109 else 0 110 rq = vdot r q 111 qq = vdot q q 112 113 114