initial interface to MKL sparse solver

5 files changed, 166 insertions, 0 deletions

diff --git a/packages/sparse/LICENSE b/packages/sparse/LICENSE
new file mode 100644
index 0000000..b86076a
--- /dev/null
+++ b/packages/sparse/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 2014 Alberto Ruiz
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    * Neither the name of the <organization> nor the
+      names of its contributors may be used to endorse or promote products
+      derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
diff --git a/packages/sparse/Setup.lhs b/packages/sparse/Setup.lhs
new file mode 100644
index 0000000..6b32049
--- /dev/null
+++ b/packages/sparse/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/packages/sparse/hmatrix-sparse.cabal b/packages/sparse/hmatrix-sparse.cabal
new file mode 100644
index 0000000..877ab1f
--- /dev/null
+++ b/packages/sparse/hmatrix-sparse.cabal
@@ -0,0 +1,39 @@
+Name:               hmatrix-sparse
+Version:            0.1.0
+License:            BSD3
+License-file:       LICENSE
+Author:             Alberto Ruiz
+Maintainer:         Alberto Ruiz <aruiz@um.es>
+Stability:          experimental
+Homepage:           https://github.com/albertoruiz/hmatrix
+Synopsis:           Sparse linear solver
+Description:        Interface to MKL direct sparse linear solver
+
+-- cabal install --extra-include-dirs=$MKL --extra-lib-dirs=$MKL
+
+Category:           Math
+tested-with:        GHC ==7.8
+
+cabal-version:      >=1.6
+build-type:         Simple
+
+
+library
+    Build-Depends:      base, hmatrix
+
+    hs-source-dirs:     src
+
+    Exposed-modules:    Numeric.LinearAlgebra.Sparse
+
+    ghc-options:        -Wall
+
+    c-sources:          src/Numeric/LinearAlgebra/sparse.c
+
+    cc-options:         -O4 -msse2 -Wall
+
+    extra-libraries:    mkl_intel mkl_sequential mkl_core
+
+source-repository head
+    type:     git
+    location: https://github.com/albertoruiz/hmatrix
+
diff --git a/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
new file mode 100644
index 0000000..ccf28b7
--- /dev/null
+++ b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
@@ -0,0 +1,32 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE RecordWildCards #-}
+
+
+
+module Numeric.LinearAlgebra.Sparse (
+    dss
+) where
+
+import Foreign.C.Types(CInt(..))
+import Data.Packed.Development
+import System.IO.Unsafe(unsafePerformIO)
+import Foreign(Ptr)
+import Numeric.HMatrix
+import Text.Printf
+import Numeric.LinearAlgebra.Util((~!~))
+
+
+type IV t = CInt -> Ptr CInt   -> t
+type  V t = CInt -> Ptr Double -> t
+type SMxV = V (IV (IV (V (V (IO CInt)))))
+
+dss :: CSR -> Vector Double -> Vector Double
+dss CSR{..} b = unsafePerformIO $ do
+    size b /= csrNRows ~!~ printf "dss: incorrect sizes: (%d,%d) x %d" csrNRows csrNCols (size b)
+    r <- createVector csrNCols
+    app5 c_dss vec csrVals vec csrCols vec csrRows vec b vec r "dss"
+    return r
+
+foreign import ccall unsafe "dss"
+  c_dss :: SMxV
+
diff --git a/packages/sparse/src/Numeric/LinearAlgebra/sparse.c b/packages/sparse/src/Numeric/LinearAlgebra/sparse.c
new file mode 100644
index 0000000..b1e257a
--- /dev/null
+++ b/packages/sparse/src/Numeric/LinearAlgebra/sparse.c
@@ -0,0 +1,65 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "mkl_dss.h"
+#include "mkl_types.h"
+#include "mkl_spblas.h"
+
+#define KIVEC(A) int A##n, const int*A##p
+#define KDVEC(A) int A##n, const double*A##p
+#define DVEC(A) int A##n, double*A##p
+#define OK return 0;
+
+
+void check_error(int error)
+{
+        if(error != MKL_DSS_SUCCESS) {
+                printf ("Solver returned error code %d\n", error);
+                exit (1);
+        }
+}
+
+int dss(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
+    MKL_INT nRows = rowsn-1, nCols = rn, nNonZeros = valsn, nRhs = 1;
+    MKL_INT *rowIndex = (MKL_INT*) rowsp;
+    MKL_INT *columns = (MKL_INT*) colsp;
+    double *values = (double*) valsp;
+    _DOUBLE_PRECISION_t *rhs = (_DOUBLE_PRECISION_t*) xp;
+//  _DOUBLE_PRECISION_t *obtrhs = (_DOUBLE_PRECISION_t*) malloc((nCols)*sizeof(_DOUBLE_PRECISION_t));
+    _DOUBLE_PRECISION_t *solValues = (_DOUBLE_PRECISION_t*) rp;
+
+    _MKL_DSS_HANDLE_t handle;
+    _INTEGER_t error;
+//      _CHARACTER_t *uplo;
+    MKL_INT opt;
+
+    opt = MKL_DSS_DEFAULTS;
+    error = dss_create(handle, opt);
+    check_error(error);
+
+    opt = MKL_DSS_NON_SYMMETRIC;
+    error = dss_define_structure(handle, opt, rowIndex, nRows, nCols, columns, nNonZeros);
+    check_error(error);
+
+    opt = MKL_DSS_DEFAULTS;
+    error = dss_reorder(handle, opt, 0);
+    check_error(error);
+
+    opt = MKL_DSS_INDEFINITE;
+    error = dss_factor_real(handle, opt, values);
+    check_error(error);
+
+    int j;
+    for (j = 0; j < nCols; j++) {
+            solValues[j] = 0.0;
+    }
+
+    // Solve system
+    opt = MKL_DSS_REFINEMENT_ON;
+    error = dss_solve_real(handle, opt, rhs, nRhs, solValues);
+    check_error(error);
+
+    OK
+}