content: \section{Cg}
  
  
     ======================================================
     =============== A LINEAR SYSTEM SOLVER ===============
     ======================================================
  
  
  
     This file implements a function to solve a system of linear
     equations using a variant of the conjugate gradient method.
  
     This method produces a series of approximate solutions along
     with the residual vector and other information associated with
     the solution at that point.
  
  
  \begin{code}
  
  module Cg(cgiters,Cg_state (..),showcg_state) 
          where
  
  import Matrix (Vector , Matrix , 
                 Block_list , Col_pos , Row_pos ,
                 vsub,vdot,svmult,vadd,mvmult,norm)
  
  import AbsDensematrix 
  
  import Utils
  
  \end{code}
  
      Note:
            AbsDensematrix is required by hbc.
            Utils provides the 'rjustify' and 'spaces' functions
  
  
  
          ============================================
       ==================== CG ===========================
          ============================================
  
  \begin{code}
  
  data Cg_state = Cg_stateC Vector Vector Vector Vector Int
  type Precond_function = Matrix -> Vector -> Vector
  type Scale_function = Matrix -> Vector -> (Matrix,Vector)
  
  \end{code}
  
  
  
    Cg_state is an approximate solution vector along with its
     associated relevent information.
  
  
           State vars   x      r     p      q     count
  
             x = solution vector
             r = residual vector
             p = ?
             q = ?
            count = number of solution guesses preceding this one
  
  
     cgiters, given a scaling function and a preconditioning function,
     and a system of linear equations to solve (a matrix & vector)
     return an infinite list of solutions along with the residual
     vector and other stuff for that solution.
  
  \begin{code}
  
  cgiters :: Scale_function -> Precond_function -> Matrix -> Vector -> [Cg_state]
  cgiters scale precond a0 b0
     = iterate cgiter (Cg_stateC x0 r0 p0 q0 0)
       where
          x0 = b0 `vsub` b0
          (a,b) = scale a0 b0
          r0 = b
          p0 = precondition r0
          q0 = a `mvmult` p0 
          precondition = precond a
          cgiter (Cg_stateC x r p q cnt) 
  	     = (Cg_stateC x' r' p' q' cnt')       
                 where 
                    qq = q `vdot` q
                    cnt' = cnt + 1
                    x' = x `vadd` (alpha `svmult` p)
                    r' = r `vsub` (alpha `svmult` q)
                    alpha = rq / qq
                    rq = r `vdot` q
                    p'' = precondition r'
                    q'' = a `mvmult` p'' 
                    beta  = qp / qq
                    qp = q `vdot` q'' 
                    p' = p'' `vsub` (beta `svmult` p)
                    q' = q'' `vsub` (beta `svmult` q) 
           
  
  
  showcg_state :: Vector -> Cg_state -> [Char]
  showcg_state soln (Cg_stateC x r p q cnt)
     = rjustify  5 (show cnt) ++ (spaces 4) ++
       rjustify 25 (show (norm err)) ++ (spaces 3) ++
       rjustify 25 (show (norm r)) ++ "\n"
       where
          err = vsub x soln
          alpha = if qq /= 0 then rq / qq
                 else  0
          rq = vdot r q
          qq = vdot q q
  
  \end{code}
  
  