diff --git a/INSTALL b/INSTALL --- a/INSTALL +++ b/INSTALL @@ -1,3 +1,3 @@ The ConjugateGradient library can be installed using haskell from hackage: - cabal install ConjugateGradient + cabal install conjugateGradient diff --git a/Math/ConjugateGradient.hs b/Math/ConjugateGradient.hs --- a/Math/ConjugateGradient.hs +++ b/Math/ConjugateGradient.hs @@ -31,10 +31,11 @@ -- @ -- -- >>> import Data.IntMap +-- >>> import System.Random -- >>> let a = (2, fromList [(0, fromList [(0, 4), (1, 1)]), (1, fromList [(0, 1), (1, 3)])]) :: SM Double -- >>> let b = fromList [(0, 1), (1, 2)] :: SV Double -- >>> let g = mkStdGen 12345 --- >>> let (_, x) = solveCG g a b +-- >>> let x = solveCG g a b -- >>> putStrLn $ showSolution 4 a b x -- A | x = b -- --------------+---------------- @@ -57,8 +58,8 @@ import Data.List (intercalate) import Data.Maybe (fromMaybe) import qualified Data.IntMap as IM (IntMap, lookup, map, unionWith, intersectionWith, fold, fromList) -import System.Random -import Numeric +import System.Random (Random, RandomGen, randomRs) +import Numeric (showFFloat) -- | A sparse vector containing elements of type 'a'. (For our purposes, the elements will be either 'Float's or 'Double's.) type SV a = IM.IntMap a @@ -70,7 +71,9 @@ -- * The matrix is implicitly assumed to be @nxn@, indexed by keys @(0, 0)@ to @(n-1, n-1)@. -- -- * When constructing a sparse-matrix, only put in rows that have a non-@0@ element in them for efficiency. --- (Nothing will break if you put in rows that have all @0@'s in them, it's just not as efficient.) +-- (Nothing will break if you put in rows that have all @0@'s in them, it's just not as efficient.) 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-@0@ elements, not just the upper (or the lower)-triangle. -- -- * Make sure the keys of the int-map is a subset of @[0 .. n-1]@, both for the row-indices and the indices of the vectors representing the sparse-rows. type SM a = (Int, IM.IntMap (SV a)) @@ -132,25 +135,28 @@ -- the current solution from the last one, as we go through the iteration. See -- -- for a discussion on the convergence properties of this algorithm. +-- +-- The solver can throw an error if it does not converge by @10^6@ iterations. This is typically an indication +-- that the input matrix is not symmetric positive definite. solveCG :: (RandomGen g, RealFloat a, Random a) => g -- ^ The seed for the random-number generator. -> SM a -- ^ The @A@ sparse matrix (@nxn@). -> SV a -- ^ The @b@ sparse vector (@nx1@). - -> (a, SV a) -- ^ The final error factor, and the @x@ sparse matrix (@nx1@), such that @Ax = b@. + -> SV a -- ^ The @x@ sparse matrix (@nx1@), such that @Ax = b@. solveCG g a@(n, _) b = cg a b x0 where rs = take n (randomRs (0, 1) g) x0 = IM.fromList [p | p@(_, j) <- zip [0..] rs, j /= 0] -- | The Conjugate-gradient algorithm. Our implementation closely follows the -- one given here: -cg :: RealFloat a => SM a -> SV a -> SV a -> (a, SV a) +cg :: RealFloat a => SM a -> SV a -> SV a -> SV a cg a b x0 = cgIter (1000000 :: Int) (norm r0) r0 r0 x0 where r0 = b `subSV` (a `mulSMV` x0) - cgIter 0 eps _ _ x = (eps, x) + cgIter 0 _ _ _ _ = error "Conjugate Gradient: No convergence after 10^6 iterations. Make sure the input matrix is symmetric positive-definite!" cgIter i eps p r x -- Stop if the square of the error is less than 1e-20, i.e., -- if the error itself is less than 1e-10. - | eps' < 1e-20 = (eps', x') + | eps' < 1e-20 = x' | True = cgIter (i-1) eps' p' r' x' where ap = a `mulSMV` p alpha = eps / ap `dotSV` p diff --git a/RELEASENOTES b/RELEASENOTES --- a/RELEASENOTES +++ b/RELEASENOTES @@ -1,7 +1,23 @@ -Hackage: -GitHub: +Hackage: +GitHub: -Latest Hackage released version: 1.1 +Latest Hackage released version: 1.3 + +====================================================================== +Version 1.3, 2013-04-16 + - Instead of returning an error-bound, throw an error if + no convergence is reached after 10^6 iterations. This is + more practical, as returning a result after that many + iterations typically indicates the input matrix is not + symmetric and positive-definite. + - Tighten import lists and the example + - Clarify that the entire matrix should be given: Even though + we assume it's symmetric, the algorithm needs all non-0 elements + to be present; not just the upper (or the lower)-triangle. + +====================================================================== +Version 1.2, 2013-04-15 + - Simplify types, clean-up example. ====================================================================== Version 1.1, 2013-04-15 diff --git a/conjugateGradient.cabal b/conjugateGradient.cabal --- a/conjugateGradient.cabal +++ b/conjugateGradient.cabal @@ -1,5 +1,5 @@ Name: conjugateGradient -Version: 1.2 +Version: 1.3 Category: Math Synopsis: Sparse matrix linear-equation solver Description: Sparse matrix linear-equation solver, using the conjugate gradient algorithm. Note that the @@ -20,8 +20,8 @@ License-file: LICENSE Stability: Experimental Author: Levent Erkok -Homepage: http://github.com/LeventErkok/ConjugateGradient -Bug-reports: http://github.com/LeventErkok/ConjugateGradient/issues +Homepage: http://github.com/LeventErkok/conjugateGradient +Bug-reports: http://github.com/LeventErkok/conjugateGradient/issues Maintainer: Levent Erkok (erkokl@gmail.com) Build-Type: Simple Cabal-Version: >= 1.14