hmatrix-sparse (empty) → 0.19.0.0

raw patch · 5 files changed

+175/−0 lines, 5 filesdep +basedep +hmatrixsetup-changed

Dependencies added: base, hmatrix

Files

• LICENSE +26/−0
• Setup.lhs +4/−0
• hmatrix-sparse.cabal +45/−0
• src/Numeric/LinearAlgebra/Sparse.hs +35/−0
• src/Numeric/LinearAlgebra/sparse.c +65/−0

@@ -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 Alberto Ruiz nor the names of other 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.+

@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain

@@ -0,0 +1,45 @@+Name:               hmatrix-sparse+Version:            0.19.0.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<5, hmatrix>=0.16++    hs-source-dirs:     src++    Exposed-modules:    Numeric.LinearAlgebra.Sparse++    ghc-options:        -Wall++    c-sources:          src/Numeric/LinearAlgebra/sparse.c++    cc-options:         -O4 -Wall++    if arch(x86_64)+       cc-options:     -msse2++    if arch(i386)+       cc-options:     -msse2++    extra-libraries:    mkl_intel mkl_sequential mkl_core++source-repository head+    type:     git+    location: https://github.com/albertoruiz/hmatrix+

@@ -0,0 +1,35 @@+{-# LANGUAGE ForeignFunctionInterface #-}+{-# LANGUAGE RecordWildCards #-}++++module Numeric.LinearAlgebra.Sparse (+    dss+) where++import Foreign.C.Types(CInt(..))+import Numeric.LinearAlgebra.Devel+import System.IO.Unsafe(unsafePerformIO)+import Foreign(Ptr)+import Numeric.LinearAlgebra.HMatrix+import Text.Printf+import Control.Monad(when)++(???) :: Bool -> String -> IO ()+infixl 0 ???+c ??? msg = when c (error msg)++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+    c_dss `apply` csrVals `apply` csrCols `apply` csrRows `apply` b `apply` r #|"dss"+    return r++foreign import ccall unsafe "dss"+  c_dss :: SMxV+

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