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