-- Hoogle documentation, generated by Haddock
-- See Hoogle, http://www.haskell.org/hoogle/


-- | Sparse matrix linear-equation solver
--   
--   Sparse matrix linear-equation solver, using the conjugate gradient
--   algorithm. Note that the technique only applies to matrices that are
--   symmetric and positive-definite. See
--   <a>http://en.wikipedia.org/wiki/Conjugate_gradient_method</a> for
--   details.
--   
--   The conjugate gradient method can handle very large sparse matrices,
--   where direct methods (such as LU decomposition) are way too expensive
--   to be useful in practice. Such large sparse matrices arise naturally
--   in many engineering problems, such as in ASIC placement algorithms and
--   when solving partial differential equations.
@package conjugateGradient
@version 1.4


-- | (The linear equation solver library is hosted at
--   <a>http://github.com/LeventErkok/conjugateGradient</a>. Comments, bug
--   reports, and patches are always welcome.)
--   
--   Sparse matrix linear-equation solver, using the conjugate gradient
--   algorithm. Note that the technique only applies to matrices that are:
--   
--   <ul>
--   <li>Symmetric</li>
--   <li>Positive-definite</li>
--   </ul>
--   
--   See <a>http://en.wikipedia.org/wiki/Conjugate_gradient_method</a> for
--   details.
--   
--   The conjugate gradient method can handle very large sparse matrices,
--   where direct methods (such as LU decomposition) are way too expensive
--   to be useful in practice. Such large sparse matrices arise naturally
--   in many engineering problems, such as in ASIC placement algorithms and
--   when solving partial differential equations.
--   
--   Here's an example usage, for the simple system:
--   
--   <pre>
--   4x +  y = 1
--    x + 3y = 2
--   </pre>
--   
--   <pre>
--   &gt;&gt;&gt; import Data.IntMap
--   
--   &gt;&gt;&gt; import System.Random
--   
--   &gt;&gt;&gt; let a = (2, fromList [(0, fromList [(0, 4), (1, 1)]), (1, fromList [(0, 1), (1, 3)])]) :: SM Double
--   
--   &gt;&gt;&gt; let b = fromList [(0, 1), (1, 2)] :: SV Double
--   
--   &gt;&gt;&gt; let g = mkStdGen 12345
--   
--   &gt;&gt;&gt; let x = solveCG g a b
--   
--   &gt;&gt;&gt; putStrLn $ showSolution 4 a b x
--         A       |   x    =   b   
--   --------------+----------------
--   4.0000 1.0000 | 0.0909 = 1.0000
--   1.0000 3.0000 | 0.6364 = 2.0000
--   </pre>
module Math.ConjugateGradient

-- | A sparse vector containing elements of type <tt>a</tt>. For our
--   purposes, the elements will be either <a>Float</a>s or <a>Double</a>s.
--   Only the indices that contain non-<tt>0</tt> elements should be given
--   for efficiency purposes. (Nothing will break if you put in elements
--   that are <tt>0</tt>'s, it's just not as efficient.)
type SV a = IntMap a

-- | A sparse matrix is an int-map containing sparse row-vectors:
--   
--   <ul>
--   <li>The first element, <tt>n</tt>, is the number of rows in the
--   matrix, including those with all <tt>0</tt> elements.</li>
--   <li>The matrix is implicitly assumed to be <tt>nxn</tt>, indexed by
--   keys <tt>(0, 0)</tt> to <tt>(n-1, n-1)</tt>.</li>
--   <li>When constructing a sparse-matrix, only put in rows that have a
--   non-<tt>0</tt> element in them for efficiency.</li>
--   <li>Note that you have to give all the non-0 elements: Even though the
--   matrix must be symmetric for the algorithm to work, the matrix should
--   contain all the non-<tt>0</tt> elements, not just the upper (or the
--   lower)-triangle.</li>
--   <li>Make sure the keys of the int-map is a subset of <tt>[0 ..
--   n-1]</tt>, both for the row-indices and the indices of the vectors
--   representing the sparse-rows.</li>
--   </ul>
type SM a = (Int, IntMap (SV a))

-- | Multiply a sparse-vector by a scalar.
sMulSV :: RealFloat a => a -> SV a -> SV a

-- | Multiply a sparse-matrix by a scalar.
sMulSM :: RealFloat a => a -> SM a -> SM a

-- | Add two sparse vectors.
addSV :: RealFloat a => SV a -> SV a -> SV a

-- | Subtract two sparse vectors.
subSV :: RealFloat a => SV a -> SV a -> SV a

-- | Dot product of two sparse vectors.
dotSV :: RealFloat a => SV a -> SV a -> a

-- | Multiply a sparse matrix (nxn) with a sparse vector (nx1), obtaining a
--   sparse vector (nx1).
mulSMV :: RealFloat a => SM a -> SV a -> SV a

-- | Norm of a sparse vector. (Square-root of its dot-product with itself.)
normSV :: RealFloat a => SV a -> a

-- | Conjugate Gradient Solver for the system <tt>Ax=b</tt>. See:
--   <a>http://en.wikipedia.org/wiki/Conjugate_gradient_method</a>.
--   
--   NB. Assumptions on the input:
--   
--   <ul>
--   <li>The <tt>A</tt> matrix is symmetric and positive definite.</li>
--   <li>The indices start from <tt>0</tt> and go consecutively up-to
--   <tt>n-1</tt>. (Only non-<tt>0</tt> value/row indices has to be
--   present, of course.)</li>
--   </ul>
--   
--   For efficiency reasons, we do not check for either property. (If these
--   assumptions are violated, the algorithm will still produce a result,
--   but not the one you expected!)
--   
--   We perform either <tt>10^6</tt> iterations of the Conjugate-Gradient
--   algorithm, or until the error factor is less than <tt>1e-10</tt>. The
--   error factor is defined as the difference of the norm of the current
--   solution from the last one, as we go through the iteration. See
--   <a>http://en.wikipedia.org/wiki/Conjugate_gradient_method#Convergence_properties_of_the_conjugate_gradient_method</a>
--   for a discussion on the convergence properties of this algorithm.
--   
--   The solver can throw an error if it does not converge by <tt>10^6</tt>
--   iterations. This is typically an indication that the input matrix is
--   not symmetric positive definite.
solveCG :: (RandomGen g, RealFloat a, Random a) => g -> SM a -> SV a -> SV a

-- | Display a solution in a human-readable form. Needless to say, only use
--   this method when the system is small enough to fit nicely on the
--   screen.
showSolution :: RealFloat a => Int -> SM a -> SV a -> SV a -> String