eigen 1.2.1 → 1.2.2

raw patch · 12 files changed

+383/−141 lines, 12 filesdep ~basedep ~primitivedep ~vectorPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, primitive, vector

API changes (from Hackage documentation)

Files

• Data/Eigen/Internal.hs +7/−0
• Data/Eigen/LA.hs +79/−35
• Data/Eigen/Matrix.hs +75/−16
• Data/Eigen/Matrix/Mutable.hs +2/−2
• cbits/eigen-proxy.cpp +107/−17
• cbits/eigen-proxy.h +23/−2
• eigen.cabal +9/−4
• test/rank.hs +16/−0
• test/regression.hs +49/−0
• test/regression1.hs +0/−49
• test/solve.hs +16/−0
• test/solve1.hs +0/−16

@@ -34,16 +34,23 @@ foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO () foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt +foreign import ccall "eigen-proxy.h eigen_random"      c_random      :: Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_add"         c_add         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_sub"         c_sub         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_mul"         c_mul         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_diagonal"    c_diagonal    :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_transpose"   c_transpose   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_inverse"     c_inverse     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_adjoint"     c_adjoint     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_conjugate"   c_conjugate   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_normalize"   c_normalize   :: Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_sum"         c_sum         :: Ptr CDouble -> CInt -> CInt -> IO CDouble+foreign import ccall "eigen-proxy.h eigen_prod"        c_prod        :: Ptr CDouble -> CInt -> CInt -> IO CDouble+foreign import ccall "eigen-proxy.h eigen_mean"        c_mean        :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_norm"        c_norm        :: Ptr CDouble -> CInt -> CInt -> IO CDouble+foreign import ccall "eigen-proxy.h eigen_trace"       c_trace       :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_squaredNorm" c_squaredNorm :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_blueNorm"    c_blueNorm    :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm   :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr CDouble -> CInt -> CInt -> IO CDouble+

