1 files changed, 75 insertions, 11 deletions

diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index 5e2ea84..50372f1 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -3,11 +3,14 @@
 
 module Numeric.LinearAlgebra.Util.CG(
     cgSolve, cgSolve',
-    CGMat, CGState(..), R, V
+    CGState(..), R, V
 ) where
 
 import Data.Packed.Numeric
+import Numeric.Sparse
 import Numeric.Vector()
+import Numeric.LinearAlgebra.Algorithms(linearSolveLS, relativeError, NormType(..))
+import Control.Arrow((***))
 
 {-
 import Util.Misc(debug, debugMat)
@@ -51,7 +54,7 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
     rdx = norm2 dx / max 1 (norm2 x)
 
 conjugrad
-  :: (Transposable m, Contraction m V V)
+  :: (Transposable m mt, Contraction m V V, Contraction mt V V)
   => Bool -> m -> V -> V -> R -> R -> [CGState]
 conjugrad sym a b = solveG (tr a ◇) (a ◇) (cg sym) b
 
@@ -82,27 +85,88 @@ takeUntil q xs = a++ take 1 b
   where
     (a,b) = break q xs
 
-class (Transposable m, Contraction m V V) => CGMat m
-
 cgSolve
-  :: CGMat m
-  => Bool     -- ^ is symmetric
-  -> m        -- ^ coefficient matrix
+  :: Bool          -- ^ is symmetric
+  -> GMatrix       -- ^ coefficient matrix
   -> Vector Double -- ^ right-hand side
-  -> Vector Double        -- ^ solution
+  -> Vector Double -- ^ solution
 cgSolve sym a b  = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0
   where
     n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double))
 
 cgSolve'
-  :: CGMat m
-  => Bool      -- ^ symmetric
+  :: Bool      -- ^ symmetric
   -> R         -- ^ relative tolerance for the residual (e.g. 1E-4)
   -> R         -- ^ relative tolerance for δx (e.g. 1E-3)
   -> Int       -- ^ maximum number of iterations
-  -> m         -- ^ coefficient matrix
+  -> GMatrix   -- ^ coefficient matrix
   -> V         -- ^ initial solution
   -> V         -- ^ right-hand side
   -> [CGState] -- ^ solution
 cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es
 
+
+--------------------------------------------------------------------------------
+
+instance Testable GMatrix
+  where
+    checkT _ = (ok,info)
+      where
+        sma = convo2 20 3
+        x1 = vect [1..20]
+        x2 = vect [1..40]
+        sm = mkSparse sma
+        dm = toDense sma
+
+        s1 = sm !#> x1
+        d1 = dm #> x1
+
+        s2 = tr sm !#> x2
+        d2 = tr dm #> x2
+
+        sdia = mkDiagR 40 20 (vect [1..10])
+        s3 =    sdia !#> x1
+        s4 = tr sdia !#> x2
+        ddia = diagRect 0 (vect [1..10])  40 20
+        d3 = ddia #> x1
+        d4 = tr ddia #> x2
+
+        v = testb 40
+        s5 = cgSolve False sm v
+        d5 = denseSolve dm v
+
+        info = do
+            print sm
+            disp (toDense sma)
+            print s1; print d1
+            print s2; print d2
+            print s3; print d3
+            print s4; print d4
+            print s5; print d5
+            print $ relativeError Infinity s5 d5
+
+        ok = s1==d1
+          && s2==d2
+          && s3==d3
+          && s4==d4
+          && relativeError Infinity s5 d5 < 1E-10
+
+        disp = putStr . dispf 2
+
+        vect = fromList :: [Double] -> Vector Double
+
+        convomat :: Int -> Int -> AssocMatrix
+        convomat n k = [ ((i,j `mod` n),1) | i<-[0..n-1], j <- [i..i+k-1]]
+
+        convo2 :: Int -> Int -> AssocMatrix
+        convo2 n k = m1 ++ m2
+          where
+            m1 = convomat n k
+            m2 = map (((+n) *** id) *** id) m1
+            
+        testb n = vect $ take n $ cycle ([0..10]++[9,8..1])
+        
+        denseSolve a = flatten . linearSolveLS a . asColumn
+
+        -- mkDiag v = mkDiagR (dim v) (dim v) v
+