conjugateGradient 1.2 → 1.3

raw patch · 4 files changed

+37/−15 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• INSTALL +1/−1
• Math/ConjugateGradient.hs +14/−8
• RELEASENOTES +19/−3
• conjugateGradient.cabal +3/−3

@@ -1,3 +1,3 @@ The ConjugateGradient library can be installed using haskell from hackage: -     cabal install ConjugateGradient+     cabal install conjugateGradient

@@ -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 -- <http://en.wikipedia.org/wiki/Conjugate_gradient_method#Convergence_properties_of_the_conjugate_gradient_method> -- 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: <http://en.wikipedia.org/wiki/Conjugate_gradient_method#Example_code_in_Matlab>-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

@@ -1,7 +1,23 @@-Hackage: <http://hackage.haskell.org/package/ConjugateGradient>-GitHub:  <http://github.com/LeventErkok/ConjugateGradient>+Hackage: <http://hackage.haskell.org/package/conjugateGradient>+GitHub:  <http://github.com/LeventErkok/conjugateGradient> -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

@@ -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