@@ -55,6 +55,9 @@ module Data.Eigen.LA (     -- * Basic linear solving     Decomposition(..),+    rank,+    kernel,+    image,     solve,     relativeError,     -- * Multiple linear regression@@ -66,29 +69,34 @@ import Foreign.C.String import Foreign.Storable import Foreign.Marshal.Alloc+import qualified Foreign.Concurrent as FC import Control.Applicative import Data.Eigen.Matrix-import Data.Eigen.Internal+import qualified Data.Eigen.Internal as I+import qualified Data.Eigen.Matrix.Mutable as M import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VSM -foreign import ccall "eigen-proxy.h eigen_solve" c_solve :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_solve"    c_solve  :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_rank"     c_rank   :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_image"    c_image  :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_kernel"   c_kernel :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h free"           c_free   :: Ptr a -> IO ()    {- | @-Decomposition           Requirements on the matrix          Speed   Accuracy+Decomposition           Requirements on the matrix          Speed   Accuracy  Rank  Kernel  Image  -PartialPivLU            Invertible                          ++      +-FullPivLU               None                                -       +++-HouseholderQR           None                                ++      +-ColPivHouseholderQR     None                                +       ++-FullPivHouseholderQR    None                                -       +++-LLT                     Positive definite                   +++     +-LDLT                    Positive or negative semidefinite   +++     ++-JacobiSVD               None                                -       ++++PartialPivLU            Invertible                          ++      +         -     -       -+FullPivLU               None                                -       +++       +     +       ++HouseholderQR           None                                ++      +         -     -       -+ColPivHouseholderQR     None                                +       ++        +     -       -+FullPivHouseholderQR    None                                -       +++       +     -       -+LLT                     Positive definite                   +++     +         -     -       -+LDLT                    Positive or negative semidefinite   +++     ++        -     -       -+JacobiSVD               None                                -       +++       +     -       -  The best way to do least squares solving for square matrices is with a SVD decomposition (JacobiSVD) @@@ -112,35 +120,71 @@     -- | Two-sided Jacobi SVD decomposition of a rectangular matrix.     | JacobiSVD deriving (Show, Enum) +rank :: Decomposition -> Matrix -> Int+rank d m = I.performIO $+    unsafeWith m $ \vals rows cols ->+        alloca $ \pr -> do+            I.call $ c_rank (I.cast $ fromEnum d) pr vals rows cols+            I.cast <$> peek pr++-- | Return matrix whose columns form a basis of the null-space of @A@+kernel :: Decomposition -> Matrix -> Matrix+kernel d m1 = I.performIO $+    alloca $ \pvals ->+    alloca $ \prows ->+    alloca $ \pcols ->+        unsafeWith m1 $ \vals1 rows1 cols1 -> do+            I.call $ c_kernel (I.cast $ fromEnum d)+                pvals prows pcols+                vals1 rows1 cols1+            vals <- peek pvals+            rows <- I.cast <$> peek prows+            cols <- I.cast <$> peek pcols+            fp <- FC.newForeignPtr vals $ c_free vals+            return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols+++-- | Return a matrix whose columns form a basis of the column-space of @A@+image :: Decomposition -> Matrix -> Matrix+image d m1 = I.performIO $ +    alloca $ \pvals ->+    alloca $ \prows ->+    alloca $ \pcols ->+        unsafeWith m1 $ \vals1 rows1 cols1 -> do+            I.call $ c_image (I.cast $ fromEnum d)+                pvals prows pcols+                vals1 rows1 cols1+            vals <- peek pvals+            rows <- I.cast <$> peek prows+            cols <- I.cast <$> peek pcols+            fp <- FC.newForeignPtr vals $ c_free vals+            return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols+ -- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@ solve :: Decomposition -> Matrix -> Matrix -> Matrix-solve d a b = performIO $ do-    let-        cols = 1-        rows = m_cols a-    vals <- VSM.new (rows * cols)-    VSM.unsafeWith vals $ \px ->-        VS.unsafeWith (m_vals a) $ \pa ->-            VS.unsafeWith (m_vals b) $ \pb ->-                call $ c_solve (cast $ fromEnum d)-                    px (cast rows) (cast cols)-                    pa (cast $ m_rows a) (cast $ m_cols a)-                    pb (cast $ m_rows b) (cast $ m_cols b)-    Matrix rows cols <$> VS.unsafeFreeze vals-+solve d a b = I.performIO $ do+    x <- M.new (m_cols a) 1+    M.unsafeWith x $ \x_vals x_rows x_cols ->+        unsafeWith a $ \a_vals a_rows a_cols ->+            unsafeWith b $ \b_vals b_rows b_cols ->+                I.call $ c_solve (I.cast $ fromEnum d)+                    x_vals x_rows x_cols+                    a_vals a_rows a_cols+                    b_vals b_rows b_cols+    unsafeFreeze x  -- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm relativeError :: Matrix -> Matrix -> Matrix -> Double-relativeError x a b = performIO $ do-    VS.unsafeWith (m_vals x) $ \px ->-        VS.unsafeWith (m_vals a) $ \pa ->-            VS.unsafeWith (m_vals b) $ \pb ->+relativeError x a b = I.performIO $+    unsafeWith x $ \x_vals x_rows x_cols ->+        unsafeWith a $ \a_vals a_rows a_cols ->+            unsafeWith b $ \b_vals b_rows b_cols ->                 alloca $ \pe -> do-                    call $ c_relativeError pe-                        px (cast $ m_rows x) (cast $ m_cols x)-                        pa (cast $ m_rows a) (cast $ m_cols a)-                        pb (cast $ m_rows b) (cast $ m_cols b)-                    cast <$> peek pe+                    I.call $ c_relativeError pe+                        x_vals x_rows x_cols+                        a_vals a_rows a_cols+                        b_vals b_rows b_cols+                    I.cast <$> peek pe  {- | [(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition

@@ -15,14 +15,13 @@     ones,     identity,     constant,+    random,     -- * Accessing matrix data     cols,     rows,     (!),     coeff,     unsafeCoeff,-    minCoeff,-    maxCoeff,     col,     row,     block,@@ -31,18 +30,31 @@     leftCols,     rightCols,     -- * Matrix properties+    sum,+    prod,+    mean,+    minCoeff,+    maxCoeff,+    trace,     norm,     squaredNorm,+    blueNorm,+    hypotNorm,     determinant,+    -- * Boolean reductions+    all,+    any,+    count,     -- * Matrix operations     add,     sub,     mul,     -- * Matrix transformations+    diagonal,+    transpose,     inverse,     adjoint,     conjugate,-    transpose,     normalize,     modify,     -- * Mutable matrices@@ -53,7 +65,7 @@     unsafeWith, ) where -import Prelude hiding (null)+import Prelude hiding (null, sum, all, any) import Data.List (intercalate) import Data.Tuple import Foreign.Ptr@@ -69,7 +81,7 @@ import qualified Data.Eigen.Internal as I import qualified Data.Eigen.Matrix.Mutable as M --- | Matrix class to be used in pure computations, uses column major memory layout+-- | Matrix to be used in pure computations, uses column major memory layout data Matrix = Matrix {     m_rows :: Int,     m_cols :: Int,@@ -123,6 +135,13 @@         VSM.write vm (n * size + n) 1     return vm +-- | The random matrix of a given size+random :: Int -> Int -> IO Matrix+random rows cols = do+    m <- M.new rows cols+    I.call $ M.unsafeWith m I.c_random+    unsafeFreeze m+ -- | Number of rows for the matrix rows :: Matrix -> Int rows = m_rows@@ -218,6 +237,34 @@             VSM.write vals (col * rows + row) (I.cast $ f row col)     return vals +-- | The sum of all coefficients of the matrix+sum :: Matrix -> Double+sum = _prop I.c_sum++-- | The product of all coefficients of the matrix+prod :: Matrix -> Double+prod = _prop I.c_prod++-- | The mean of all coefficients of the matrix+mean :: Matrix -> Double+mean = _prop I.c_prod++-- | The trace of a matrix is the sum of the diagonal coefficients and can also be computed as sum (diagonal m)+trace :: Matrix -> Double+trace = _prop I.c_trace++-- | Returns true if all of the coefficients in a given matrix evaluate to true+all :: (Double -> Bool) -> Matrix -> Bool+all f = VS.all (f . I.cast) . m_vals++-- | Returns true if at least one of the coefficients in a given matrix evaluates to true+any :: (Double -> Bool) -> Matrix -> Bool+any f = VS.any (f . I.cast) . m_vals++-- | Returns the number of coefficients in a given matrix that evaluate to true+count :: (Double -> Bool) -> Matrix -> Int+count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . m_vals+ -- | For vectors, the l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the square root of the sum of the square of all the matrix entries. For vectors, this is also equals to the square root of the dot product of this with itself. norm :: Matrix -> Double norm = _prop I.c_norm@@ -226,6 +273,14 @@ squaredNorm :: Matrix -> Double squaredNorm = _prop I.c_squaredNorm +-- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.+blueNorm :: Matrix -> Double+blueNorm = _prop I.c_blueNorm++-- | The l2 norm of the matrix avoiding undeflow and overflow. This version use a concatenation of hypot calls, and it is very slow.+hypotNorm :: Matrix -> Double+hypotNorm = _prop I.c_hypotNorm+ -- | The determinant of the matrix determinant :: Matrix -> Double determinant m@@ -253,9 +308,13 @@     | m_cols m1 == m_rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.c_mul m1 m2     | otherwise = error "Matrix.mul: number of columns for lhs matrix should be the same as number of rows for rhs matrix" +-- | Diagonal of the matrix+diagonal :: Matrix -> Matrix+diagonal = _unop swap I.c_diagonal+ {- | Inverse of the matrix -For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses class 'PartialPivLU'+For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition -} inverse :: Matrix -> Matrix inverse m@Matrix{..}@@ -303,31 +362,31 @@ unsafeThaw Matrix{..} = VS.unsafeThaw m_vals >>= return . M.MMatrix m_rows m_cols  -- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.-unsafeWith :: Matrix -> (CInt -> CInt -> Ptr CDouble -> IO a) -> IO a+unsafeWith :: Matrix -> (Ptr CDouble -> CInt -> CInt -> IO a) -> IO a unsafeWith m@Matrix{..} f     | not (valid m) = fail "matrix layout is invalid"-    | otherwise = VS.unsafeWith m_vals $ f (I.cast m_rows) (I.cast m_cols)+    | otherwise = VS.unsafeWith m_vals $ \p -> f p (I.cast m_rows) (I.cast m_cols)  _prop :: (Ptr CDouble -> CInt -> CInt -> IO CDouble) -> Matrix -> Double-_prop f m = I.performIO $ unsafeWith m $ \rows cols vals -> I.cast <$> f vals rows cols+_prop f m = I.cast $ I.performIO $ unsafeWith m f  _binop :: ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix -> Matrix _binop f g m1 m2 = I.performIO $ do-    mm <- uncurry M.new $ f (dims m1) (dims m2)-    M.unsafeWith mm $ \rows0 cols0 vals0 ->-        unsafeWith m1 $ \rows1 cols1 vals1 ->-            unsafeWith m2 $ \rows2 cols2 vals2 ->+    m0 <- uncurry M.new $ f (dims m1) (dims m2)+    M.unsafeWith m0 $ \vals0 rows0 cols0 ->+        unsafeWith m1 $ \vals1 rows1 cols1 ->+            unsafeWith m2 $ \vals2 rows2 cols2 ->                 I.call $ g                     vals0 rows0 cols0                     vals1 rows1 cols1                     vals2 rows2 cols2-    unsafeFreeze mm+    unsafeFreeze m0  _unop :: ((Int,Int) -> (Int,Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix _unop f g m1 = I.performIO $ do     m0 <- uncurry M.new $ f (dims m1)-    M.unsafeWith m0 $ \rows0 cols0 vals0 ->-        unsafeWith m1 $ \rows1 cols1 vals1 ->+    M.unsafeWith m0 $ \vals0 rows0 cols0 ->+        unsafeWith m1 $ \vals1 rows1 cols1 ->             I.call $ g                 vals0 rows0 cols0                 vals1 rows1 cols1

@@ -85,8 +85,8 @@ unsafeWrite MMatrix{..} row col val = VSM.unsafeWrite mm_vals (col * mm_rows + row) (I.cast val)  -- | Pass a pointer to the matrix's data to the IO action. Modifying data through the pointer is unsafe if the matrix could have been frozen before the modification.-unsafeWith :: IOMatrix -> (CInt -> CInt -> Ptr CDouble -> IO a) -> IO a+unsafeWith :: IOMatrix -> (Ptr CDouble -> CInt -> CInt -> IO a) -> IO a unsafeWith mm@MMatrix{..} f     | not (valid mm) = fail "mutable matrix layout is invalid"-    | otherwise = VSM.unsafeWith mm_vals $ \p -> f (I.cast mm_rows) (I.cast mm_cols) p+    | otherwise = VSM.unsafeWith mm_vals $ \p -> f p (I.cast mm_rows) (I.cast mm_cols) 

@@ -21,16 +21,19 @@ }  typedef Map< Matrix<double,Dynamic,Dynamic> > MapMatrix;+typedef Map< Matrix<double,1,Dynamic> > MapVector;+typedef Map< Vector3d > MapVector3d;+typedef Map< Vector4d > MapVector4d;  extern "C" {  const char* eigen_add(-    double* data, int rows, int cols, +    double* data, int rows, int cols,     const double* data1, int rows1, int cols1,     const double* data2, int rows2, int cols2) {     GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) + MapMatrix(data2, rows2, cols2);;+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) + MapMatrix(data2, rows2, cols2);     GUARD_END } @@ -40,7 +43,7 @@     const double* data2, int rows2, int cols2) {     GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) - MapMatrix(data2, rows2, cols2);;+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) - MapMatrix(data2, rows2, cols2);     GUARD_END } @@ -50,7 +53,7 @@     const double* data2, int rows2, int cols2) {     GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) * MapMatrix(data2, rows2, cols2);;+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) * MapMatrix(data2, rows2, cols2);     GUARD_END } @@ -67,7 +70,6 @@ const char* eigen_inverse(double* data, int rows, int cols, const double* data1, int rows1, int cols1) {     GUARD_START-    puts("inverse");     MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).inverse();     GUARD_END }@@ -75,7 +77,6 @@ const char* eigen_adjoint(double* data, int rows, int cols, const double* data1, int rows1, int cols1) {     GUARD_START-    puts("adjoint");     MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).adjoint();     GUARD_END }@@ -83,32 +84,121 @@ const char* eigen_conjugate(double* data, int rows, int cols, const double* data1, int rows1, int cols1) {     GUARD_START-    puts("conjugate");     MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).conjugate();     GUARD_END } +const char* eigen_diagonal(double* data, int rows, int cols, const double* data1, int rows1, int cols1)+{+    GUARD_START+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).diagonal();+    GUARD_END+}+ const char* eigen_transpose(double* data, int rows, int cols, const double* data1, int rows1, int cols1) {     GUARD_START-    puts("transpose");     MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).transpose();     GUARD_END } +const char* eigen_dot(double* retval, +    const double* data1, int size1,+    const double* data2, int size2)+{+    GUARD_START+    *retval = MapVector(data1, size1).dot(MapVector(data2, size2));+    GUARD_END+}++const char* eigen_cross(+    double* data,+    const double* data1,+    const double* data2)+{+    GUARD_START+    MapVector3d retval(data);+    retval = MapVector3d(data1).cross(MapVector3d(data2));+    GUARD_END+}++const char* eigen_cross3(+    double* data,+    const double* data1,+    const double* data2)+{+    GUARD_START+    MapVector4d retval(data);+    retval = MapVector4d(data1).cross3(MapVector4d(data2));+    GUARD_END+}+ const char* eigen_normalize(double* data, int rows, int cols) {     GUARD_START-    puts("normalize");     MapMatrix(data, rows, cols).normalize();     GUARD_END } +const char* eigen_random(double* data, int rows, int cols)+{+    GUARD_START+    MapMatrix(data, rows, cols) = MatrixXd::Random(rows, cols);+    GUARD_END+} +const char* eigen_rank(Decomposition d, int* r, const double* data, int rows, int cols) {+    GUARD_START+    MapMatrix A(data, rows, cols);+    switch (d) {+        case ::FullPivLU:+            *r = A.fullPivLu().rank();+            break;+        case ::ColPivHouseholderQR:+            *r = A.colPivHouseholderQr().rank();+            break;+        case ::FullPivHouseholderQR:+            *r = A.fullPivHouseholderQr().rank();+            break;+        case ::JacobiSVD:+            *r = A.jacobiSvd(ComputeThinU | ComputeThinV).rank();+            break;+        default:+            return strdup("Selected decomposition doesn't support rank revealing.");+    }+    GUARD_END+}++const char* eigen_kernel(Decomposition d, double** data0, int* rows0, int* cols0, const double* data1, int rows1, int cols1) {+    GUARD_START+    if (d != ::FullPivLU)+        return strdup("Selected decomposition doesn't support kernel revealing.");+    MapMatrix A(data1, rows1, cols1);+    MatrixXd B = A.fullPivLu().kernel();+    *rows0 = B.rows();+    *cols0 = B.cols();+    *data0 = (double*)malloc(*rows0 * *cols0 * sizeof(double));+    MapMatrix(*data0, *rows0, *cols0) = B;+    GUARD_END+}++const char* eigen_image(Decomposition d, double** data0, int* rows0, int* cols0, const double* data1, int rows1, int cols1) {+    GUARD_START+    if (d != ::FullPivLU)+        return strdup("Selected decomposition doesn't support image revealing.");+    MapMatrix A(data1, rows1, cols1);+    MatrixXd B = A.fullPivLu().image(A);+    *rows0 = B.rows();+    *cols0 = B.cols();+    *data0 = (double*)malloc(*rows0 * *cols0 * sizeof(double));+    MapMatrix(*data0, *rows0, *cols0) = B;+    GUARD_END+}+ const char* eigen_solve(Decomposition d,-    double* px, int rx, int cx, // x-    const double* pa, int ra, int ca, // A-    const double* pb, int rb, int cb) // b+    double* px, int rx, int cx,+    const double* pa, int ra, int ca,+    const double* pb, int rb, int cb) {     GUARD_START     MapMatrix x(px, rx, cx);@@ -144,16 +234,16 @@ }  -const char* eigen_relativeError(double& e,-    const double* px, int rx, int cx, // x-    const double* pa, int ra, int ca, // A-    const double* pb, int rb, int cb) // b+const char* eigen_relativeError(double* e,+    const double* px, int rx, int cx,+    const double* pa, int ra, int ca,+    const double* pb, int rb, int cb) {     GUARD_START     MapMatrix x(px, rx, cx);     MapMatrix A(pa, ra, ca);     MapMatrix b(pb, rb, cb);-    e = (A*x - b).norm() / b.norm();+    *e = (A*x - b).norm() / b.norm();     GUARD_END } 

@@ -24,8 +24,13 @@ double eigen_mean(const double*, int, int); double eigen_trace(const double*, int, int); double eigen_determinant(const double*, int, int);+const char* eigen_dot(double*, const double*, int, const double*, int);+const char* eigen_cross(double*, const double*, const double*);+const char* eigen_cross3(double*, const double*, const double*);+const char* eigen_diagonal(double*, int, int, const double*, int, int); const char* eigen_transpose(double*, int, int, const double*, int, int); const char* eigen_normalize(double*, int, int);+const char* eigen_random(double*, int, int); const char* eigen_inverse(double*, int, int, const double*, int, int); const char* eigen_conjugate(double*, int, int, const double*, int, int); const char* eigen_adjoint(double*, int, int, const double*, int, int);@@ -34,15 +39,31 @@ void eigen_setNbThreads(int);  enum Decomposition {-	PartialPivLU, FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, LLT, LDLT, JacobiSVD+	PartialPivLU,+	FullPivLU,+	HouseholderQR,+	ColPivHouseholderQR,+	FullPivHouseholderQR,+	LLT,+	LDLT,+	JacobiSVD,+	// SelfAdjointEigenSolver,+	// ComplexEigenSolver,+	// EigenSolver,+	// GeneralizedSelfAdjointEigenSolver }; +const char* eigen_rank(Decomposition d, int*, const double*, int, int);+const char* eigen_kernel(Decomposition d, double**, int*, int*, const double*, int, int);+const char* eigen_image(Decomposition d, double**, int*, int*, const double*, int, int);+ const char* eigen_solve(Decomposition d, // Ax=b  	double*, int, int, // x 	const double*, int, int, // A 	const double*, int, int); // b -const char* eigen_relativeError(double& e, ++const char* eigen_relativeError(double* e,  	const double*, int, int, // x 	const double*, int, int, // A 	const double*, int, int); // b;

@@ -1,5 +1,5 @@ name:           eigen-version:        1.2.1+version:        1.2.2 homepage:       https://github.com/osidorkin/haskell-eigen synopsis:       Eigen C++ library (linear algebra: matrices, vectors, numerical solvers). description:    This module provides Haskell binding for Eigen C++ library.@@ -372,17 +372,22 @@                         vector >= 0.5 && < 0.11,                         primitive >= 0.1 && < 0.6 -    include-dirs:       eigen3+    include-dirs:       eigen3, cbits     c-sources:          cbits/eigen-proxy.cpp     extra-libraries:    stdc++  Test-Suite test-solve     type:               exitcode-stdio-1.0-    main-is:            test/solve1.hs+    main-is:            test/solve.hs     build-depends:      base, primitive, vector, eigen +Test-Suite test-rank+    type:               exitcode-stdio-1.0+    main-is:            test/rank.hs+    build-depends:      base, primitive, vector, eigen+ Test-Suite test-regression     type:               exitcode-stdio-1.0-    main-is:            test/regression1.hs+    main-is:            test/regression.hs     build-depends:      base, primitive, vector, eigen 

@@ -0,0 +1,16 @@+import Data.Eigen.Matrix+import Data.Eigen.LA++main = do+    let a = fromList [[1,2,5],[2,1,4],[3,0,3]]+    putStrLn "Here is the matrix A:"+    print a++    putStrLn "The rank of A is:"+    print $ rank FullPivLU a+   +    putStrLn "Here is a matrix whose columns form a basis of the null-space of A:"+    print $ kernel FullPivLU a++    putStrLn "Here is a matrix whose columns form a basis of the column-space of A:"+    print $ image FullPivLU a

@@ -0,0 +1,49 @@+{-# LANGUAGE RecordWildCards #-}+import Data.Eigen.Matrix as M+import Data.Eigen.LA+import Data.List as L+import Control.Monad++main = do+    let+        a = fromList [+            [1,3.02, 6.89],+            [1,2.01, 5.39],+            [1,2.41, 6.01],+            [1,2.09, 5.55],+            [1,2.58, 6.32]]++        b = fromList $ map return [-4.32,-3.79,-4.01,-3.86,-4.10]++    print a+    print b+    forM_ [FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, JacobiSVD] $ \d -> do+        let x = solve d a b+            e = relativeError x a b+            e' = norm (a*x - b) / norm b+        putStrLn $ replicate 20 '*'+        print d+        print x+        print e+        print e'+    putStrLn "\n-2.34666 - 0.25349 x1 - 0.174965 x2"+    putStrLn "done"++    print $ identity 4+    print $ M.normalize a+    print $ M.transpose a++    let+        a = M.fromList [[0.68,  0.597,  -0.33],[-0.211,  0.823,  0.536],[ 0.566, -0.605, -0.444]]+        b = M.inverse a+    print a+    print b+    print $ a * b+    print $ linearRegression [+            [-4.32, 3.02, 6.89],+            [-3.79, 2.01, 5.39],+            [-4.01, 2.41, 6.01],+            [-3.86, 2.09, 5.55],+            [-4.10, 2.58, 6.32]]++

@@ -1,49 +0,0 @@-{-# LANGUAGE RecordWildCards #-}-import Data.Eigen.Matrix as M-import Data.Eigen.LA-import Data.List as L-import Control.Monad--main = do-    let-        a = fromList [-            [1,3.02, 6.89],-            [1,2.01, 5.39],-            [1,2.41, 6.01],-            [1,2.09, 5.55],-            [1,2.58, 6.32]]--        b = fromList $ map return [-4.32,-3.79,-4.01,-3.86,-4.10]--    print a-    print b-    forM_ [FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, JacobiSVD] $ \d -> do-        let x = solve d a b-            e = relativeError x a b-            e' = norm (a*x - b) / norm b-        putStrLn $ replicate 20 '*'-        print d-        print x-        print e-        print e'-    putStrLn "\n-2.34666 - 0.25349 x1 - 0.174965 x2"-    putStrLn "done"--    print $ identity 4-    print $ M.normalize a-    print $ M.transpose a--    let-        a = M.fromList [[0.68,  0.597,  -0.33],[-0.211,  0.823,  0.536],[ 0.566, -0.605, -0.444]]-        b = M.inverse a-    print a-    print b-    print $ a * b-    print $ linearRegression [-            [-4.32, 3.02, 6.89],-            [-3.79, 2.01, 5.39],-            [-4.01, 2.41, 6.01],-            [-3.86, 2.09, 5.55],-            [-4.10, 2.58, 6.32]]--

@@ -0,0 +1,16 @@+import Data.Eigen.Matrix+import Data.Eigen.LA++main = do+    let+        a = fromList [[1,2,3], [4,5,6], [7,8,10]]+        b = fromList [[3],[3],[4]]+        x = solve ColPivHouseholderQR a b+    putStrLn "Here is the matrix A:"+    print a++    putStrLn "Here is the vector b:"+    print b++    putStrLn "The solution is:"+    print x

@@ -1,16 +0,0 @@-import Data.Eigen.Matrix-import Data.Eigen.LA--main = do-	let-		a = fromList [[1,2,3], [4,5,6], [7,8,10]]-		b = fromList [[3],[3],[4]]-		x = solve ColPivHouseholderQR a b-	putStrLn "Here is the matrix A:"-	print a--	putStrLn "Here is the vector b:"-	print b--	putStrLn "The solution is:"-	print x