diff --git a/.gitignore b/.gitignore
index ede5318e4378fb440f0bcaa516296f04aa1f6217..3e2a3aa9ba503d55c6dfe5a8b71957d49457d855 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,7 @@ _CPack_Packages
 install_manifest.txt
 TAGS
 build*/
+wrappers/python/spm/enum.py
 
 #################################################################
 #  https://github.com/github/gitignore/blob/master/C.gitignore
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 31fe917389dfacd7c17682c3ca2d0bb8d707af7b..f25b1b2b92f3821a5f2f534a70e6e2ef2ad79d7d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -34,6 +34,7 @@ test_spm:
       - spm.lcov
       - spm-gcov.log
   script:
+    - source install/bin/spm_env.sh
     - (cd build &&
        eval "ctest
              $TESTS_RESTRICTION
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fdcd443e09423452246d517baaca8961403809ea..1f8fea6b75c26c51890bcddf97b108d7caa9ae68 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -93,7 +93,6 @@ endif()
 
 
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake_modules")
-include(GenSPMPkgConfig)
 include(AddSourceFiles)
 
 # The current version number
@@ -114,12 +113,12 @@ include(CheckSystem)
 
 # SPM depends on Lapacke and CBLAS
 #---------------------------------
-find_package(CBLAS) # Should be REQUIRED for BLAS sequential only
+find_package(CBLAS REQUIRED)
 if(CBLAS_FOUND)
   include_directories(${CBLAS_INCLUDE_DIRS})
 endif()
 
-find_package(LAPACKE) # Should be also REQUIRED
+find_package(LAPACKE REQUIRED)
 if(LAPACKE_FOUND)
   include_directories(${LAPACKE_INCLUDE_DIRS})
 endif()
@@ -246,6 +245,11 @@ install(TARGETS spm
 install(FILES include/spm.h
   DESTINATION include )
 
+### Build pkg-config and environment file
+include(GenSPMPkgConfig)
+spm_generate_pkgconfig_file()
+spm_generate_env_file()
+
 ### Add documented files to the global property
 add_documented_files(
   DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
@@ -267,3 +271,7 @@ add_documented_files(
 
 # Testing executables
 add_subdirectory(tests)
+
+### Wrappers
+add_subdirectory(wrappers)
+
diff --git a/cmake_modules/GenSPMPkgConfig.cmake b/cmake_modules/GenSPMPkgConfig.cmake
index 0a1fc3357f858da57a46796620dc495080a3d9b1..935b432f3b3df4aa1eb9fd34509c4c162e815cfb 100644
--- a/cmake_modules/GenSPMPkgConfig.cmake
+++ b/cmake_modules/GenSPMPkgConfig.cmake
@@ -78,7 +78,7 @@
 # GENERATE_PKGCONFIG_FILE: generate files spm.pc
 #
 ###
-macro(GENERATE_SPM_PKGCONFIG_FILE)
+macro(spm_generate_pkgconfig_file)
 
   set(SPM_PKGCONFIG_LIBS "-lspm")
   set(SPM_PKGCONFIG_LIBS_PRIVATE "-lm")
@@ -89,7 +89,7 @@ macro(GENERATE_SPM_PKGCONFIG_FILE)
 
   set(_output_spm_file "${CMAKE_BINARY_DIR}/spm.pc")
   configure_file(
-    "${CMAKE_CURRENT_SOURCE_DIR}/lib/pkgconfig/spm.pc.in"
+    "${CMAKE_CURRENT_SOURCE_DIR}/tools/spm.pc.in"
     "${_output_spm_file}"
     @ONLY
     )
@@ -98,19 +98,19 @@ macro(GENERATE_SPM_PKGCONFIG_FILE)
     DESTINATION lib/pkgconfig
     )
 
-endmacro(GENERATE_SPM_PKGCONFIG_FILE)
+endmacro(spm_generate_pkgconfig_file)
 
 ###
 #
 # generate_env_file: generate files pastix.pc
 #
 ###
-macro(generate_env_file)
+macro(spm_generate_env_file)
 
     # Create .sh file
     # ---------------
     configure_file(
-      "${CMAKE_CURRENT_SOURCE_DIR}/spm_env.sh.in"
+      "${CMAKE_CURRENT_SOURCE_DIR}/tools/spm_env.sh.in"
       "${CMAKE_BINARY_DIR}/bin/spm_env.sh" @ONLY)
 
     # installation
@@ -118,7 +118,7 @@ macro(generate_env_file)
     install(FILES "${CMAKE_BINARY_DIR}/bin/spm_env.sh"
       DESTINATION bin)
 
-endmacro(generate_env_file)
+endmacro(spm_generate_env_file)
 
 ##
 ## @end file GenPkgConfig.cmake
diff --git a/cmake_modules/morse_cmake b/cmake_modules/morse_cmake
index 1deddb2781f62dbbf0ee9199f569e49f7346397a..8ddf01ff80d9f9fb87465b75c1bfe083bf4e1849 160000
--- a/cmake_modules/morse_cmake
+++ b/cmake_modules/morse_cmake
@@ -1 +1 @@
-Subproject commit 1deddb2781f62dbbf0ee9199f569e49f7346397a
+Subproject commit 8ddf01ff80d9f9fb87465b75c1bfe083bf4e1849
diff --git a/tools/gen_wrappers.py b/tools/gen_wrappers.py
new file mode 100755
index 0000000000000000000000000000000000000000..36fb00b52bdfe978d614042e2396e92ad49ba656
--- /dev/null
+++ b/tools/gen_wrappers.py
@@ -0,0 +1,545 @@
+#!/usr/bin/env python
+"""
+ @file gen_wrappers.py
+
+ Python and Fortran 90 wrapper generator for some of the solverstack
+ libraries, inspired from the PLASMA-OMP fortran generator.
+
+ @copyright 2016-2017 University of Tennessee, US, University of
+                      Manchester, UK. All rights reserved.
+ @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @date 2017-05-04
+
+"""
+import os
+import re
+import argparse
+import wrappers
+
+description = '''\
+Generates Fortran 90 and Python wrappers from the spm header files.'''
+
+help = '''\
+----------------------------------------------------------------------
+Example uses:
+
+   $SPM_SRC_DIR/gen_wrappers.py
+
+----------------------------------------------------------------------
+'''
+
+# ------------------------------------------------------------
+# command line options
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.RawDescriptionHelpFormatter,
+    description=description,
+    epilog=help )
+parser.add_argument('--prefix',        action='store', help='Prefix for variables in Makefile', default='./')
+parser.add_argument('args', nargs='*', action='store', help='Files to process')
+opts = parser.parse_args()
+
+# exclude inline functions from the interface
+exclude_list = [ "inline", "spmIntSort", "spmIntMSort" ]
+
+def polish_file(whole_file):
+    """Preprocessing and cleaning of the header file.
+       Do not change the order of the regular expressions !
+       Works with a long string."""
+
+    clean_file = whole_file
+
+    # borrowed from cfwrapper.py
+    # Remove C comments:
+    clean_file = re.sub(r"(?s)/\*.*?\*/", "", clean_file)
+    clean_file = re.sub(r"//.*", "", clean_file)
+    # Remove C directives (multilines then monoline):
+    clean_file = re.sub(r"(?m)^#(.*[\\][\n])+.*?$", "", clean_file)
+    clean_file = re.sub("(?m)^#.*$", "", clean_file)
+    clean_file = re.sub("(?m)#.*", "", clean_file)
+    # Remove TABs and overnumerous spaces:
+    clean_file = clean_file.replace("\t", " ")
+    clean_file = re.sub("[ ]{2,}", " ", clean_file)
+    # Remove extern C statement:
+    clean_file = re.sub("(?m)^(extern).*$", "", clean_file)
+    # Remove empty lines:
+    clean_file = re.sub(r"(?m)^\n$", "", clean_file)
+    # Merge structs
+    clean_file = re.sub(r"(?m)$", "", clean_file)
+
+    # Merge string into single line
+    clean_file = re.sub(r"\n", "", clean_file)
+
+    # Split the line based on ";" and "}"
+    clean_file = re.sub(r";", "\n", clean_file)
+    clean_file = re.sub(r"}", "}\n", clean_file)
+
+    return clean_file
+
+def preprocess_list(initial_list):
+    """Preprocessing and cleaning of the header file.
+       Works with a list of strings.
+       Produces a new list in which each function, enum or struct
+       corresponds to a single item."""
+
+    # merge braces
+    list1 = []
+    merged_line = ""
+    nopen = 0
+    inStruct = False
+    for line in initial_list:
+
+        if (line.find("struct") > -1):
+            inStruct = True
+
+        if (inStruct):
+            split_character = ","
+        else:
+            split_character = ""
+
+        nopen += line.count("{") - line.count("}")
+        merged_line += line + split_character
+
+        if (nopen <= 0):
+            list1.append(merged_line)
+            merged_line = ""
+            isOpen   = False
+            inStruct = False
+            nopen = 0
+
+    # merge structs
+    list2 = []
+    merged_line = ""
+    for line in list1:
+
+        merged_line += line
+
+        if (line.find("struct") == -1):
+            list2.append(merged_line)
+            merged_line = ""
+
+    # clean orphan braces
+    list3 = []
+    for line in list2:
+
+        if (line.strip() != "}"):
+            list3.append(line)
+
+    #print '\n\n'.join(list3)
+
+    return list3
+
+def parse_triple(string):
+    """Parse string of
+       type (*)name
+       into triple of [type, pointer, name, const]"""
+
+    if "const" in string:
+        const=1
+        string = re.sub(r"const", "", string)
+    else:
+        const=0
+
+    parts = string.split()
+    if (len(parts) < 2 or len(parts) > 3):
+        print("Error: Cannot detect type for ", string)
+
+    type_part = str.strip(parts[0])
+
+    if (len(parts) == 2):
+        name_with_pointer = str.strip(parts[1])
+        if (name_with_pointer.find("**") > -1):
+            pointer_part = "**"
+            name_part = name_with_pointer.replace("**", "")
+        elif (name_with_pointer.find("*") > -1):
+            pointer_part = "*"
+            name_part    = name_with_pointer.replace("*", "")
+        else:
+            pointer_part = ""
+            name_part    = name_with_pointer
+
+    elif (len(parts) == 3):
+        if (str.strip(parts[1]) == "**"):
+            pointer_part = "**"
+            name_part    = str.strip(parts[2])
+        elif (str.strip(parts[1]) == "*"):
+            pointer_part = "*"
+            name_part    = str.strip(parts[2])
+        else:
+            print("Error: Too many parts for ", string)
+
+    name_part = name_part.strip()
+
+    return [type_part, pointer_part, name_part, const]
+
+
+def parse_enums(preprocessed_list):
+    """Each enum will be parsed into a list of its arguments."""
+
+    enum_list = []
+    for proto in preprocessed_list:
+
+        # extract the part of the function from the prototype
+        fun_parts = proto.split("{")
+
+        split_fun = fun_parts[0].strip().split()
+        if len(split_fun) == 0:
+            continue
+
+        if ((split_fun[0] == "enum") or
+            (split_fun[0] == "typedef" and split_fun[1] == "enum")):
+
+
+            if split_fun[0] == "enum":
+                enumname = split_fun[1]
+            else:
+                enumname = split_fun[2]
+            enumname = re.sub(r"_e$", "", enumname)
+            enumname = re.sub(r"^pastix_", "", enumname)
+            enumname = re.sub(r"^spm_", "", enumname)
+
+            args_string = fun_parts[1];
+            args_string = re.sub(r"}", "", args_string)
+            args_list = args_string.split(",")
+            params = [];
+            for args in args_list:
+                args = args.strip();
+                if (args != ""):
+                    values = args.split("=")
+
+                    name = values[0].strip()
+                    if (len(values) > 1):
+                        value = values[1].strip()
+                    else:
+                        if (len(params) > 0):
+                            value = params[len(params)-1][1] + 1
+                        else:
+                            value = 0
+
+                    params.append([name, value])
+
+            enum_list.append([enumname, params])
+
+    return enum_list
+
+
+def parse_structs(preprocessed_list):
+    """Each struct will be parsed into a list of its arguments."""
+
+    struct_list = []
+    for proto in preprocessed_list:
+
+        # extract the part of the function from the prototype
+        fun_parts = proto.split("{")
+
+        if (fun_parts[0].find("struct") > -1):
+            args_string = fun_parts[1]
+            parts = args_string.split("}")
+            args_string = parts[0].strip()
+            args_string = re.sub(r"volatile", "", args_string)
+            if (len(parts) > 1):
+                name_string = parts[1]
+                name_string = re.sub(r"(?m),", "", name_string)
+                name_string = name_string.strip()
+            else:
+                print("Error: Cannot detect name for ", proto)
+                name_string = "name_not_recognized"
+
+            args_list = args_string.split(",")
+            params = [];
+            params.append(["struct","",name_string])
+            for arg in args_list:
+                if (not (arg == "" or arg == " ")):
+                    params.append(parse_triple(arg))
+
+            struct_list.append(params)
+            wrappers.derived_types.append(name_string)
+
+    # reorder the list so that only defined types are exported
+    goAgain = True
+    while (goAgain):
+        goAgain = False
+        for istruct in range(0,len(struct_list)-1):
+            struct = struct_list[istruct]
+            for j in range(1,len(struct)-1):
+                type_name = struct_list[istruct][j][0]
+
+                if (type_name in wrappers.derived_types):
+
+                    # try to find the name in the registered types
+                    definedEarlier = False
+                    for jstruct in range(0,istruct):
+                        struct2 = struct_list[jstruct]
+                        that_name = struct2[0][2]
+                        if (that_name == type_name):
+                            definedEarlier = True
+
+                    # if not found, try to find it behind
+                    if (not definedEarlier):
+                        definedLater = False
+                        for jstruct in range(istruct+1,len(struct_list)-1):
+                            struct2 = struct_list[jstruct]
+                            that_name = struct2[0][2]
+                            if (that_name == type_name):
+                                index = jstruct
+                                definedLater = True
+
+                        # swap the entries
+                        if (definedLater):
+                            print("Swapping " + struct_list[istruct][0][2] + " and " + struct_list[index][0][2])
+                            tmp = struct_list[index]
+                            struct_list[index] = struct_list[istruct]
+                            struct_list[istruct] = tmp
+                            goAgain = True
+                        else:
+                            print("Error: Cannot find a derived type " + type_name + " in imported structs.")
+
+    return struct_list
+
+
+def parse_prototypes(preprocessed_list):
+    """Each prototype will be parsed into a list of its arguments."""
+
+    function_list = []
+    for proto in preprocessed_list:
+
+        if (proto.find("(") == -1):
+            continue
+
+        # extract the part of the function from the prototype
+        fun_parts = proto.split("(")
+        fun_def  = str.strip(fun_parts[0])
+
+        exclude_this_function = False
+        for exclude in exclude_list:
+            if (fun_def.find(exclude) != -1):
+                exclude_this_function = True
+
+        if (exclude_this_function):
+            continue
+
+        # clean keywords
+        fun_def = fun_def.replace("static", "")
+
+        # extract arguments from the prototype and make a list from them
+        if (len(fun_parts) > 1):
+            fun_args = fun_parts[1]
+        else:
+            fun_args = ""
+
+        fun_args = fun_args.split(")")[0]
+        fun_args = fun_args.replace(";", "")
+        fun_args = re.sub(r"volatile", "", fun_args)
+        fun_args = fun_args.replace("\n", "")
+        fun_args_list = fun_args.split(",")
+
+        # generate argument list
+        argument_list = []
+        # function itself on the first position
+        argument_list.append(parse_triple(fun_def))
+        # append arguments
+        for arg in fun_args_list:
+            if (not (arg == "" or arg == " ")):
+                argument_list.append(parse_triple(arg))
+
+        # add it only if there is no duplicity with previous one
+        is_function_already_present = False
+        fun_name = argument_list[0][2]
+        for fun in function_list:
+            if (fun_name == fun[0][2]):
+                is_function_already_present = True
+
+        if (not is_function_already_present):
+            function_list.append(argument_list)
+
+    return function_list
+
+def write_module(f, wrapper, generator, enum_list, struct_list, function_list):
+    """Generate a single Fortran module. Its structure will be:
+       enums converted to constants
+       structs converted to derived types
+       interfaces of all C functions
+       Fortran wrappers of the C functions"""
+
+    modulefile = open( f['filename'], "w" )
+
+    f_interface = generator.header( f )
+    modulefile.write(f_interface + "\n")
+
+    # enums
+    if (len(enum_list) > 0):
+        for enum in enum_list:
+            f_interface = generator.enum( f, enum )
+            modulefile.write(f_interface + "\n")
+
+    # derived types
+    if (len(struct_list) > 0):
+        for struct in struct_list:
+            f_interface = generator.struct(struct)
+            modulefile.write(f_interface + "\n")
+
+    # functions
+    if (len(function_list) > 0):
+        for function in function_list:
+            f_interface = generator.function(function)
+            modulefile.write(f_interface + "\n")
+
+        if wrapper:
+            modulefile.write("contains\n\n")
+            modulefile.write("  " + "! Wrappers of the C functions.\n")
+
+            for function in function_list:
+                f_wrapper = generator.wrapper(function)
+                modulefile.write(f_wrapper + "\n")
+
+    footer=generator.footer( f )
+    modulefile.write(footer + "\n")
+
+    modulefile.close()
+
+    return f['filename']
+
+
+enums_python_coeftype='''
+    @staticmethod
+    def getptype ( dtype ):
+        np_dict = {
+            np.dtype('float32')    : coeftype.Float,
+            np.dtype('float64')    : coeftype.Double,
+            np.dtype('complex64')  : coeftype.Complex32,
+            np.dtype('complex128') : coeftype.Complex64,
+        }
+        if dtype in np_dict:
+            return np_dict[dtype]
+        else:
+            return -1
+
+    @staticmethod
+    def getctype ( flttype ):
+        class c_float_complex(Structure):
+            _fields_ = [("r",c_float),("i", c_float)]
+        class c_double_complex(Structure):
+            _fields_ = [("r",c_double),("i", c_double)]
+
+        np_dict = {
+            coeftype.Float     : c_float,
+            coeftype.Double    : c_double,
+            coeftype.Complex32 : c_float_complex,
+            coeftype.Complex64 : c_double_complex
+        }
+        if flttype in np_dict:
+            return np_dict[flttype]
+        else:
+            return -1
+
+    @staticmethod
+    def getnptype ( flttype ):
+        np_dict = {
+            coeftype.Float     : np.dtype('float32'),
+            coeftype.Double    : np.dtype('float64'),
+            coeftype.Complex32 : np.dtype('complex64'),
+            coeftype.Complex64 : np.dtype('complex128')
+        }
+        if flttype in np_dict:
+            return np_dict[flttype]
+        else:
+            return -1
+'''
+
+enums_fortran_footer='''
+  integer, parameter :: spm_int_t = SPM_INT_KIND
+
+contains
+
+  function spm_getintsize()
+    integer :: spm_getintsize
+    spm_getintsize = SPM_INT_KIND
+    return
+  end function spm_getintsize
+'''
+
+spm_enums = {
+    'filename' : [ "include/spm_const.h" ],
+    'python'   : { 'filename'    : "wrappers/python/spm/enum.py.in",
+                   'description' : "SPM python wrapper to define enums and datatypes",
+                   'header'      : "# Start with __ to prevent broadcast to file importing enum\n__spm_int__ = @SPM_PYTHON_INTEGER@\n",
+                   'footer'      : "",
+                   'enums'       : { 'coeftype' : enums_python_coeftype,
+                                     'mtxtype'  : "    SymPosDef = trans.ConjTrans + 1\n    HerPosDef = trans.ConjTrans + 2\n" }
+    },
+    'fortran'  : { 'filename'    : "wrappers/fortran90/src/spm_enums.F90",
+                   'description' : "SPM fortran 90 wrapper to define enums and datatypes",
+                   'header'      : "  implicit none\n",
+                   'footer'      : enums_fortran_footer,
+                   'enums'       : { 'mtxtype'  : "    enumerator :: SpmSymPosDef = SpmConjTrans + 1\n    enumerator :: HerPosDef    = SpmConjTrans + 2\n" }
+    },
+}
+
+spm = {
+    'filename' : [ "include/spm.h" ],
+    'python'   : { 'filename'    : "wrappers/python/spm/__spm__.py",
+                   'description' : "SPM python wrapper",
+                   'header'      : "from . import libspm\nfrom .enum import __spm_int__\n",
+                   'footer'      : "",
+                   'enums'       : {}
+    },
+    'fortran'  : { 'filename'    : "wrappers/fortran90/src/spmf.f90",
+                   'description' : "SPM Fortran 90 wrapper",
+                   'header'      : "  use spm_enums\n  implicit none\n",
+                   'footer'      : "",
+                   'enums'       : {}
+    },
+}
+
+def main():
+
+    # common cleaned header files
+    preprocessed_list = []
+
+    # source header files
+    for f in [ spm_enums, spm ]:
+        preprocessed_list = []
+        for filename in f['filename']:
+
+            # source a header file
+            c_header_file = open(filename, 'r').read()
+
+            # clean the string (remove comments, macros, etc.)
+            clean_file = polish_file(c_header_file)
+
+            # convert the string to a list of strings
+            initial_list = clean_file.split("\n")
+
+            # process the list so that each enum, struct or function
+            # are just one item
+            nice_list = preprocess_list(initial_list)
+
+            # compose all files into one big list
+            preprocessed_list += nice_list
+
+            # register all enums
+            enum_list = parse_enums(preprocessed_list)
+
+            # register all structs
+            struct_list = parse_structs(preprocessed_list)
+
+            # register all individual functions and their signatures
+            function_list = parse_prototypes(preprocessed_list)
+
+        # export the module
+        modulefilename = write_module( f['fortran'], True,
+                                       wrappers.wrap_fortran,
+                                       enum_list, struct_list, function_list)
+        print( "Exported file: " + modulefilename )
+
+        modulefilename = write_module( f['python'], False,
+                                       wrappers.wrap_python,
+                                       enum_list, struct_list, function_list)
+        print( "Exported file: " + modulefilename )
+
+# execute the program
+main()
diff --git a/lib/pkgconfig/spm.pc.in b/tools/spm.pc.in
similarity index 100%
rename from lib/pkgconfig/spm.pc.in
rename to tools/spm.pc.in
diff --git a/tools/spm_env.sh.in b/tools/spm_env.sh.in
new file mode 100644
index 0000000000000000000000000000000000000000..1ae1a3813e0cf805866a2bc7a4db061d4661174b
--- /dev/null
+++ b/tools/spm_env.sh.in
@@ -0,0 +1,36 @@
+#
+#  @file spm_env.sh
+#
+#  @copyright 2016-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                       Univ. Bordeaux. All rights reserved.
+#
+#  @version 6.0.0
+#  @author Mathieu Faverge
+#  @date 2017-06-24
+#
+#!/bin/sh
+
+LIB=spm
+
+export SPM_DIR=@CMAKE_INSTALL_PREFIX@
+
+for i in PATH DYLD_LIBRARY_PATH LD_LIBRARY_PATH LIBRARY_PATH LD_RUN_PATH INCLUDE INCLUDE_PATH PKG_CONFIG_PATH PYTHONPATH
+do
+
+  for j in /spm
+  do
+    cmd1="echo \\\"\$$i\\\" | sed -E 's+^(\(.*:|\))[^:]*${j}[^:]*(\(|:.*\))$+\1\2+' | sed 's/::/:/' | sed 's/^://' | sed 's/:$//' "
+    temp=`eval $cmd1`;
+    eval "$i=$temp";
+  done
+done
+
+export PATH=$PATH:$SPM_DIR/bin
+export LD_RUN_PATH=$LD_RUN_PATH:$SPM_DIR/lib
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$SPM_DIR/lib
+export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:$SPM_DIR/lib
+export LIBRARY_PATH=$LIBRARY_PATH:$SPM_DIR/lib
+export PYTHONPATH=$PYTHONPATH:$SPM_DIR/lib/python
+export INCLUDE=$INCLUDE:$SPM_DIR/include
+export INCLUDE_PATH=$INCLUDE_PATH:$SPM_DIR/include
+export PKG_CONFIG_PATH=$PKG_CONFIG_PATH:$SPM_DIR/lib/pkgconfig
diff --git a/tools/wrappers/__init__.py b/tools/wrappers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c910afd6ba36daa575037c3bf59db421fcd61dcf
--- /dev/null
+++ b/tools/wrappers/__init__.py
@@ -0,0 +1,41 @@
+"""
+Wrappers
+========
+
+ @file wrappers/__init__.py
+
+ PaStiX wrapper generators module intialization
+
+ @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Mathieu Faverge
+ @date 2017-05-04
+
+"""
+
+# translation_table with names of auxiliary variables
+return_variables_dict = {
+    "double":            ("value"),
+    "float":             ("value"),
+    "pastix_int_t":      ("value"),
+    "pastix_spm_t":      ("spmo"),
+    "pastix_order_t":    ("order"),
+    "spm_int_t"   :      ("value"),
+    "spmatrix_t":        ("spmo"),
+}
+
+# global list used to determine derived types
+derived_types = [ 'spmatrix_t', 'spm_int_t', 'pastix_int_t', 'pastix_data_t', 'pastix_order_t']
+
+# name arrays which will be translated to assumed-size arrays, e.g. pA(*)
+arrays_names_2D = ["pA", "pB", "pC", "pAB", "pQ", "pX", "pAs"]
+arrays_names_1D = ["colptr", "rowptr", "loc2glob", "dofs", "row",
+                   "iparm", "dparm", "bindtab", "perm", "invp", "schur_list",
+                   "rang", "tree" ]
+
+__all__ = [ 'return_variables_dict', 'derived_types', 'arrays_names_1D', 'arrays_names_2D' ]
+
+from .wrap_python  import *
+from .wrap_fortran import *
diff --git a/tools/wrappers/wrap_fortran.py b/tools/wrappers/wrap_fortran.py
new file mode 100644
index 0000000000000000000000000000000000000000..e4250f832399f6fab099d1d1086cc0807713a482
--- /dev/null
+++ b/tools/wrappers/wrap_fortran.py
@@ -0,0 +1,474 @@
+#!/usr/bin/env python
+"""
+Wrapper Fortran 90
+==================
+
+ @file wrappers/wrap_fortran.py
+
+ PaStiX generator for the Fortran 90 wrapper
+
+ @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Mathieu Faverge
+ @date 2017-05-04
+
+"""
+import os
+import re
+import argparse
+from . import *
+
+# set indentation in the f90 file
+tab = "  "
+indent = "   "
+
+itab=2
+iindent=3
+
+# translation_table of types
+types_dict = {
+    "int":            ("integer(kind=c_int)"),
+    "spm_coeftype_t": ("integer(c_int)"),
+    "spm_dir_t":      ("integer(c_int)"),
+    "spm_trans_t":    ("integer(c_int)"),
+    "spm_uplo_t":     ("integer(c_int)"),
+    "spm_diag_t":     ("integer(c_int)"),
+    "spm_side_t":     ("integer(c_int)"),
+    "spm_driver_t":   ("integer(c_int)"),
+    "spm_fmttype_t":  ("integer(c_int)"),
+    "spm_layout_t":   ("integer(c_int)"),
+    "spm_normtype_t": ("integer(c_int)"),
+    "spm_rhstype_t":  ("integer(c_int)"),
+    "spm_mtxtype_t":  ("integer(c_int)"),
+    "spmatrix_t":     ("type(spmatrix_t)"),
+    "spm_int_t":      ("integer(kind=spm_int_t)"),
+    "size_t":         ("integer(kind=c_size_t)"),
+    "char":           ("character(kind=c_char)"),
+    "double":         ("real(kind=c_double)"),
+    "float":          ("real(kind=c_float)"),
+    "spm_complex64_t":("complex(kind=c_double_complex)"),
+    "spm_complex32_t":("complex(kind=c_float_complex)"),
+    "void":           ("type(c_ptr)"),
+    "MPI_Comm":       ("integer(kind=c_int)"),
+    "FILE":           ("type(c_ptr)"),
+}
+
+def iso_c_interface_type(arg, return_value, list):
+    """Generate a declaration for a variable in the interface."""
+
+    if (arg[1] == "*" or arg[1] == "**"):
+        is_pointer = True
+    else:
+        is_pointer = False
+
+    if (is_pointer):
+        f_type = "type(c_ptr)"
+    else:
+        f_type = types_dict[arg[0]]
+
+    if (not return_value and arg[1] != "**"):
+        f_type += ", "
+        f_pointer = "value"
+    else:
+        f_pointer = ""
+
+    f_name = arg[2]
+
+    list.append( [ f_type, f_pointer, f_name ] );
+    return len(f_type + f_pointer)
+
+def iso_c_wrapper_type(arg, args_list, args_size):
+    """Generate a declaration for a variable in the Fortran wrapper."""
+
+    if (arg[1] == "*" or arg[1] == "**"):
+        is_pointer = True
+    else:
+        is_pointer = False
+
+    if (is_pointer and arg[0].strip() == "void"):
+        f_type = "type(c_ptr), "
+    else:
+        f_type = types_dict[arg[0]] + ", "
+
+    if is_pointer and not arg[3]:
+        f_intent = "intent(inout)"
+    else:
+        f_intent = "intent(in)"
+
+    if (is_pointer):
+        f_intent += ", "
+        if (arg[1] == "*"):
+           f_target = "target"
+        else:
+           f_target = "pointer"
+    else:
+        f_target = ""
+
+    f_name = arg[2]
+
+    # Force case of myorder
+    if f_name == "myorder":
+        f_intent = "intent(in), "
+        f_target = "pointer"
+
+    # detect array argument
+    if   (is_pointer and f_name in arrays_names_2D):
+        f_array = "(:,:)"
+    elif (is_pointer and f_name in arrays_names_1D):
+        f_array = "(:)"
+    else:
+        f_array = ""
+
+    f_line = f_type + f_intent + f_target + " :: " + f_name + f_array
+
+    args_size[0] = max(args_size[0], len(f_type))
+    args_size[1] = max(args_size[1], len(f_intent))
+    args_size[2] = max(args_size[2], len(f_target))
+    args_list.append( (f_type, f_intent, f_target, f_name, f_array ) )
+    return f_line
+
+class wrap_fortran:
+
+    @staticmethod
+    def header( f ):
+        filename = os.path.basename( f['filename'] )
+        modname = re.sub(r".f90", "", filename, flags=re.IGNORECASE)
+        header = '''
+!
+! @file '''+ filename +'''
+!
+! ''' + f['description'] + '''
+!
+! @copyright 2017      Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+!                      Univ. Bordeaux. All rights reserved.
+!
+! @version 6.0.0
+! @author Mathieu Faverge
+! @date 2017-01-01
+!
+! This file has been automatically generated with gen_wrappers.py
+!
+module ''' + modname + '''
+  use iso_c_binding
+'''
+
+        if f['header'] != "":
+            header += f['header']
+
+        return header
+
+    @staticmethod
+    def footer( f ):
+        filename = os.path.basename( f['filename'] )
+        modname = re.sub(r".f90", "", filename, flags=re.IGNORECASE)
+        footer = f['footer'] + '''
+end module ''' + modname
+        return footer
+
+    @staticmethod
+    def enum( f, enum ):
+        """Generate an interface for an enum.
+           Translate it into constants."""
+
+        ename  = enum[0]
+        params = enum[1]
+
+        # initialize a string with the fortran interface
+        f_interface  = tab + "! enum " + ename + "\n"
+        f_interface += tab + "enum, bind(C)\n"
+
+        # loop over the arguments of the enum to get max param length
+        length=0
+        for param in params:
+            length= max( length, len(param[0]) )
+        fmt="%-"+ str(length) + "s"
+
+        # Increment for index array enums
+        inc = 0
+        if ename[1:5] == "parm":
+            inc=1
+
+        # loop over the arguments of the enum
+        for param in params:
+            name  = param[0]
+            if isinstance(param[1],int):
+                if name[1:10] == "PARM_SIZE":
+                    value = str(param[1])
+                else:
+                    value = str(param[1] + inc)
+            else:
+                value = str(param[1])
+            f_interface += tab + "   enumerator :: " + format(fmt % name) + " = " + value + "\n"
+
+        f_interface += tab + "end enum\n"
+        return f_interface
+
+    @staticmethod
+    def struct(struct):
+        """Generate an interface for a struct.
+           Translate it into a derived type."""
+
+        # initialize a string with the fortran interface
+        f_interface = ""
+
+        s = itab
+        name = struct[0][2]
+        f_interface += s*" " + "type, bind(c) :: " + name + "\n"
+
+        # loop over the arguments of the struct to get the length
+        s += iindent
+        slist = []
+        length= 0
+        for j in range(1,len(struct)):
+            length = max( length, iso_c_interface_type(struct[j], True, slist) )
+        fmt = s*" " + "%-"+ str(length) + "s :: %s\n"
+
+        # loop over the arguments of the struct
+        for j in range(0,len(slist)):
+            f_interface += format( fmt % (slist[j][0], slist[j][2]) )
+
+        s -= iindent
+        f_interface += s*" " + "end type " + name + "\n"
+
+        return f_interface
+
+    @staticmethod
+    def function(function):
+        """Generate an interface for a function."""
+
+        # is it a function or a subroutine
+        if (function[0][0] == "void"):
+            is_function = False
+            symbol="subroutine"
+        else:
+            is_function = True
+            symbol="function"
+
+        c_symbol = function[0][2]
+        f_symbol = function[0][2] + "_c"
+
+        used_derived_types = set([])
+        for arg in function:
+            type_name = arg[0]
+            if (type_name in derived_types):
+                used_derived_types.add(type_name)
+
+        # initialize a string with the fortran interface
+        s = itab
+        f_interface = s*" " + "interface\n"
+        s += iindent
+
+        f_interface += s*" " + symbol + " " + f_symbol + "("
+        s += itab
+
+        initial_len = s + len(symbol + " " + f_symbol + "(" )
+
+        # loop over the arguments to compose the first line
+        s += iindent
+        for j in range(1,len(function)):
+            if (j != 1):
+                f_interface += ", "
+                initial_len += 2
+
+            l = len(function[j][2])
+            if ((initial_len + l) > 77):
+                f_interface += "&\n" + s*" "
+                initial_len = s
+
+            f_interface += function[j][2]
+            initial_len += l
+
+        f_interface += ") &\n"
+        f_interface += s*" " + "bind(c, name='" + c_symbol +"')\n"
+        s -= iindent
+
+        # add common header
+        f_interface += s*" " + "use iso_c_binding\n"
+        # import derived types
+        for derived_type in used_derived_types:
+            f_interface += s*" " + "import " + derived_type +"\n"
+        f_interface += s*" " + "implicit none\n"
+
+        plist = []
+        length = 0
+        # add the return value of the function
+        if (is_function):
+            l = iso_c_interface_type(function[0], True, plist ) + 2
+            length = max( length, l );
+            plist[0][2] += "_c"
+
+        # loop over the arguments to describe them
+        for j in range(1,len(function)):
+            length = max( length, iso_c_interface_type(function[j], False, plist ) )
+
+        # loop over the parameters
+        for j in range(0,len(plist)):
+            fmt = s*" " + "%s%" + str(length-len(plist[j][0])) + "s :: %s\n"
+            f_interface += format( fmt % (plist[j][0], plist[j][1], plist[j][2]) )
+
+        s -= itab
+        f_interface += s*" " + "end " + symbol + " " + f_symbol + "\n"
+
+        s -= iindent
+        f_interface += s*" " + "end interface\n"
+
+        return f_interface
+
+    @staticmethod
+    def wrapper(function):
+        """Generate a wrapper for a function.
+           void functions in C will be called as subroutines,
+           functions in C will be turned to subroutines by appending
+           the return value as the last argument."""
+
+        return_type = function[0][0]
+        return_pointer = function[0][1]
+        return_var = ""
+
+        # is it a function or a subroutine
+        if (return_type == "void"):
+            is_function = False
+            f_return = "call "
+        else:
+            is_function = True
+            if (return_type == "int"):
+                return_var = "info"
+            else:
+                return_var = return_variables_dict[return_type]
+            f_return = return_var + " = "
+
+        c_symbol = function[0][2]
+        f_symbol = c_symbol + "_c"
+
+        # loop over the arguments to compose the first line and call line
+        s = itab
+        signature_line = s*" " + "subroutine " + c_symbol + "("
+        call_line      = ""
+        signature_line_length = len(signature_line)
+        call_line_length      = s + itab + len(f_return + f_symbol + "(" )
+        double_pointers = []
+        args_list = []
+        args_size = [ 0, 0, 0 ]
+        length = 0
+        for j in range(1,len(function)):
+            # pointers
+            arg_type    = function[j][0]
+            arg_pointer = function[j][1]
+            arg_name    = function[j][2]
+
+            isfile = (arg_type == "FILE")
+
+            if (j != 1):
+                if not isfile:
+                    signature_line += ", "
+                    signature_line_length += 2
+                call_line += ", "
+                call_line_length += 2
+
+            # Signature line (do not make the argument list too long)
+            if not isfile:
+                line = iso_c_wrapper_type(function[j], args_list, args_size)
+
+                l = len(arg_name)
+                if ((signature_line_length + l) > 77):
+                    signature_line_length = s+itab+iindent
+                    signature_line += "&\n" + signature_line_length*" "
+                signature_line += arg_name
+                signature_line_length += l
+
+            # Call line
+            if isfile:
+                call_param = "c_null_ptr"
+            elif (arg_pointer == "**"):
+                aux_name = arg_name + "_aux"
+                call_param = aux_name
+                double_pointers.append(arg_name)
+            elif (arg_pointer == "*") and (not "type(c_ptr), " in args_list[len(args_list)-1]):
+                call_param = "c_loc(" + arg_name + ")"
+            else:
+                call_param = arg_name
+
+            # Do not make the call line too long
+            l = len(call_param)
+            if ((call_line_length + l) > 77):
+                call_line_length = s+itab+iindent+itab
+                call_line += "&\n" + call_line_length*" "
+            call_line += call_param
+            call_line_length += l
+
+
+        # initialize a string with the fortran interface
+        f_wrapper = signature_line
+        if (is_function):
+            if (len(function) > 1):
+                f_wrapper += ", "
+                signature_line_length += 2
+
+            return_type = function[0][0]
+            return_pointer = function[0][1]
+
+            l = len(return_var)
+            if ((signature_line_length + l) > 78):
+                signature_line_length = s+itab+iindent
+                f_wrapper += "&\n" + signature_line_length*" "
+
+            f_wrapper += return_var
+
+        f_wrapper += ")\n"
+
+        s += itab
+        # add common header
+        f_wrapper += s*" " + "use iso_c_binding\n"
+        f_wrapper += s*" " + "implicit none\n"
+
+        # add the return info value of the function
+        if (is_function):
+            f_type = types_dict[return_type] + ", "
+            f_intent = "intent(out)"
+            if (function[0][1] == "*"):
+                f_intent += ", "
+                f_target = "pointer"
+            else:
+                f_target = ""
+
+            args_size[0] = max(args_size[0], len(f_type))
+            args_size[1] = max(args_size[1], len(f_intent))
+            args_size[2] = max(args_size[2], len(f_target))
+            args_list.append( (f_type, f_intent, f_target, return_var, "" ) )
+
+        # loop over potential double pointers and generate auxiliary variables for them
+        init_line = ""
+        for double_pointer in double_pointers:
+            aux_name = double_pointer + "_aux"
+
+            args_size[0] = max(args_size[0], len("type(c_ptr)"))
+            args_list.append( ("type(c_ptr)", "", "", aux_name, "") )
+            init_line += s*" " + aux_name + " = c_loc(" + double_pointer + ")\n"
+
+        # loop over the arguments to describe them
+        fmt = s*" " + "%-" + str(args_size[0]) + "s%-" + str(args_size[1]) + "s%-" + str(args_size[2]) + "s :: %s%s\n"
+        for j in range(0,len(args_list)):
+            f_wrapper += format( fmt % args_list[j] )
+
+        f_wrapper += "\n"
+        f_wrapper += init_line
+        if init_line != "":
+            f_wrapper += "\n"
+
+        # generate the call to the C function
+        if (is_function and return_pointer == "*"):
+            f_wrapper += s*" " + "call c_f_pointer(" + f_symbol + "(" + call_line + "), " + return_var + ")\n"
+        else:
+            f_wrapper += s*" " + f_return + f_symbol + "(" + call_line + ")\n"
+
+        # loop over potential double pointers and translate them to Fortran pointers
+        for double_pointer in double_pointers:
+            aux_name = double_pointer + "_aux"
+            f_wrapper += s*" " + "call c_f_pointer(" + aux_name + ", " + double_pointer + ")\n"
+
+        s -= itab
+        f_wrapper += s*" " + "end subroutine " + c_symbol + "\n"
+
+        return f_wrapper
diff --git a/tools/wrappers/wrap_python.py b/tools/wrappers/wrap_python.py
new file mode 100644
index 0000000000000000000000000000000000000000..60945858ffae065efb42db29682bb9d4d9952044
--- /dev/null
+++ b/tools/wrappers/wrap_python.py
@@ -0,0 +1,521 @@
+#!/usr/bin/env python
+"""
+Wrapper Python
+==============
+
+ @file wrappers/wrap_python.py
+
+ PaStiX generator for the python wrapper
+
+ @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Mathieu Faverge
+ @date 2017-05-04
+
+"""
+import os
+import re
+import argparse
+from . import *
+
+indent="    "
+iindent=4
+
+# translation_table of types
+types_dict = {
+    "int":            ("c_int"),
+    "spm_coeftype_t": ("c_int"),
+    "spm_dir_t":      ("c_int"),
+    "spm_trans_t":    ("c_int"),
+    "spm_uplo_t":     ("c_int"),
+    "spm_diag_t":     ("c_int"),
+    "spm_side_t":     ("c_int"),
+    "spm_driver_t":   ("c_int"),
+    "spm_fmttype_t":  ("c_int"),
+    "spm_layout_t":   ("c_int"),
+    "spm_normtype_t": ("c_int"),
+    "spm_rhstype_t":  ("c_int"),
+    "spm_mtxtype_t":  ("c_int"),
+    "spm_int_t":      ("__spm_int__"),
+    "spmatrix_t":     ("pyspm_spmatrix_t"),
+    "size_t":         ("c_size_t"),
+    "char":           ("c_char"),
+    "double":         ("c_double"),
+    "float":          ("c_float"),
+    "void":           ("c_void"),
+    "MPI_Comm":       ("c_int"),
+    "FILE":           ("c_void"),
+}
+
+def iso_c_interface_type(arg, return_value, args_list, args_size):
+    """Generate a declaration for a variable in the interface."""
+
+    if (arg[1] == "*" or arg[1] == "**"):
+        is_pointer = True
+    else:
+        is_pointer = False
+
+    f_type = types_dict[arg[0]]
+    if is_pointer:
+        if f_type == "c_void":
+            f_type = "c_void_p"
+        elif f_type == "c_char":
+            f_type = "c_char_p"
+        elif f_type == "c_int":
+            f_type = "c_int_p"
+        else:
+            f_type = "POINTER("+f_type+")"
+
+    if (not return_value and arg[1] != "**"):
+        f_pointer = "value"
+    else:
+        f_pointer = ""
+
+    f_name = format("\"%s\", " % arg[2] )
+
+    args_size[0] = max(args_size[0], len(f_name))
+    args_size[1] = max(args_size[1], len(f_type))
+    args_list.append( [ f_name, f_type ] );
+
+def iso_c_wrapper_type(arg, args_list, args_size):
+    """Generate a declaration for a variable in the Fortran wrapper."""
+
+    if (arg[1] == "*" or arg[1] == "**"):
+        is_pointer = True
+    else:
+        is_pointer = False
+
+    isfile = (arg[0] == "FILE")
+    f_type = types_dict[arg[0]]
+
+    if is_pointer:
+        if f_type == "c_void":
+            f_type = "c_void_p"
+        elif f_type == "c_char":
+            f_type = "c_char_p"
+        elif f_type == "c_int":
+            f_type = "c_int_p"
+        else:
+            f_type = "POINTER("+f_type+")"
+
+    f_name = arg[2]
+    f_call = f_name
+
+    if isfile:
+        f_name = ""
+        f_call = "None"
+
+    # detect array argument
+    if (is_pointer and f_name in arrays_names_2D):
+        f_call = f_name + ".ctypes.data_as( " + f_type + " )"
+    elif (is_pointer and f_name in arrays_names_1D):
+        f_call = f_name + ".ctypes.data_as( " + f_type + " )"
+
+    if arg[1] == "**":
+        f_call = "pointer( " + f_call + " )"
+
+    args_size[0] = max(args_size[0], len(f_name))
+    args_size[1] = max(args_size[1], len(f_type))
+    args_size[2] = max(args_size[2], len(f_call))
+    args_list.append( [f_name, f_type, f_call ] )
+
+class wrap_python:
+
+    @staticmethod
+    def header( f ):
+        filename = os.path.basename( f['filename'] )
+        filename = re.sub(r"\.in", "", filename)
+        header = '''"""
+
+ @file ''' + filename + '''
+
+ ''' + f['description'] + '''
+
+ @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @author Louis Poirel
+ @date 2017-05-04
+
+This file has been automatically generated with gen_wrappers.py
+
+"""
+from ctypes import *
+import numpy as np
+'''
+        if f['header'] != "":
+            header += "\n" + f['header']
+        return header;
+
+    @staticmethod
+    def footer( f ):
+        filename = os.path.basename( f['filename'] )
+        modname = re.sub(r".f90", "", filename, flags=re.IGNORECASE)
+        footer = f['footer']
+        return footer
+
+    @staticmethod
+    def enum( f, enum ):
+        """Generate an interface for an enum.
+           Translate it into constants."""
+
+        ename  = enum[0]
+        params = enum[1]
+
+        # initialize a string with the fortran interface
+        py_interface = "class " + ename + ":\n"
+
+        # loop over the arguments of the enum to get max param length
+        length=0
+        for param in params:
+            # Convert IPARM/DPARM to lower case
+            if param[0][1:5] == "PARM":
+                param[0] = re.sub(r"[ID]PARM_", "", param[0])
+                param[0] = param[0].lower()
+
+            # Remove Pastix from everything
+            param[0] = re.sub(r"Pastix", "", param[0])
+            param[0] = re.sub(r"Spm", "", param[0])
+
+            if ename == "error":
+                param[0] = re.sub(r"PASTIX_", "", param[0])
+                param[0] = re.sub(r"SPM_", "", param[0])
+                param[0] = re.sub(r"ERR_", "", param[0])
+            elif ename == "fact_mode" or ename == "factotype" or ename == "solv_mode":
+                param[0] = re.sub(r"Fact", "", param[0])
+                param[0] = re.sub(r"Solv", "", param[0])
+                param[0] = re.sub(r"Mode", "", param[0])
+            elif ename == "scheduler":
+                param[0] = re.sub(r"Sched", "", param[0])
+            elif ename == "threadmode":
+                param[0] = re.sub(r"Thread", "", param[0])
+            elif ename[0:8] == "compress":
+                param[0] = re.sub(r"Compress", "", param[0])
+                param[0] = re.sub(r"When", "", param[0])
+                param[0] = re.sub(r"Method", "", param[0])
+            elif ename == "rhstype":
+                param[0] = re.sub(r"Rhs", "", param[0])
+            elif ename == "trans":
+                param[0] = param[0]
+            elif ename == "mtxtype":
+                param[1] = re.sub(r"Pastix", "trans.", param[1])
+                param[1] = re.sub(r"Spm", "trans.", param[1])
+            elif ename == "normtype":
+                param[0] = re.sub(r"Norm", "", param[0])
+            else:
+                param[0] = re.sub(ename, "", param[0], flags=re.IGNORECASE)
+            length = max( length, len(param[0]))
+        fmt="%-"+ str(length) + "s"
+
+        # loop over the arguments of the enum
+        for param in params:
+            name  = param[0]
+            value = str(param[1])
+
+            py_interface += indent + format(fmt % name) + " = " + value + "\n"
+
+        if ename in f['enums']:
+            py_interface += f['enums'][ename]
+
+        return py_interface
+
+    @staticmethod
+    def struct(struct):
+        """Generate an interface for a struct.
+           Translate it into a derived type."""
+
+        # initialize a string with the fortran interface
+        py_interface = ""
+
+        s = 0
+        name = struct[0][2]
+        name = re.sub(r"pastix_", "", name)
+        py_interface +=  "class pyspm_" + name + "(Structure):\n"
+
+        s = iindent
+        py_interface +=  s*" " + "_fields_ = ["
+        headline = (s+len("_fields_ = ["))*" "
+
+        slist = []
+        ssize = [ 0, 0 ]
+
+        # loop over the arguments of the enum
+        for j in range(1,len(struct)):
+            iso_c_interface_type(struct[j], True, slist, ssize)
+
+        s += iindent
+        fmt = "(%-"+ str(ssize[0]) + "s %-"+ str(ssize[1]) +"s)"
+
+        for j in range(0,len(slist)):
+            if (j > 0):
+                py_interface += ",\n" + headline
+
+            py_interface += format( fmt % (slist[j][0], slist[j][1]) )
+
+        py_interface += " ]\n"
+
+        return py_interface
+
+    @staticmethod
+    def function(function):
+        """Generate an interface for a function."""
+
+        return_type    = function[0][0]
+        return_pointer = function[0][1]
+
+        # is it a function or a subroutine
+        if (return_type == "void"):
+            is_function = False
+        else:
+            is_function = True
+
+        c_symbol = function[0][2]
+        if "pastix" in c_symbol:
+            libname = "libpastix"
+            prefix  = "pypastix_"
+        elif "spm" in c_symbol:
+            libname = "libspm"
+            prefix  = "pyspm_"
+        else:
+            print("ERROR: function name without pastix nor spm")
+            return
+
+        # loop over the arguments to compose the different call lines
+        func_line = "def " + prefix + c_symbol + "("
+        func_line_length = len(func_line)
+        func_line_indent = len(func_line)
+
+        args_line = iindent*" " + libname + "." + c_symbol + ".argtypes = ["
+        args_line_length = len(args_line)
+        args_line_indent = len(args_line)
+
+        if is_function:
+            call_line = iindent*" " + "return " + libname + "." + c_symbol + "("
+
+            py_type = types_dict[return_type]
+            if return_pointer == "*":
+                if py_type == "c_void":
+                    py_type = "c_void_p"
+                else:
+                    py_type = "POINTER("+py_type+")"
+            elif  return_pointer == "**":
+                py_type = "c_void_p"
+
+            retv_line = iindent*" " + libname + "." + c_symbol + ".restype = " + py_type
+        else:
+            call_line = iindent*" " + libname + "." + c_symbol + "("
+            retv_line = ""
+        call_line_length = len(call_line)
+        call_line_indent = len(call_line)
+
+        args_list = []
+        args_size = [ 0, 0, 0 ]
+
+        # loop over the arguments to compose the first line
+        for j in range(1,len(function)):
+            iso_c_wrapper_type(function[j], args_list, args_size)
+
+        for j in range(0, len(args_list)):
+            # pointers
+            arg_name = args_list[j][0]
+            arg_type = args_list[j][1]
+            arg_call = args_list[j][2]
+
+            isfile = (len(arg_name) == 0)
+
+            if j > 0:
+                if not isfile:
+                    func_line += ","
+                    func_line_length += 2
+                args_line += ","
+                args_line_length += 2
+                call_line += ","
+                call_line_length += 2
+
+            # func_line
+            l = len(arg_name)
+            if not isfile:
+                if ((func_line_length + l) > 78):
+                    func_line_length = func_line_indent
+                    func_line += "\n" + func_line_indent*" "
+                func_line += " " + arg_name
+                func_line_length += l
+
+            # args_line
+            l = len(arg_type)
+            if ((args_line_length + l) > 78):
+                args_line_length = args_line_indent
+                args_line += "\n" + args_line_indent*" "
+            args_line += " " + arg_type
+            args_line_length += l
+
+            # call_line
+            l = len(arg_call)
+            if ((call_line_length + l) > 78):
+                call_line_length = call_line_indent
+                call_line += "\n" + call_line_indent*" "
+            call_line += " " + arg_call
+            call_line_length += l
+
+        py_interface  = func_line + " ):\n"
+        py_interface += args_line + " ]\n"
+        if len(retv_line) > 0:
+            py_interface += retv_line + "\n"
+        py_interface += call_line + " )\n"
+
+        # # add common header
+        # py_interface += indent + 2*indent + "use iso_c_binding\n"
+        # # import derived types
+        # for derived_type in used_derived_types:
+        #     py_interface += indent + 2*indent + "import " + derived_type +"\n"
+        # py_interface += indent + 2*indent + "implicit none\n"
+
+
+        # # add the return value of the function
+        # if (is_function):
+        #     py_interface +=  indent + 2*indent + iso_c_interface_type(function[0], True) + "_c"
+        #     py_interface += "\n"
+
+        # # loop over the arguments to describe them
+        # for j in range(1,len(function)):
+        #     py_interface += indent + 2*indent + iso_c_interface_type(function[j], False)
+        #     py_interface += "\n"
+
+        # if (is_function):
+        #     py_interface += indent + indent + "end function\n"
+        # else:
+        #     py_interface += indent + indent + "end subroutine\n"
+
+        # py_interface += indent + "end interface\n"
+
+        return py_interface
+
+    @staticmethod
+    def wrapper(function):
+        """Generate a wrapper for a function.
+           void functions in C will be called as subroutines,
+           functions in C will be turned to subroutines by appending
+           the return value as the last argument."""
+
+        # is it a function or a subroutine
+        if (function[0][0] == "void"):
+            is_function = False
+        else:
+            is_function = True
+
+        c_symbol = function[0][2]
+        f_symbol = c_symbol + "_c"
+
+        if (is_function):
+            initial_indent_signature = len(indent + "subroutine " + c_symbol + "(") * " "
+            initial_indent_call      = len(indent + indent + "info = " + f_symbol + "(") * " "
+        else:
+            initial_indent_signature = len(indent + "subroutine " + c_symbol + "(") * " "
+            initial_indent_call      = len(indent + indent + "call " + f_symbol + "(") * " "
+
+        # loop over the arguments to compose the first line and call line
+        signature_line = ""
+        call_line = ""
+        double_pointers = []
+        for j in range(1,len(function)):
+            if (j != 1):
+                signature_line += ", "
+                call_line += ", "
+
+            # do not make the argument list too long
+            if (j%9 == 0):
+                call_line      += "&\n" + initial_indent_call
+                signature_line += "&\n" + initial_indent_signature
+
+            # pointers
+            arg_type    = function[j][0]
+            arg_pointer = function[j][1]
+            arg_name    = function[j][2]
+
+            signature_line += arg_name
+            if (arg_pointer == "**"):
+                aux_name = arg_name + "_aux"
+                call_line += aux_name
+                double_pointers.append(arg_name)
+            elif (arg_pointer == "*"):
+                call_line += "c_loc(" + arg_name + ")"
+            else:
+                call_line += arg_name
+
+        contains_derived_types = False
+        for arg in function:
+            if (arg[0] in derived_types):
+                contains_derived_types = True
+
+        # initialize a string with the fortran interface
+        f_wrapper = ""
+        f_wrapper += indent + "subroutine "
+        f_wrapper += c_symbol + "("
+
+        # add the info argument at the end
+        f_wrapper += signature_line
+        if (is_function):
+            if (len(function) > 1):
+                f_wrapper += ", "
+
+            return_type = function[0][0]
+            return_pointer = function[0][1]
+            if (return_type == "int"):
+                return_var = "info"
+            else:
+                return_var = return_variables_dict[return_type]
+
+            f_wrapper += return_var
+
+        f_wrapper += ")\n"
+
+        # add common header
+        f_wrapper += indent + indent + "use iso_c_binding\n"
+        f_wrapper += indent + indent + "implicit none\n"
+
+        # loop over the arguments to describe them
+        for j in range(1,len(function)):
+            f_wrapper += indent + indent + iso_c_wrapper_type(function[j]) + "\n"
+
+        # add the return info value of the function
+        if (is_function):
+            if (function[0][1] == "*"):
+                f_target = ", pointer"
+            else:
+                f_target = ""
+
+            f_wrapper += indent + indent + types_dict[return_type] + ", intent(out)" + f_target + " :: " + return_var + "\n"
+
+        f_wrapper += "\n"
+
+        # loop over potential double pointers and generate auxiliary variables for them
+        for double_pointer in double_pointers:
+            aux_name = double_pointer + "_aux"
+            f_wrapper += indent + indent + "type(c_ptr) :: " + aux_name + "\n"
+            f_wrapper += "\n"
+
+        if (is_function):
+            f_return = return_var
+            f_return += " = "
+        else:
+            f_return = "call "
+
+        # generate the call to the C function
+        if (is_function and return_pointer == "*"):
+            f_wrapper += indent + indent + "call c_f_pointer(" + f_symbol + "(" + call_line + "), " + return_var + ")\n"
+        else:
+            f_wrapper += indent + indent + f_return + f_symbol + "(" + call_line + ")\n"
+
+        # loop over potential double pointers and translate them to Fortran pointers
+        for double_pointer in double_pointers:
+            aux_name = double_pointer + "_aux"
+            f_wrapper += indent + indent + "call c_f_pointer(" + aux_name + ", " + double_pointer + ")\n"
+
+        f_wrapper += indent + "end subroutine\n"
+
+        return f_wrapper
diff --git a/wrappers/CMakeLists.txt b/wrappers/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cf763161f61066f9b67c1f1e33f2ea7c39193b84
--- /dev/null
+++ b/wrappers/CMakeLists.txt
@@ -0,0 +1,20 @@
+###
+#
+#  @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                       Univ. Bordeaux. All rights reserved.
+#
+#  @version 6.0.0
+#  @author Mathieu Faverge
+#  @date 2017-05-22
+#
+###
+
+if (SPM_WITH_FORTRAN)
+  add_subdirectory( fortran90 )
+endif()
+
+if (BUILD_SHARED_LIBS)
+  add_subdirectory( python )
+else()
+  message(STATUS "--- Python wrapper is disabled with static libraries")
+endif()
diff --git a/wrappers/fortran90/CMakeLists.txt b/wrappers/fortran90/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..dd80adc3ae8ba3c1e84ffe8e6242a5335df6732a
--- /dev/null
+++ b/wrappers/fortran90/CMakeLists.txt
@@ -0,0 +1,56 @@
+###
+#
+#  @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                       Univ. Bordeaux. All rights reserved.
+#
+#  @version 6.0.0
+#  @author Mathieu Faverge
+#  @date 2017-05-22
+#
+###
+cmake_minimum_required (VERSION 3.1)
+
+# Coherce CMake to install the generated .mod files
+set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/mod_files)
+install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}/ DESTINATION include)
+
+add_library( spmf
+  src/spm_enums.F90
+  src/spmf.f90 )
+
+if ( SPM_INT64 )
+  set_source_files_properties(
+    src/spm_enums.F90
+    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=8")
+else()
+  set_source_files_properties(
+    src/spm_enums.F90
+    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=4")
+endif()
+
+target_link_libraries( spmf spm )
+install(TARGETS spmf
+  RUNTIME DESTINATION bin
+  ARCHIVE DESTINATION lib
+  LIBRARY DESTINATION lib )
+
+#
+# Add examples
+#
+set (EXAMPLES
+  spm_driver.f90
+  spm_user.f90
+  )
+
+foreach (_file ${EXAMPLES})
+  get_filename_component(_name_we ${_file} NAME_WE)
+  add_executable(${_name_we} examples/${_file})
+  target_link_libraries(${_name_we} spmf)
+
+  install(TARGETS ${_name_we}       RUNTIME DESTINATION examples )
+  install(FILES   examples/${_file}         DESTINATION examples )
+
+  add_test(fortran_${_name_we} ./${_name_we})
+
+endforeach()
+
diff --git a/wrappers/fortran90/examples/spm_driver.f90 b/wrappers/fortran90/examples/spm_driver.f90
new file mode 100644
index 0000000000000000000000000000000000000000..fc7f05e450898a1801c52a72aebcf95cd2121c48
--- /dev/null
+++ b/wrappers/fortran90/examples/spm_driver.f90
@@ -0,0 +1,70 @@
+!
+! @file spm_driver.f90
+!
+! Fortran 90 example using a matrix read with the spm driver.
+!
+! @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+!                      Univ. Bordeaux. All rights reserved.
+!
+! @version 6.0.0
+! @author Mathieu Faverge
+! @date 2017-01-01
+!
+program spm_driver
+  use iso_c_binding
+  use spmf
+  ! use mpi_f08
+  implicit none
+
+  type(spmatrix_t),           target                       :: spm
+  type(spmatrix_t),           pointer                      :: spm2
+  real(kind=c_double)                                      :: normA
+  real(kind=c_double)                                      :: eps = 1.e-15
+  integer(c_int)                                           :: info
+  integer(kind=spm_int_t)                                  :: nrhs
+  real(kind=c_double), dimension(:,:), allocatable, target :: x0, x, b
+  type(c_ptr)                                              :: x0_ptr, x_ptr, b_ptr
+
+  !
+  ! Initialize the problem
+  !   1- The matrix
+  call spmReadDriver( SpmDriverLaplacian, "d:10:10:10:4.", spm, info )
+
+  call spmCheckAndCorrect( spm, spm2 )
+  if (.not. c_associated(c_loc(spm), c_loc(spm2))) then
+     call spmExit( spm )
+     spm = spm2
+  end if
+
+  call spmPrintInfo( spm )
+
+  ! Scale A for better stability with low-rank computations
+  call spmNorm( SpmFrobeniusNorm, spm, normA )
+  call spmScalMatrix( 1. / normA, spm )
+
+  !   2- The right hand side
+  nrhs = 10
+  allocate(x0(spm%n, nrhs))
+  allocate(x( spm%n, nrhs))
+  allocate(b( spm%n, nrhs))
+  x0_ptr = c_loc(x0)
+  x_ptr  = c_loc(x)
+  b_ptr  = c_loc(b)
+
+  ! Compute b = A * x, with x random
+  call spmGenRHS( SpmRhsRndX, nrhs, spm, x0_ptr, spm%n, b_ptr, spm%n, info )
+
+  ! Copy x0 into x
+  x = x0
+
+  !
+  ! Check the solution
+  !
+  call spmCheckAxb( eps, nrhs, spm, x0_ptr, spm%n, b_ptr, spm%n, x_ptr, spm%n, info )
+
+  call spmExit( spm )
+  deallocate(x0)
+  deallocate(x)
+  deallocate(b)
+
+end program spm_driver
diff --git a/wrappers/fortran90/examples/spm_user.f90 b/wrappers/fortran90/examples/spm_user.f90
new file mode 100644
index 0000000000000000000000000000000000000000..cdf8df9bb84dd30ed8f0de8d9221f460661fb0a6
--- /dev/null
+++ b/wrappers/fortran90/examples/spm_user.f90
@@ -0,0 +1,156 @@
+!
+! @file spm_user.f90
+!
+! Fortran 90 example using a laplacian matrix.
+!
+! @copyright 2015-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+!                      Univ. Bordeaux. All rights reserved.
+!
+! @version 6.0.0
+! @author Mathieu Faverge
+! @date 2017-01-01
+!
+program spm_user
+  use iso_c_binding
+  use spmf
+  implicit none
+
+  integer(kind=spm_int_t), dimension(:), allocatable, target :: rowptr
+  integer(kind=spm_int_t), dimension(:), allocatable, target :: colptr
+  real(kind=c_double),  dimension(:), allocatable, target    :: values
+  real(kind=c_double),  dimension(:,:), allocatable, target  :: x0, x, b
+  type(c_ptr)                                                :: x0_ptr, x_ptr, b_ptr
+  real(kind=c_double)                                        :: eps = 1.e-15
+  type(spmatrix_t),        target                            :: spm
+  type(spmatrix_t),        pointer                           :: spm2
+  integer(kind=spm_int_t)                                    :: dim1, dim2, dim3, n, nnz
+  integer(kind=spm_int_t)                                    :: i, j, k, l, nrhs
+  integer(c_int)                                             :: info
+
+  !
+  ! Generate a 10x10x10 complex Laplacian in IJV format
+  !
+  dim1 = 10
+  dim2 = 10
+  dim3 = 10
+  n    = dim1 * dim2 * dim3
+  nnz  = (2*(dim1)-1) * dim2 * dim3 + (dim2-1)*dim1*dim3 + dim2*dim1*(dim3-1)
+
+  allocate(rowptr(nnz))
+  allocate(colptr(nnz))
+  allocate(values(nnz))
+
+  l = 1
+  do i=1,dim1
+     do j=1,dim2
+        do k=1,dim3
+           rowptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+           colptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+           values(l) = 6.
+
+           if (i == 1) then
+              values(l) = values(l) - 1.
+           end if
+           if (i == dim1) then
+              values(l) = values(l) - 1.
+           end if
+           if (j == 1) then
+              values(l) = values(l) - 1.
+           end if
+           if (j == dim2) then
+              values(l) = values(l) - 1.
+           end if
+           if (k == 1) then
+              values(l) = values(l) - 1.
+           end if
+           if (k == dim3) then
+              values(l) = values(l) - 1.
+           end if
+
+           values(l) = values(l) * 8.
+           l = l + 1
+
+           if (i < dim1) then
+              rowptr(l) =  i    + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+              colptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+              values(l) = - 1. - 1. * I
+              l = l + 1
+           end if
+           if (j < dim2) then
+              rowptr(l) = (i-1) + dim1 *  j    + dim1 * dim2 * (k-1) + 1
+              colptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+              values(l) = - 1. - 1. * I
+              l = l + 1
+           end if
+           if (k < dim3) then
+              rowptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 *  k    + 1
+              colptr(l) = (i-1) + dim1 * (j-1) + dim1 * dim2 * (k-1) + 1
+              values(l) = -1. - 1. * I
+              l = l + 1
+           end if
+        end do
+     end do
+  end do
+
+  if ( l .ne. nnz+1 ) then
+     write(6,*) 'l ', l, " nnz ", nnz
+  end if
+
+  !
+  ! Create the spm out of the internal data
+  !
+  call spmInit( spm )
+  spm%mtxtype = SpmSymmetric
+  spm%flttype = SpmDouble
+  spm%fmttype = SpmIJV
+  spm%n       = n
+  spm%nnz     = nnz
+  spm%dof     = 1
+  spm%rowptr  = c_loc(rowptr)
+  spm%colptr  = c_loc(colptr)
+  spm%values  = c_loc(values)
+
+  call spmUpdateComputedFields( spm )
+
+  call spmCheckAndCorrect( spm, spm2 )
+  if (.not. c_associated(c_loc(spm), c_loc(spm2))) then
+     deallocate(rowptr)
+     deallocate(colptr)
+     deallocate(values)
+
+     spm%rowptr = c_null_ptr
+     spm%colptr = c_null_ptr
+     spm%values = c_null_ptr
+
+     call spmExit( spm )
+     spm = spm2
+  end if
+
+  call spmPrintInfo( spm )
+
+  !   2- The right hand side
+  nrhs = 10
+  allocate(x0(spm%n, nrhs))
+  allocate(x( spm%n, nrhs))
+  allocate(b( spm%n, nrhs))
+  x0_ptr = c_loc(x0)
+  x_ptr  = c_loc(x)
+  b_ptr  = c_loc(b)
+
+  ! Compute b = A * x, with x random
+  call spmGenRHS( SpmRhsRndX, nrhs, spm, x0_ptr, spm%n, b_ptr, spm%n, info )
+
+  ! Copy x0 into x
+  x = x0
+
+  !
+  ! Check the solution
+  !
+  call spmCheckAxb( eps, nrhs, spm, x0_ptr, spm%n, b_ptr, spm%n, x_ptr, spm%n, info )
+
+  call spmExit( spm )
+  deallocate(x0)
+  deallocate(x)
+  deallocate(b)
+
+end program spm_user
diff --git a/wrappers/fortran90/src/spm_enums.F90 b/wrappers/fortran90/src/spm_enums.F90
new file mode 100644
index 0000000000000000000000000000000000000000..5b7f9c93d6cb68a2aa2c0b5f54e95b03a24d5330
--- /dev/null
+++ b/wrappers/fortran90/src/spm_enums.F90
@@ -0,0 +1,143 @@
+
+!
+! @file spm_enums.F90
+!
+! SPM fortran 90 wrapper to define enums and datatypes
+!
+! @copyright 2017      Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+!                      Univ. Bordeaux. All rights reserved.
+!
+! @version 6.0.0
+! @author Mathieu Faverge
+! @date 2017-01-01
+!
+! This file has been automatically generated with gen_wrappers.py
+!
+module spm_enums
+  use iso_c_binding
+  implicit none
+
+  ! enum verbose
+  enum, bind(C)
+     enumerator :: SpmVerboseNot = 0
+     enumerator :: SpmVerboseNo  = 1
+     enumerator :: SpmVerboseYes = 2
+  end enum
+
+  ! enum coeftype
+  enum, bind(C)
+     enumerator :: SpmPattern   = 0
+     enumerator :: SpmFloat     = 2
+     enumerator :: SpmDouble    = 3
+     enumerator :: SpmComplex32 = 4
+     enumerator :: SpmComplex64 = 5
+  end enum
+
+  ! enum fmttype
+  enum, bind(C)
+     enumerator :: SpmCSC = 0
+     enumerator :: SpmCSR = 1
+     enumerator :: SpmIJV = 2
+  end enum
+
+  ! enum error
+  enum, bind(C)
+     enumerator :: SPM_SUCCESS            = 0
+     enumerator :: SPM_ERR_UNKNOWN        = 1
+     enumerator :: SPM_ERR_ALLOC          = 2
+     enumerator :: SPM_ERR_NOTIMPLEMENTED = 3
+     enumerator :: SPM_ERR_OUTOFMEMORY    = 4
+     enumerator :: SPM_ERR_THREAD         = 5
+     enumerator :: SPM_ERR_INTERNAL       = 6
+     enumerator :: SPM_ERR_BADPARAMETER   = 7
+     enumerator :: SPM_ERR_FILE           = 8
+     enumerator :: SPM_ERR_INTEGER_TYPE   = 9
+     enumerator :: SPM_ERR_IO             = 10
+     enumerator :: SPM_ERR_MPI            = 11
+  end enum
+
+  ! enum driver
+  enum, bind(C)
+     enumerator :: SpmDriverRSA        = 0
+     enumerator :: SpmDriverHB         = 1
+     enumerator :: SpmDriverIJV        = 2
+     enumerator :: SpmDriverMM         = 3
+     enumerator :: SpmDriverLaplacian  = 4
+     enumerator :: SpmDriverXLaplacian = 5
+     enumerator :: SpmDriverGraph      = 6
+     enumerator :: SpmDriverSPM        = 7
+  end enum
+
+  ! enum rhstype
+  enum, bind(C)
+     enumerator :: SpmRhsOne  = 0
+     enumerator :: SpmRhsI    = 1
+     enumerator :: SpmRhsRndX = 2
+     enumerator :: SpmRhsRndB = 3
+  end enum
+
+  ! enum layout
+  enum, bind(C)
+     enumerator :: SpmRowMajor = 101
+     enumerator :: SpmColMajor = 102
+  end enum
+
+  ! enum trans
+  enum, bind(C)
+     enumerator :: SpmNoTrans   = 111
+     enumerator :: SpmTrans     = 112
+     enumerator :: SpmConjTrans = 113
+  end enum
+
+  ! enum mtxtype
+  enum, bind(C)
+     enumerator :: SpmGeneral   = SpmNoTrans
+     enumerator :: SpmSymmetric = SpmTrans
+     enumerator :: SpmHermitian = SpmConjTrans
+  end enum
+
+  ! enum uplo
+  enum, bind(C)
+     enumerator :: SpmUpper      = 121
+     enumerator :: SpmLower      = 122
+     enumerator :: SpmUpperLower = 123
+  end enum
+
+  ! enum diag
+  enum, bind(C)
+     enumerator :: SpmNonUnit = 131
+     enumerator :: SpmUnit    = 132
+  end enum
+
+  ! enum side
+  enum, bind(C)
+     enumerator :: SpmLeft  = 141
+     enumerator :: SpmRight = 142
+  end enum
+
+  ! enum normtype
+  enum, bind(C)
+     enumerator :: SpmOneNorm       = 171
+     enumerator :: SpmFrobeniusNorm = 174
+     enumerator :: SpmInfNorm       = 175
+     enumerator :: SpmMaxNorm       = 177
+  end enum
+
+  ! enum dir
+  enum, bind(C)
+     enumerator :: SpmDirForward  = 391
+     enumerator :: SpmDirBackward = 392
+  end enum
+
+
+  integer, parameter :: spm_int_t = SPM_INT_KIND
+
+contains
+
+  function spm_getintsize()
+    integer :: spm_getintsize
+    spm_getintsize = SPM_INT_KIND
+    return
+  end function spm_getintsize
+
+end module spm_enums
diff --git a/wrappers/fortran90/src/spmf.f90 b/wrappers/fortran90/src/spmf.f90
new file mode 100644
index 0000000000000000000000000000000000000000..2b552dd50832ba8ff1fc45c23d0162eb1504fa6d
--- /dev/null
+++ b/wrappers/fortran90/src/spmf.f90
@@ -0,0 +1,711 @@
+
+!
+! @file spmf.f90
+!
+! SPM Fortran 90 wrapper
+!
+! @copyright 2017      Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+!                      Univ. Bordeaux. All rights reserved.
+!
+! @version 6.0.0
+! @author Mathieu Faverge
+! @date 2017-01-01
+!
+! This file has been automatically generated with gen_wrappers.py
+!
+module spmf
+  use iso_c_binding
+  use spm_enums
+  implicit none
+
+  type, bind(c) :: spmatrix_t
+     integer(c_int)          :: mtxtype
+     integer(c_int)          :: flttype
+     integer(c_int)          :: fmttype
+     integer(kind=spm_int_t) :: gN
+     integer(kind=spm_int_t) :: n
+     integer(kind=spm_int_t) :: gnnz
+     integer(kind=spm_int_t) :: nnz
+     integer(kind=spm_int_t) :: gNexp
+     integer(kind=spm_int_t) :: nexp
+     integer(kind=spm_int_t) :: gnnzexp
+     integer(kind=spm_int_t) :: nnzexp
+     integer(kind=spm_int_t) :: dof
+     type(c_ptr)             :: dofs
+     integer(c_int)          :: layout
+     type(c_ptr)             :: colptr
+     type(c_ptr)             :: rowptr
+     type(c_ptr)             :: loc2glob
+     type(c_ptr)             :: values
+  end type spmatrix_t
+
+  interface
+     subroutine spmInit_c(spm) &
+          bind(c, name='spmInit')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spmInit_c
+  end interface
+
+  interface
+     subroutine spmExit_c(spm) &
+          bind(c, name='spmExit')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spmExit_c
+  end interface
+
+  interface
+     function spmCopy_c(spm) &
+          bind(c, name='spmCopy')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr)        :: spmCopy_c
+       type(c_ptr), value :: spm
+     end function spmCopy_c
+  end interface
+
+  interface
+     subroutine spmBase_c(spm, baseval) &
+          bind(c, name='spmBase')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr),         value :: spm
+       integer(kind=c_int), value :: baseval
+     end subroutine spmBase_c
+  end interface
+
+  interface
+     function spmFindBase_c(spm) &
+          bind(c, name='spmFindBase')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=spm_int_t)   :: spmFindBase_c
+       type(c_ptr),        value :: spm
+     end function spmFindBase_c
+  end interface
+
+  interface
+     function spmConvert_c(ofmttype, ospm) &
+          bind(c, name='spmConvert')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)        :: spmConvert_c
+       integer(kind=c_int), value :: ofmttype
+       type(c_ptr),         value :: ospm
+     end function spmConvert_c
+  end interface
+
+  interface
+     subroutine spmUpdateComputedFields_c(spm) &
+          bind(c, name='spmUpdateComputedFields')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spmUpdateComputedFields_c
+  end interface
+
+  interface
+     subroutine spmGenFakeValues_c(spm) &
+          bind(c, name='spmGenFakeValues')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spmGenFakeValues_c
+  end interface
+
+  interface
+     function spmNorm_c(ntype, spm) &
+          bind(c, name='spmNorm')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       real(kind=c_double)   :: spmNorm_c
+       integer(c_int), value :: ntype
+       type(c_ptr),    value :: spm
+     end function spmNorm_c
+  end interface
+
+  interface
+     function spmMatVec_c(trans, alpha, spm, x, beta, y) &
+          bind(c, name='spmMatVec')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)   :: spmMatVec_c
+       integer(c_int), value :: trans
+       type(c_ptr),    value :: alpha
+       type(c_ptr),    value :: spm
+       type(c_ptr),    value :: x
+       type(c_ptr),    value :: beta
+       type(c_ptr),    value :: y
+     end function spmMatVec_c
+  end interface
+
+  interface
+     function spmMatMat_c(trans, n, alpha, A, B, ldb, beta, C, ldc) &
+          bind(c, name='spmMatMat')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)            :: spmMatMat_c
+       integer(c_int),          value :: trans
+       integer(kind=spm_int_t), value :: n
+       type(c_ptr),             value :: alpha
+       type(c_ptr),             value :: A
+       type(c_ptr),             value :: B
+       integer(kind=spm_int_t), value :: ldb
+       type(c_ptr),             value :: beta
+       type(c_ptr),             value :: C
+       integer(kind=spm_int_t), value :: ldc
+     end function spmMatMat_c
+  end interface
+
+  interface
+     subroutine spmScalMatrix_c(alpha, spm) &
+          bind(c, name='spmScalMatrix')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       real(kind=c_double), value :: alpha
+       type(c_ptr),         value :: spm
+     end subroutine spmScalMatrix_c
+  end interface
+
+  interface
+     subroutine spmScalVector_c(flt, alpha, n, x, incx) &
+          bind(c, name='spmScalVector')
+       use iso_c_binding
+       import spm_int_t
+       implicit none
+       integer(c_int),          value :: flt
+       real(kind=c_double),     value :: alpha
+       integer(kind=spm_int_t), value :: n
+       type(c_ptr),             value :: x
+       integer(kind=spm_int_t), value :: incx
+     end subroutine spmScalVector_c
+  end interface
+
+  interface
+     function spmSort_c(spm) &
+          bind(c, name='spmSort')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)   :: spmSort_c
+       type(c_ptr),    value :: spm
+     end function spmSort_c
+  end interface
+
+  interface
+     function spmMergeDuplicate_c(spm) &
+          bind(c, name='spmMergeDuplicate')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=spm_int_t)   :: spmMergeDuplicate_c
+       type(c_ptr),        value :: spm
+     end function spmMergeDuplicate_c
+  end interface
+
+  interface
+     function spmSymmetrize_c(spm) &
+          bind(c, name='spmSymmetrize')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=spm_int_t)   :: spmSymmetrize_c
+       type(c_ptr),        value :: spm
+     end function spmSymmetrize_c
+  end interface
+
+  interface
+     function spmCheckAndCorrect_c(spm) &
+          bind(c, name='spmCheckAndCorrect')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr)        :: spmCheckAndCorrect_c
+       type(c_ptr), value :: spm
+     end function spmCheckAndCorrect_c
+  end interface
+
+  interface
+     function spmGenRHS_c(type, nrhs, spm, x, ldx, b, ldb) &
+          bind(c, name='spmGenRHS')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)            :: spmGenRHS_c
+       integer(c_int),          value :: type
+       integer(kind=spm_int_t), value :: nrhs
+       type(c_ptr),             value :: spm
+       type(c_ptr),             value :: x
+       integer(kind=spm_int_t), value :: ldx
+       type(c_ptr),             value :: b
+       integer(kind=spm_int_t), value :: ldb
+     end function spmGenRHS_c
+  end interface
+
+  interface
+     function spmCheckAxb_c(eps, nrhs, spm, x0, ldx0, b, ldb, x, ldx) &
+          bind(c, name='spmCheckAxb')
+       use iso_c_binding
+       import spm_int_t
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)            :: spmCheckAxb_c
+       real(kind=c_double),     value :: eps
+       integer(kind=spm_int_t), value :: nrhs
+       type(c_ptr),             value :: spm
+       type(c_ptr),             value :: x0
+       integer(kind=spm_int_t), value :: ldx0
+       type(c_ptr),             value :: b
+       integer(kind=spm_int_t), value :: ldb
+       type(c_ptr),             value :: x
+       integer(kind=spm_int_t), value :: ldx
+     end function spmCheckAxb_c
+  end interface
+
+  interface
+     function spmIntConvert_c(n, input) &
+          bind(c, name='spmIntConvert')
+       use iso_c_binding
+       import spm_int_t
+       implicit none
+       type(c_ptr)                    :: spmIntConvert_c
+       integer(kind=spm_int_t), value :: n
+       type(c_ptr),             value :: input
+     end function spmIntConvert_c
+  end interface
+
+  interface
+     function spmLoad_c(spm, infile) &
+          bind(c, name='spmLoad')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)   :: spmLoad_c
+       type(c_ptr),    value :: spm
+       type(c_ptr),    value :: infile
+     end function spmLoad_c
+  end interface
+
+  interface
+     function spmSave_c(spm, outfile) &
+          bind(c, name='spmSave')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)   :: spmSave_c
+       type(c_ptr),    value :: spm
+       type(c_ptr),    value :: outfile
+     end function spmSave_c
+  end interface
+
+  interface
+     function spmReadDriver_c(driver, filename, spm) &
+          bind(c, name='spmReadDriver')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       integer(kind=c_int)   :: spmReadDriver_c
+       integer(c_int), value :: driver
+       type(c_ptr),    value :: filename
+       type(c_ptr),    value :: spm
+     end function spmReadDriver_c
+  end interface
+
+  interface
+     function spmParseLaplacianInfo_c(filename, flttype, dim1, dim2, dim3, &
+          alpha, beta) &
+          bind(c, name='spmParseLaplacianInfo')
+       use iso_c_binding
+       import spm_int_t
+       implicit none
+       integer(kind=c_int)   :: spmParseLaplacianInfo_c
+       type(c_ptr),    value :: filename
+       type(c_ptr),    value :: flttype
+       type(c_ptr),    value :: dim1
+       type(c_ptr),    value :: dim2
+       type(c_ptr),    value :: dim3
+       type(c_ptr),    value :: alpha
+       type(c_ptr),    value :: beta
+     end function spmParseLaplacianInfo_c
+  end interface
+
+  interface
+     subroutine spm2Dense_c(spm) &
+          bind(c, name='spm2Dense')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spm2Dense_c
+  end interface
+
+  interface
+     subroutine spmPrint_c(spm, f) &
+          bind(c, name='spmPrint')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+       type(c_ptr), value :: f
+     end subroutine spmPrint_c
+  end interface
+
+  interface
+     subroutine spmPrintInfo_c(spm, f) &
+          bind(c, name='spmPrintInfo')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+       type(c_ptr), value :: f
+     end subroutine spmPrintInfo_c
+  end interface
+
+  interface
+     function spmExpand_c(spm) &
+          bind(c, name='spmExpand')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr)        :: spmExpand_c
+       type(c_ptr), value :: spm
+     end function spmExpand_c
+  end interface
+
+  interface
+     function spmDofExtend_c(spm, type, dof) &
+          bind(c, name='spmDofExtend')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr)                :: spmDofExtend_c
+       type(c_ptr),         value :: spm
+       integer(kind=c_int), value :: type
+       integer(kind=c_int), value :: dof
+     end function spmDofExtend_c
+  end interface
+
+contains
+
+  ! Wrappers of the C functions.
+  subroutine spmInit(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target :: spm
+
+    call spmInit_c(c_loc(spm))
+  end subroutine spmInit
+
+  subroutine spmExit(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target :: spm
+
+    call spmExit_c(c_loc(spm))
+  end subroutine spmExit
+
+  subroutine spmCopy(spm, spmo)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(in),  target  :: spm
+    type(spmatrix_t), intent(out), pointer :: spmo
+
+    call c_f_pointer(spmCopy_c(c_loc(spm)), spmo)
+  end subroutine spmCopy
+
+  subroutine spmBase(spm, baseval)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),    intent(inout), target :: spm
+    integer(kind=c_int), intent(in)            :: baseval
+
+    call spmBase_c(c_loc(spm), baseval)
+  end subroutine spmBase
+
+  subroutine spmFindBase(spm, value)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),        intent(in), target :: spm
+    integer(kind=spm_int_t), intent(out)        :: value
+
+    value = spmFindBase_c(c_loc(spm))
+  end subroutine spmFindBase
+
+  subroutine spmConvert(ofmttype, ospm, info)
+    use iso_c_binding
+    implicit none
+    integer(kind=c_int), intent(in)            :: ofmttype
+    type(spmatrix_t),    intent(inout), target :: ospm
+    integer(kind=c_int), intent(out)           :: info
+
+    info = spmConvert_c(ofmttype, c_loc(ospm))
+  end subroutine spmConvert
+
+  subroutine spmUpdateComputedFields(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target :: spm
+
+    call spmUpdateComputedFields_c(c_loc(spm))
+  end subroutine spmUpdateComputedFields
+
+  subroutine spmGenFakeValues(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target :: spm
+
+    call spmGenFakeValues_c(c_loc(spm))
+  end subroutine spmGenFakeValues
+
+  subroutine spmNorm(ntype, spm, value)
+    use iso_c_binding
+    implicit none
+    integer(c_int),      intent(in)         :: ntype
+    type(spmatrix_t),    intent(in), target :: spm
+    real(kind=c_double), intent(out)        :: value
+
+    value = spmNorm_c(ntype, c_loc(spm))
+  end subroutine spmNorm
+
+  subroutine spmMatVec(trans, alpha, spm, x, beta, y, info)
+    use iso_c_binding
+    implicit none
+    integer(c_int),      intent(in)            :: trans
+    type(c_ptr),         intent(in),    target :: alpha
+    type(spmatrix_t),    intent(in),    target :: spm
+    type(c_ptr),         intent(in),    target :: x
+    type(c_ptr),         intent(in),    target :: beta
+    type(c_ptr),         intent(inout), target :: y
+    integer(kind=c_int), intent(out)           :: info
+
+    info = spmMatVec_c(trans, alpha, c_loc(spm), x, beta, y)
+  end subroutine spmMatVec
+
+  subroutine spmMatMat(trans, n, alpha, A, B, ldb, beta, C, ldc, info)
+    use iso_c_binding
+    implicit none
+    integer(c_int),          intent(in)            :: trans
+    integer(kind=spm_int_t), intent(in)            :: n
+    type(c_ptr),             intent(in),    target :: alpha
+    type(spmatrix_t),        intent(in),    target :: A
+    type(c_ptr),             intent(in),    target :: B
+    integer(kind=spm_int_t), intent(in)            :: ldb
+    type(c_ptr),             intent(in),    target :: beta
+    type(c_ptr),             intent(inout), target :: C
+    integer(kind=spm_int_t), intent(in)            :: ldc
+    integer(kind=c_int),     intent(out)           :: info
+
+    info = spmMatMat_c(trans, n, alpha, c_loc(A), B, ldb, beta, C, ldc)
+  end subroutine spmMatMat
+
+  subroutine spmScalMatrix(alpha, spm)
+    use iso_c_binding
+    implicit none
+    real(kind=c_double), intent(in)            :: alpha
+    type(spmatrix_t),    intent(inout), target :: spm
+
+    call spmScalMatrix_c(alpha, c_loc(spm))
+  end subroutine spmScalMatrix
+
+  subroutine spmScalVector(flt, alpha, n, x, incx)
+    use iso_c_binding
+    implicit none
+    integer(c_int),          intent(in)            :: flt
+    real(kind=c_double),     intent(in)            :: alpha
+    integer(kind=spm_int_t), intent(in)            :: n
+    type(c_ptr),             intent(inout), target :: x
+    integer(kind=spm_int_t), intent(in)            :: incx
+
+    call spmScalVector_c(flt, alpha, n, x, incx)
+  end subroutine spmScalVector
+
+  subroutine spmSort(spm, info)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),    intent(inout), target :: spm
+    integer(kind=c_int), intent(out)           :: info
+
+    info = spmSort_c(c_loc(spm))
+  end subroutine spmSort
+
+  subroutine spmMergeDuplicate(spm, value)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),        intent(inout), target :: spm
+    integer(kind=spm_int_t), intent(out)           :: value
+
+    value = spmMergeDuplicate_c(c_loc(spm))
+  end subroutine spmMergeDuplicate
+
+  subroutine spmSymmetrize(spm, value)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),        intent(inout), target :: spm
+    integer(kind=spm_int_t), intent(out)           :: value
+
+    value = spmSymmetrize_c(c_loc(spm))
+  end subroutine spmSymmetrize
+
+  subroutine spmCheckAndCorrect(spm, spmo)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target  :: spm
+    type(spmatrix_t), intent(out),   pointer :: spmo
+
+    call c_f_pointer(spmCheckAndCorrect_c(c_loc(spm)), spmo)
+  end subroutine spmCheckAndCorrect
+
+  subroutine spmGenRHS(type, nrhs, spm, x, ldx, b, ldb, info)
+    use iso_c_binding
+    implicit none
+    integer(c_int),          intent(in)            :: type
+    integer(kind=spm_int_t), intent(in)            :: nrhs
+    type(spmatrix_t),        intent(in),    target :: spm
+    type(c_ptr),             intent(inout), target :: x
+    integer(kind=spm_int_t), intent(in)            :: ldx
+    type(c_ptr),             intent(inout), target :: b
+    integer(kind=spm_int_t), intent(in)            :: ldb
+    integer(kind=c_int),     intent(out)           :: info
+
+    info = spmGenRHS_c(type, nrhs, c_loc(spm), x, ldx, b, ldb)
+  end subroutine spmGenRHS
+
+  subroutine spmCheckAxb(eps, nrhs, spm, x0, ldx0, b, ldb, x, ldx, info)
+    use iso_c_binding
+    implicit none
+    real(kind=c_double),     intent(in)            :: eps
+    integer(kind=spm_int_t), intent(in)            :: nrhs
+    type(spmatrix_t),        intent(in),    target :: spm
+    type(c_ptr),             intent(inout), target :: x0
+    integer(kind=spm_int_t), intent(in)            :: ldx0
+    type(c_ptr),             intent(inout), target :: b
+    integer(kind=spm_int_t), intent(in)            :: ldb
+    type(c_ptr),             intent(in),    target :: x
+    integer(kind=spm_int_t), intent(in)            :: ldx
+    integer(kind=c_int),     intent(out)           :: info
+
+    info = spmCheckAxb_c(eps, nrhs, c_loc(spm), x0, ldx0, b, ldb, x, ldx)
+  end subroutine spmCheckAxb
+
+  subroutine spmIntConvert(n, input, value)
+    use iso_c_binding
+    implicit none
+    integer(kind=spm_int_t), intent(in)             :: n
+    integer(kind=c_int),     intent(inout), target  :: input
+    integer(kind=spm_int_t), intent(out),   pointer :: value
+
+    call c_f_pointer(spmIntConvert_c(n, c_loc(input)), value)
+  end subroutine spmIntConvert
+
+  subroutine spmLoad(spm, info)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),    intent(inout), target :: spm
+    integer(kind=c_int), intent(out)           :: info
+
+    info = spmLoad_c(c_loc(spm), c_null_ptr)
+  end subroutine spmLoad
+
+  subroutine spmSave(spm, info)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),    intent(in), target :: spm
+    integer(kind=c_int), intent(out)        :: info
+
+    info = spmSave_c(c_loc(spm), c_null_ptr)
+  end subroutine spmSave
+
+  subroutine spmReadDriver(driver, filename, spm, info)
+    use iso_c_binding
+    implicit none
+    integer(c_int),         intent(in)            :: driver
+    character(kind=c_char), intent(in),    target :: filename
+    type(spmatrix_t),       intent(inout), target :: spm
+    integer(kind=c_int),    intent(out)           :: info
+
+    info = spmReadDriver_c(driver, c_loc(filename), c_loc(spm))
+  end subroutine spmReadDriver
+
+  subroutine spmParseLaplacianInfo(filename, flttype, dim1, dim2, dim3, alpha, &
+       beta, info)
+    use iso_c_binding
+    implicit none
+    character(kind=c_char),  intent(in),    target :: filename
+    integer(c_int),          intent(inout), target :: flttype
+    integer(kind=spm_int_t), intent(inout), target :: dim1
+    integer(kind=spm_int_t), intent(inout), target :: dim2
+    integer(kind=spm_int_t), intent(inout), target :: dim3
+    real(kind=c_double),     intent(inout), target :: alpha
+    real(kind=c_double),     intent(inout), target :: beta
+    integer(kind=c_int),     intent(out)           :: info
+
+    info = spmParseLaplacianInfo_c(c_loc(filename), c_loc(flttype), &
+         c_loc(dim1), c_loc(dim2), c_loc(dim3), c_loc(alpha), c_loc(beta))
+  end subroutine spmParseLaplacianInfo
+
+  subroutine spm2Dense(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(in), target :: spm
+
+    call spm2Dense_c(c_loc(spm))
+  end subroutine spm2Dense
+
+  subroutine spmPrint(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(in), target :: spm
+
+    call spmPrint_c(c_loc(spm), c_null_ptr)
+  end subroutine spmPrint
+
+  subroutine spmPrintInfo(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(in), target :: spm
+
+    call spmPrintInfo_c(c_loc(spm), c_null_ptr)
+  end subroutine spmPrintInfo
+
+  subroutine spmExpand(spm, spmo)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(in),  target  :: spm
+    type(spmatrix_t), intent(out), pointer :: spmo
+
+    call c_f_pointer(spmExpand_c(c_loc(spm)), spmo)
+  end subroutine spmExpand
+
+  subroutine spmDofExtend(spm, type, dof, spmo)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t),    intent(in),  target  :: spm
+    integer(kind=c_int), intent(in)           :: type
+    integer(kind=c_int), intent(in)           :: dof
+    type(spmatrix_t),    intent(out), pointer :: spmo
+
+    call c_f_pointer(spmDofExtend_c(c_loc(spm), type, dof), spmo)
+  end subroutine spmDofExtend
+
+
+end module spmf
diff --git a/wrappers/python/CMakeLists.txt b/wrappers/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5266b5086b97cadc26015e44b7070559033eca45
--- /dev/null
+++ b/wrappers/python/CMakeLists.txt
@@ -0,0 +1,48 @@
+###
+#
+#  @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                       Univ. Bordeaux. All rights reserved.
+#
+#  @version 6.0.0
+#  @author Mathieu Faverge
+#  @date 2018-05-14
+#
+###
+
+# Configure enum.py
+if (SPM_INT64)
+  set(SPM_PYTHON_INTEGER c_int64)
+else()
+  set(SPM_PYTHON_INTEGER c_int)
+endif()
+
+configure_file(
+  "${CMAKE_CURRENT_SOURCE_DIR}/spm/enum.py.in"
+  "${CMAKE_CURRENT_SOURCE_DIR}/spm/enum.py" @ONLY)
+
+# Install python package
+install(FILES
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm/__init__.py
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm/__spm__.py
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm/spm.py
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm/enum.py
+  DESTINATION lib/python/spm )
+
+# Install python examples
+install(FILES
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm_driver.py
+  ${CMAKE_CURRENT_SOURCE_DIR}/spm_scipy.py
+  DESTINATION examples
+  )
+
+## CTest execution
+find_package(PythonInterp QUIET)
+if (PYTHONINTERP_FOUND)
+  set( PYTHON_TESTS
+    spm_driver spm_scipy )
+
+  foreach(example ${PYTHON_TESTS} )
+    add_test(python_${example} ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/${example}.py)
+  endforeach()
+endif()
+
diff --git a/wrappers/python/spm/__init__.py b/wrappers/python/spm/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..caf07917da51da1e91b27984eaa0d30c98c0d3bf
--- /dev/null
+++ b/wrappers/python/spm/__init__.py
@@ -0,0 +1,48 @@
+#
+# @file __init__.py
+#
+# SParse Matrix package python module intialization
+#
+# @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                      Univ. Bordeaux. All rights reserved.
+#
+# @version 6.0.0
+# @author Pierre Ramet
+# @author Mathieu Faverge
+# @author Louis Poirel
+# @date 2017-05-04
+#
+"""
+PySpm
+=====
+
+Provides
+  1. A sparse matrix structure
+  2. Mathematical operations over this sparse matrix structure
+  3. Driver to read from different file formats and to convert from Scipy package
+
+"""
+import ctypes
+import ctypes.util
+
+# Load the SPM library
+libspm_name = ctypes.util.find_library('spm')
+if libspm_name == None:
+    raise EnvironmentError("Could not find shared library: spm."
+                           "The path to libpastix_spm.so should be in "
+                           "$LIBRARY_PATH")
+
+try:
+    libspm = ctypes.cdll.LoadLibrary(libspm_name)
+except:
+    raise EnvironmentError("Could not load shared library: spm."
+                           "The path to libpastix_spm.so should be in "
+                           "$LD_LIBRARY_PATH or $DYLD_LIBRARY_PATH on MacOS");
+
+__all__ = [ 'libspm' ]
+
+from .enum   import *
+from .spm    import *
+
+#__all__.extend(enum.__all__)
+#__all__.extend(spm.__all__)
diff --git a/wrappers/python/spm/__spm__.py b/wrappers/python/spm/__spm__.py
new file mode 100644
index 0000000000000000000000000000000000000000..a575b2f7ec77c7a9aeded96b868280a4bd27d25c
--- /dev/null
+++ b/wrappers/python/spm/__spm__.py
@@ -0,0 +1,196 @@
+"""
+
+ @file __spm__.py
+
+ SPM python wrapper
+
+ @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @author Louis Poirel
+ @date 2017-05-04
+
+This file has been automatically generated with gen_wrappers.py
+
+"""
+from ctypes import *
+import numpy as np
+
+from . import libspm
+from .enum import __spm_int__
+
+class pyspm_spmatrix_t(Structure):
+    _fields_ = [("mtxtype",   c_int               ),
+                ("flttype",   c_int               ),
+                ("fmttype",   c_int               ),
+                ("gN",        __spm_int__         ),
+                ("n",         __spm_int__         ),
+                ("gnnz",      __spm_int__         ),
+                ("nnz",       __spm_int__         ),
+                ("gNexp",     __spm_int__         ),
+                ("nexp",      __spm_int__         ),
+                ("gnnzexp",   __spm_int__         ),
+                ("nnzexp",    __spm_int__         ),
+                ("dof",       __spm_int__         ),
+                ("dofs",      POINTER(__spm_int__)),
+                ("layout",    c_int               ),
+                ("colptr",    POINTER(__spm_int__)),
+                ("rowptr",    POINTER(__spm_int__)),
+                ("loc2glob",  POINTER(__spm_int__)),
+                ("values",    c_void_p            ) ]
+
+def pyspm_spmInit( spm ):
+    libspm.spmInit.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmInit( spm )
+
+def pyspm_spmExit( spm ):
+    libspm.spmExit.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmExit( spm )
+
+def pyspm_spmCopy( spm ):
+    libspm.spmCopy.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmCopy.restype = POINTER(pyspm_spmatrix_t)
+    return libspm.spmCopy( spm )
+
+def pyspm_spmBase( spm, baseval ):
+    libspm.spmBase.argtypes = [ POINTER(pyspm_spmatrix_t), c_int ]
+    libspm.spmBase( spm, baseval )
+
+def pyspm_spmFindBase( spm ):
+    libspm.spmFindBase.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmFindBase.restype = __spm_int__
+    return libspm.spmFindBase( spm )
+
+def pyspm_spmConvert( ofmttype, ospm ):
+    libspm.spmConvert.argtypes = [ c_int, POINTER(pyspm_spmatrix_t) ]
+    libspm.spmConvert.restype = c_int
+    return libspm.spmConvert( ofmttype, ospm )
+
+def pyspm_spmUpdateComputedFields( spm ):
+    libspm.spmUpdateComputedFields.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmUpdateComputedFields( spm )
+
+def pyspm_spmGenFakeValues( spm ):
+    libspm.spmGenFakeValues.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmGenFakeValues( spm )
+
+def pyspm_spmNorm( ntype, spm ):
+    libspm.spmNorm.argtypes = [ c_int, POINTER(pyspm_spmatrix_t) ]
+    libspm.spmNorm.restype = c_double
+    return libspm.spmNorm( ntype, spm )
+
+def pyspm_spmMatVec( trans, alpha, spm, x, beta, y ):
+    libspm.spmMatVec.argtypes = [ c_int, c_void_p, POINTER(pyspm_spmatrix_t),
+                                  c_void_p, c_void_p, c_void_p ]
+    libspm.spmMatVec.restype = c_int
+    return libspm.spmMatVec( trans, alpha, spm, x, beta, y )
+
+def pyspm_spmMatMat( trans, n, alpha, A, B, ldb, beta, C, ldc ):
+    libspm.spmMatMat.argtypes = [ c_int, __spm_int__, c_void_p,
+                                  POINTER(pyspm_spmatrix_t), c_void_p,
+                                  __spm_int__, c_void_p, c_void_p, __spm_int__ ]
+    libspm.spmMatMat.restype = c_int
+    return libspm.spmMatMat( trans, n, alpha, A, B, ldb, beta, C, ldc )
+
+def pyspm_spmScalMatrix( alpha, spm ):
+    libspm.spmScalMatrix.argtypes = [ c_double, POINTER(pyspm_spmatrix_t) ]
+    libspm.spmScalMatrix( alpha, spm )
+
+def pyspm_spmScalVector( flt, alpha, n, x, incx ):
+    libspm.spmScalVector.argtypes = [ c_int, c_double, __spm_int__, c_void_p,
+                                      __spm_int__ ]
+    libspm.spmScalVector( flt, alpha, n, x, incx )
+
+def pyspm_spmSort( spm ):
+    libspm.spmSort.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmSort.restype = c_int
+    return libspm.spmSort( spm )
+
+def pyspm_spmMergeDuplicate( spm ):
+    libspm.spmMergeDuplicate.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmMergeDuplicate.restype = __spm_int__
+    return libspm.spmMergeDuplicate( spm )
+
+def pyspm_spmSymmetrize( spm ):
+    libspm.spmSymmetrize.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmSymmetrize.restype = __spm_int__
+    return libspm.spmSymmetrize( spm )
+
+def pyspm_spmCheckAndCorrect( spm ):
+    libspm.spmCheckAndCorrect.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmCheckAndCorrect.restype = POINTER(pyspm_spmatrix_t)
+    return libspm.spmCheckAndCorrect( spm )
+
+def pyspm_spmGenRHS( type, nrhs, spm, x, ldx, b, ldb ):
+    libspm.spmGenRHS.argtypes = [ c_int, __spm_int__, POINTER(pyspm_spmatrix_t),
+                                  c_void_p, __spm_int__, c_void_p, __spm_int__ ]
+    libspm.spmGenRHS.restype = c_int
+    return libspm.spmGenRHS( type, nrhs, spm, x, ldx, b, ldb )
+
+def pyspm_spmCheckAxb( eps, nrhs, spm, x0, ldx0, b, ldb, x, ldx ):
+    libspm.spmCheckAxb.argtypes = [ c_double, __spm_int__,
+                                    POINTER(pyspm_spmatrix_t), c_void_p,
+                                    __spm_int__, c_void_p, __spm_int__,
+                                    c_void_p, __spm_int__ ]
+    libspm.spmCheckAxb.restype = c_int
+    return libspm.spmCheckAxb( eps, nrhs, spm, x0, ldx0, b, ldb, x, ldx )
+
+def pyspm_spmIntConvert( n, input ):
+    libspm.spmIntConvert.argtypes = [ __spm_int__, c_int_p ]
+    libspm.spmIntConvert.restype = POINTER(__spm_int__)
+    return libspm.spmIntConvert( n, input )
+
+def pyspm_spmLoad( spm ):
+    libspm.spmLoad.argtypes = [ POINTER(pyspm_spmatrix_t), c_void_p ]
+    libspm.spmLoad.restype = c_int
+    return libspm.spmLoad( spm, None )
+
+def pyspm_spmSave( spm ):
+    libspm.spmSave.argtypes = [ POINTER(pyspm_spmatrix_t), c_void_p ]
+    libspm.spmSave.restype = c_int
+    return libspm.spmSave( spm, None )
+
+def pyspm_spmReadDriver( driver, filename, spm ):
+    libspm.spmReadDriver.argtypes = [ c_int, c_char_p,
+                                      POINTER(pyspm_spmatrix_t) ]
+    libspm.spmReadDriver.restype = c_int
+    return libspm.spmReadDriver( driver, filename, spm )
+
+def pyspm_spmParseLaplacianInfo( filename, flttype, dim1, dim2, dim3, alpha,
+                                 beta ):
+    libspm.spmParseLaplacianInfo.argtypes = [ c_char_p, c_int_p,
+                                              POINTER(__spm_int__),
+                                              POINTER(__spm_int__),
+                                              POINTER(__spm_int__),
+                                              POINTER(c_double),
+                                              POINTER(c_double) ]
+    libspm.spmParseLaplacianInfo.restype = c_int
+    return libspm.spmParseLaplacianInfo( filename, flttype, dim1, dim2, dim3,
+                                         alpha, beta )
+
+def pyspm_spm2Dense( spm ):
+    libspm.spm2Dense.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spm2Dense( spm )
+
+def pyspm_spmPrint( spm ):
+    libspm.spmPrint.argtypes = [ POINTER(pyspm_spmatrix_t), c_void_p ]
+    libspm.spmPrint( spm, None )
+
+def pyspm_spmPrintInfo( spm ):
+    libspm.spmPrintInfo.argtypes = [ POINTER(pyspm_spmatrix_t), c_void_p ]
+    libspm.spmPrintInfo( spm, None )
+
+def pyspm_spmExpand( spm ):
+    libspm.spmExpand.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmExpand.restype = POINTER(pyspm_spmatrix_t)
+    return libspm.spmExpand( spm )
+
+def pyspm_spmDofExtend( spm, type, dof ):
+    libspm.spmDofExtend.argtypes = [ POINTER(pyspm_spmatrix_t), c_int, c_int ]
+    libspm.spmDofExtend.restype = POINTER(pyspm_spmatrix_t)
+    return libspm.spmDofExtend( spm, type, dof )
+
+
diff --git a/wrappers/python/spm/enum.py.in b/wrappers/python/spm/enum.py.in
new file mode 100644
index 0000000000000000000000000000000000000000..b73be1c74bf3eb52e0dfaa79f0588d98fa12d44b
--- /dev/null
+++ b/wrappers/python/spm/enum.py.in
@@ -0,0 +1,155 @@
+"""
+
+ @file enum.py
+
+ SPM python wrapper to define enums and datatypes
+
+ @copyright 2017-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @author Louis Poirel
+ @date 2017-05-04
+
+This file has been automatically generated with gen_wrappers.py
+
+"""
+from ctypes import *
+import numpy as np
+
+# Start with __ to prevent broadcast to file importing enum
+__spm_int__ = @SPM_PYTHON_INTEGER@
+
+class verbose:
+    Not = 0
+    No  = 1
+    Yes = 2
+
+class coeftype:
+    Pattern   = 0
+    Float     = 2
+    Double    = 3
+    Complex32 = 4
+    Complex64 = 5
+
+    @staticmethod
+    def getptype ( dtype ):
+        np_dict = {
+            np.dtype('float32')    : coeftype.Float,
+            np.dtype('float64')    : coeftype.Double,
+            np.dtype('complex64')  : coeftype.Complex32,
+            np.dtype('complex128') : coeftype.Complex64,
+        }
+        if dtype in np_dict:
+            return np_dict[dtype]
+        else:
+            return -1
+
+    @staticmethod
+    def getctype ( flttype ):
+        class c_float_complex(Structure):
+            _fields_ = [("r",c_float),("i", c_float)]
+        class c_double_complex(Structure):
+            _fields_ = [("r",c_double),("i", c_double)]
+
+        np_dict = {
+            coeftype.Float     : c_float,
+            coeftype.Double    : c_double,
+            coeftype.Complex32 : c_float_complex,
+            coeftype.Complex64 : c_double_complex
+        }
+        if flttype in np_dict:
+            return np_dict[flttype]
+        else:
+            return -1
+
+    @staticmethod
+    def getnptype ( flttype ):
+        np_dict = {
+            coeftype.Float     : np.dtype('float32'),
+            coeftype.Double    : np.dtype('float64'),
+            coeftype.Complex32 : np.dtype('complex64'),
+            coeftype.Complex64 : np.dtype('complex128')
+        }
+        if flttype in np_dict:
+            return np_dict[flttype]
+        else:
+            return -1
+
+class fmttype:
+    CSC = 0
+    CSR = 1
+    IJV = 2
+
+class error:
+    SUCCESS        = 0
+    UNKNOWN        = 1
+    ALLOC          = 2
+    NOTIMPLEMENTED = 3
+    OUTOFMEMORY    = 4
+    THREAD         = 5
+    INTERNAL       = 6
+    BADPARAMETER   = 7
+    FILE           = 8
+    INTEGER_TYPE   = 9
+    IO             = 10
+    MPI            = 11
+
+class driver:
+    RSA        = 0
+    HB         = 1
+    IJV        = 2
+    MM         = 3
+    Laplacian  = 4
+    XLaplacian = 5
+    Graph      = 6
+    SPM        = 7
+
+class rhstype:
+    One  = 0
+    I    = 1
+    RndX = 2
+    RndB = 3
+
+class layout:
+    RowMajor = 101
+    ColMajor = 102
+
+class trans:
+    NoTrans   = 111
+    Trans     = 112
+    ConjTrans = 113
+
+class mtxtype:
+    General   = trans.NoTrans
+    Symmetric = trans.Trans
+    Hermitian = trans.ConjTrans
+    SymPosDef = trans.ConjTrans + 1
+    HerPosDef = trans.ConjTrans + 2
+
+class uplo:
+    Upper      = 121
+    Lower      = 122
+    UpperLower = 123
+
+class diag:
+    NonUnit = 131
+    Unit    = 132
+
+class side:
+    Left  = 141
+    Right = 142
+
+class normtype:
+    One       = 171
+    Frobenius = 174
+    Inf       = 175
+    Max       = 177
+
+class dir:
+    Forward  = 391
+    Backward = 392
+
+
diff --git a/wrappers/python/spm/spm.py b/wrappers/python/spm/spm.py
new file mode 100644
index 0000000000000000000000000000000000000000..875d9763c35598520fdfcdcb2725eeae07641fb6
--- /dev/null
+++ b/wrappers/python/spm/spm.py
@@ -0,0 +1,227 @@
+"""
+ @file spm.py
+
+ SPM python wrapper
+
+ @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @date 2017-05-04
+
+"""
+from ctypes import *
+import numpy as np
+
+__use_sps__ = True
+try:
+    import scipy.sparse as sps
+except ImportError:
+    __use_sps__ = False
+
+from .enum    import *
+from .enum    import __spm_int__
+from .__spm__ import *
+
+class spmatrix():
+
+    dtype = None
+
+    def __init__( self, A=None, mtxtype_=mtxtype.General, driver=None, filename="" ):
+        """
+        Initialize the SPM wrapper by loading the libraries
+        """
+        if mtxtype_ == mtxtype.SymPosDef:
+            mtxtype_ = mtxtype.Symmetric
+        if mtxtype_ == mtxtype.HerPosDef:
+            mtxtype_ = mtxtype.Hermitian
+
+        self.spm_c = pyspm_spmatrix_t( mtxtype_,
+                                     coeftype.Double,
+                                     fmttype.CSC,
+                                     0, 0, 0, 0, 0, 0, 0, 0,
+                                     1, None,
+                                     layout.ColMajor,
+                                     None, None, None, None )
+        self.id_ptr = pointer( self.spm_c )
+        if A is not None:
+            self.fromsps( A, mtxtype_ )
+        elif driver is not None:
+            self.fromdriver( driver, filename )
+
+    if __use_sps__:
+        def fromsps( self, A, mtxtype_=mtxtype.General ):
+            """
+            Initialize the SPM wrapper by loading the libraries
+            """
+
+            if not sps.isspmatrix(A):
+                raise TypeError( "A must be of type scipy.sparse matrix" )
+
+            # Assume A is already in Scipy sparse format
+            self.dtype = A.dtype
+            flttype = coeftype.getptype( A.dtype )
+            if flttype == -1:
+                raise TypeError( "Invalid data type. Must be part of (f4, f8, c8 or c16)" )
+
+            A = sps.csc_matrix( A )
+            A.sort_indices()
+
+            # Pointer variables
+            self.py_colptr = np.array( A.indptr[:],  dtype=__spm_int__ )
+            self.py_rowptr = np.array( A.indices[:], dtype=__spm_int__ )
+            self.py_values = np.array( A.data[:] )
+
+            if mtxtype_ == mtxtype.SymPosDef:
+                mtxtype_ = mtxtype.Symmetric
+            if mtxtype_ == mtxtype.HerPosDef:
+                mtxtype_ = mtxtype.Hermitian
+
+            self.spm_c.mtxtype  = mtxtype_
+            self.spm_c.flttype  = flttype
+            self.spm_c.fmttype  = fmttype.CSC
+            self.spm_c.n        = A.shape[0]
+            self.spm_c.nnz      = A.getnnz()
+            self.spm_c.dof      = 1
+            self.spm_c.dofs     = None
+            self.spm_c.layout   = layout.ColMajor
+            self.spm_c.colptr   = self.py_colptr.ctypes.data_as(POINTER(__spm_int__))
+            self.spm_c.rowptr   = self.py_rowptr.ctypes.data_as(POINTER(__spm_int__))
+            self.spm_c.loc2glob = None
+            self.spm_c.values   = self.py_values.ctypes.data_as(c_void_p)
+
+            self.id_ptr = pointer( self.spm_c )
+
+            self.updateComputedField()
+            self.checkAndCorrect()
+
+        def tosps( self ):
+            """
+            Return a Scipy sparse matrix
+            """
+            n      = int( self.spm_c.n )
+            nnz    = int( self.spm_c.nnz )
+            cflt   = coeftype.getctype(  self.spm_c.flttype )
+            nflt   = coeftype.getnptype( self.spm_c.flttype )
+            colptr = self.py_colptr.copy()
+            rowptr = self.py_rowptr.copy()
+            values = self.py_values.copy()
+
+            baseval = colptr[0]
+            colptr = colptr - baseval
+            rowptr = rowptr - baseval
+
+            return sps.csc_matrix((values, rowptr, colptr), shape=(n, n))
+
+    def fromdriver( self, driver=driver.Laplacian, filename="10:10:10" ):
+        """
+        Initialize the SPM wrapper by loading the libraries
+        """
+        if filename == "":
+            raise ValueError("filename must be prodived")
+
+        pyspm_spmReadDriver( driver, filename.encode('utf-8'), self.id_ptr )
+
+        self.dtype = coeftype.getnptype( self.spm_c.flttype )
+        self.checkAndCorrect()
+
+    def scale( self, alpha ):
+        return pyspm_spmScalMatrix( float(alpha), self.id_ptr )
+
+    def norm( self, ntype=normtype.Frobenius ):
+        return pyspm_spmNorm( ntype, self.id_ptr )
+
+    def base( self, baseval ):
+        pyspm_spmBase( self.id_ptr, baseval )
+
+    def findBase( self ):
+        return pyspm_spmFindBase( self.id_ptr )
+
+    def updateComputedField( self ):
+        pyspm_spmUpdateComputedFields( self.id_ptr )
+
+    def printInfo( self ):
+        pyspm_spmPrintInfo( self.id_ptr )
+
+    def printSpm( self ):
+        pyspm_spmPrint( self.id_ptr )
+
+    def checkAndCorrect( self ):
+        spm1 = self.id_ptr
+        spm2 = pyspm_spmCheckAndCorrect( self.id_ptr )
+        if (( spm1.contents.fmttype == spm2.contents.fmttype ) and
+            ( spm1.contents.nnzexp  == spm2.contents.nnzexp  ) ):
+            return
+        self.spm_c = cast( spm2, POINTER(pyspm_spmatrix_t) ).contents
+        self.id_ptr = pointer( self.spm_c )
+
+        self.py_colptr = np.frombuffer( (__spm_int__ * (n+1)).from_address( cast(self.spm_c.colptr, c_void_p).value ), __spm_int__ ).copy()
+        self.py_rowptr = np.frombuffer( (__spm_int__ *  nnz ).from_address( cast(self.spm_c.rowptr, c_void_p).value ), __spm_int__ ).copy()
+        self.py_values = np.frombuffer( (cflt       *  nnz ).from_address( self.spm_c.values ), nflt ).copy()
+
+    def __checkVector( self, n, nrhs, x ):
+        if x.dtype != self.dtype:
+            raise TypeError( "Vectors must use the same arithmetic as the spm" )
+
+        if x.ndim > 2:
+            raise TypeError( "Vectors must be of dimension 1 or 2" )
+
+        if x.shape[0] < n:
+            raise TypeError( "Vectors must be of size at least ", n )
+
+        if (x.ndim == 1 and nrhs > 1) or (x.ndim>1 and x.shape[1] < nrhs):
+            raise TypeError( "At least nrhs vectors must be stored in the vector" )
+
+    def checkAxb( self, x0, b, x, nrhs=-1 ):
+        # if libspm == None:
+        #     raise EnvironmentError( "SPM Instance badly instanciated" )
+
+        b = np.array(b, self.dtype)
+        x = np.asarray(x, self.dtype)
+
+        n = self.spm_c.n
+        if nrhs == -1:
+            if x.ndim == 1:
+                nrhs = 1
+            else:
+                nrhs = x.shape[1]
+
+        if x0 is not None:
+            x0 = np.asarray(x0, self.dtype)
+            self.__checkVector( n, nrhs, x0 )
+            ldx0 = x0.shape[0]
+            x0ptr = x0.ctypes.data_as(c_void_p)
+        else:
+            ldx0 = 1
+            x0ptr = None
+        self.__checkVector( n, nrhs, b )
+        self.__checkVector( n, nrhs, x )
+
+        ldb  = b.shape[0]
+        ldx  = x.shape[0]
+
+        pyspm_spmCheckAxb( -1., nrhs, self.id_ptr,
+                           x0ptr, ldx0,
+                           b.ctypes.data_as(c_void_p), ldb,
+                           x.ctypes.data_as(c_void_p), ldx )
+
+    def genRHS( self, rhstype=rhstype.One, nrhs=1 ):
+        # if libspm == None:
+        #     raise EnvironmentError( "SPM Instance badly instanciated" )
+
+        n = self.spm_c.n
+        b = np.zeros((n, nrhs), self.dtype)
+        x = np.zeros((n, nrhs), self.dtype)
+
+        ldb = b.shape[0]
+        ldx = x.shape[0]
+
+        self.__checkVector( n, nrhs, x )
+        self.__checkVector( n, nrhs, b )
+
+        pyspm_spmGenRHS( rhstype, nrhs, self.id_ptr,
+                         x.ctypes.data_as(c_void_p), ldx,
+                         b.ctypes.data_as(c_void_p), ldb )
+        return x, b
diff --git a/wrappers/python/spm_driver.py b/wrappers/python/spm_driver.py
new file mode 100755
index 0000000000000000000000000000000000000000..fa88e594052daa326843ce7c7741aa2288392d64
--- /dev/null
+++ b/wrappers/python/spm_driver.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+"""
+ @file spm_driver.py
+
+ SPM example to generate a sparse matrix from the spm drivers
+
+ @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @author Louis Poirel
+ @date 2017-05-04
+
+"""
+import spm
+import numpy as np
+
+# Hack to make sure that the mkl is loaded
+tmp = np.eye(2).dot(np.ones(2))
+
+# Load a sparse matrix from the Laplacian driver
+A = spm.spmatrix( None, driver=spm.driver.Laplacian, filename="10:10:10:2.:1." )
+
+# Example from a RSA file
+#A = spm( None, driver=driver.RSA, filename="$PASTIX_DIR/test/matrix/oilpan.rsa" )
+
+A.printInfo()
+
+# Scale A for low-rank: A / ||A||_f
+norm = A.norm()
+A.scale( 1. / norm )
+
+# Generate b and x0 vectors such that A * x0 = b
+nrhs = 10
+x0, b = A.genRHS( spm.rhstype.RndX, nrhs )
+
+# Check that A * x = b
+x = x0.copy()
+A.checkAxb( x0, b, x )
+
diff --git a/wrappers/python/spm_scipy.py b/wrappers/python/spm_scipy.py
new file mode 100755
index 0000000000000000000000000000000000000000..6a574b6a5100de854189502139aac8db74e3fb97
--- /dev/null
+++ b/wrappers/python/spm_scipy.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+"""
+ @file spm_scipy.py
+
+ SPM example to geneate a sparse matrix from Scipy to SPM
+
+ @copyright 2017-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+                      Univ. Bordeaux. All rights reserved.
+
+ @version 6.0.0
+ @author Pierre Ramet
+ @author Mathieu Faverge
+ @author Louis Poirel
+ @date 2017-05-04
+
+"""
+import spm
+import scipy.sparse as sps
+import numpy as np
+
+# Hack to make sure that the mkl is loaded
+tmp = np.eye(2).dot(np.ones(2))
+
+# Set the problem
+n = 9
+A = sps.spdiags([np.ones(n)*i for i in [4, -1, -1, -1, -1]],
+                [0, 1, 3, -1, -3], n, n)
+x0 = np.arange(n)
+b  = A.dot(x0)
+
+spmA = spm.spmatrix( A )
+
+# Check that A * x = b
+x = x0.copy()
+spmA.checkAxb( x0, b, x )
+