maybe variant of linearSolve

4 files changed, 44 insertions, 22 deletions

diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
index 786fb6d..024e462 100644
--- a/packages/base/src/Numeric/HMatrix.hs
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -152,7 +152,8 @@ import Numeric.LinearAlgebra.Data
 import Numeric.Matrix()
 import Numeric.Vector()
 import Data.Packed.Numeric hiding ((<>))
-import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Algorithms hiding (linearSolve)
+import qualified Numeric.LinearAlgebra.Algorithms as A
 import Numeric.LinearAlgebra.Util
 import Numeric.LinearAlgebra.Random
 import Numeric.Sparse((!#>))
@@ -163,3 +164,6 @@ import Numeric.LinearAlgebra.Util.CG
 (<>) = mXm
 infixr 8 <>
 
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. 
+linearSolve m b = A.mbLinearSolve m b
+
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index bbcc513..7e36978 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -26,6 +26,7 @@ module Numeric.LinearAlgebra.Algorithms (
     Field(),
 -- * Linear Systems
     linearSolve,
+    mbLinearSolve,
     luSolve,
     cholSolve,
     linearSolveLS,
@@ -102,6 +103,7 @@ class (Product t,
     sv'          :: Matrix t -> Vector Double
     luPacked'    :: Matrix t -> (Matrix t, [Int])
     luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
+    mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (Matrix t)
     linearSolve' :: Matrix t -> Matrix t -> Matrix t
     cholSolve'   :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
@@ -125,6 +127,7 @@ instance Field Double where
     luPacked' = luR
     luSolve' (l_u,perm) = lusR l_u perm
     linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
+    mbLinearSolve' = mbLinearSolveR
     cholSolve' = cholSolveR
     linearSolveLS' = linearSolveLSR
     linearSolveSVD' = linearSolveSVDR Nothing
@@ -151,6 +154,7 @@ instance Field (Complex Double) where
     luPacked' = luC
     luSolve' (l_u,perm) = lusC l_u perm
     linearSolve' = linearSolveC
+    mbLinearSolve' = mbLinearSolveC
     cholSolve' = cholSolveC
     linearSolveLS' = linearSolveLSC
     linearSolveSVD' = linearSolveSVDC Nothing
@@ -234,6 +238,10 @@ luSolve = {-# SCC "luSolve" #-} luSolve'
 linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolve = {-# SCC "linearSolve" #-} linearSolve'
 
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. 
+mbLinearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
+mbLinearSolve = {-# SCC "linearSolve" #-} mbLinearSolve'
+
 -- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
 cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
 cholSolve = {-# SCC "cholSolve" #-} cholSolve'
diff --git a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
index 40fef45..b32b67f 100644
--- a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
@@ -17,6 +17,7 @@ module Numeric.LinearAlgebra.LAPACK (
     multiplyR, multiplyC, multiplyF, multiplyQ,
     -- * Linear systems
     linearSolveR, linearSolveC,
+    mbLinearSolveR, mbLinearSolveC,
     lusR, lusC,
     cholSolveR, cholSolveC,
     linearSolveLSR, linearSolveLSC,
@@ -329,8 +330,8 @@ foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM
 foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM
 foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM
 
-linearSolveSQAux f st a b
-    | n1==n2 && n1==r = unsafePerformIO $ do
+linearSolveSQAux g f st a b
+    | n1==n2 && n1==r = unsafePerformIO . g $ do
         s <- createMatrix ColumnMajor r c
         app3 f mat a mat b mat s st
         return s
@@ -342,20 +343,26 @@ linearSolveSQAux f st a b
 
 -- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
 linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
+linearSolveR a b = linearSolveSQAux id dgesv "linearSolveR" (fmat a) (fmat b)
+
+mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
+mbLinearSolveR a b = linearSolveSQAux mbCatch dgesv "linearSolveR" (fmat a) (fmat b)
+
 
 -- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.
 linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
+linearSolveC a b = linearSolveSQAux id zgesv "linearSolveC" (fmat a) (fmat b)
 
+mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
+mbLinearSolveC a b = linearSolveSQAux mbCatch zgesv "linearSolveC" (fmat a) (fmat b)
 
 -- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
 cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b)
+cholSolveR a b = linearSolveSQAux id dpotrs "cholSolveR" (fmat a) (fmat b)
 
 -- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.
 cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b)
+cholSolveC a b = linearSolveSQAux id zpotrs "cholSolveC" (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
 foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM
diff --git a/packages/base/src/Numeric/LinearAlgebra/Real.hs b/packages/base/src/Numeric/LinearAlgebra/Real.hs
index 0e54555..8627084 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Real.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Real.hs
@@ -29,7 +29,7 @@ module Numeric.LinearAlgebra.Real(
     -- * Vector
     R, C,
     vec2, vec3, vec4, (&), (#),
-    vect,
+    vector,
     linspace, range, dim,
     -- * Matrix
     L, Sq, M,
@@ -39,7 +39,7 @@ module Numeric.LinearAlgebra.Real(
     eye,
     diagR, diag,
     blockAt,
-    mat,
+    matrix,
     -- * Products
     (<>),(#>),(<·>),
     -- * Linear Systems
@@ -64,6 +64,7 @@ import Data.Proxy(Proxy)
 import Numeric.LinearAlgebra.Static
 import Text.Printf
 
+
 𝑖 :: Sized ℂ s c => s
 𝑖 = konst i_C
 
@@ -71,7 +72,7 @@ instance forall n . KnownNat n => Show (R n)
   where
     show (ud1 -> v)
       | singleV v = "("++show (v!0)++" :: R "++show d++")"
-      | otherwise   = "(vect"++ drop 8 (show v)++" :: R "++show d++")"
+      | otherwise   = "(vector"++ drop 8 (show v)++" :: R "++show d++")"
       where
         d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
 
@@ -79,7 +80,7 @@ instance forall n . KnownNat n => Show (C n)
   where
     show (C (Dim v))
       | singleV v = "("++show (v!0)++" :: C "++show d++")"
-      | otherwise   = "(fromList"++ drop 8 (show v)++" :: C "++show d++")"
+      | otherwise   = "(vector"++ drop 8 (show v)++" :: C "++show d++")"
       where
         d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
 
@@ -117,8 +118,11 @@ vec3 a b c = R (gvec3 a b c)
 vec4 :: ℝ -> ℝ -> ℝ -> ℝ -> R 4
 vec4 a b c d = R (gvec4 a b c d)
 
-vect :: forall n . KnownNat n => [ℝ] -> R n
-vect xs = R (gvect "R" xs)
+vector :: KnownNat n => [ℝ] -> R n
+vector = fromList
+
+matrix :: (KnownNat m, KnownNat n) => [ℝ] -> L m n
+matrix = fromList
 
 linspace :: forall n . KnownNat n => (ℝ,ℝ) -> R n
 linspace (a,b) = mkR (LA.linspace d (a,b))
@@ -157,7 +161,7 @@ instance forall m n . (KnownNat m, KnownNat n) => Show (L m n)
     show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'
     show (ud2 -> x)
        | singleM x = printf "(%s :: L %d %d)" (show (x `atIndex` (0,0))) m' n'
-       | otherwise = "(mat"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")"
+       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")"
       where
         m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
         n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
@@ -167,7 +171,7 @@ instance forall m n . (KnownNat m, KnownNat n) => Show (M m n)
     show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'
     show (M (Dim (Dim x)))
        | singleM x = printf "(%s :: M %d %d)" (show (x `atIndex` (0,0))) m' n'
-       | otherwise = "(fromList"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")"
+       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")"
       where
         m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
         n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
@@ -191,7 +195,7 @@ instance forall n. KnownNat n => Sized ℝ (R n) (Vector ℝ)
   where
     konst x = mkR (LA.scalar x)
     unwrap = ud1
-    fromList = vect
+    fromList xs = R (gvect "R" xs)
     extract (unwrap -> v)
       | singleV v = LA.konst (v!0) d
       | otherwise = v
@@ -203,7 +207,7 @@ instance forall n. KnownNat n => Sized ℝ (R n) (Vector ℝ)
 instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) (Matrix ℝ)
   where
     konst x = mkL (LA.scalar x)
-    fromList = mat
+    fromList xs = L (gmat "L" xs)
     unwrap = ud2
     extract (isDiag -> Just (z,y,(m',n'))) = diagRect z y m' n'
     extract (unwrap -> a)
@@ -260,8 +264,7 @@ blockAt x r c a = mkL res
 
 
 
-mat :: forall m n . (KnownNat m, KnownNat n) => [ℝ] -> L m n
-mat xs = L (gmat "L" xs)
+
 
 --------------------------------------------------------------------------------
 
@@ -505,8 +508,8 @@ instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
 
 --------------------------------------------------------------------------------
 
-linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> L m n
-linSolve (extract -> a) (extract -> b) = mkL (LA.linearSolve a b)
+linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)
+linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b)
 
 (<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r
 (extract -> a) <\> (extract -> b) = mkL (a LA.<\> b)
@@ -617,7 +620,7 @@ test = (ok,info)
 
     u = vec2 3 5
 
-    𝕧 x = vect [x] :: R 1
+    𝕧 x = vector [x] :: R 1
 
     v = 𝕧 2 & 4 & 7