diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f25b1b2b92f3821a5f2f534a70e6e2ef2ad79d7d..acde4028783e462b1f6a2adb1f75554e616bdbf9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ stages:
   - build
   - test
   - sonar
-#  - doc
+#  - deploy
 
 build_spm:
   stage: build
@@ -16,7 +16,7 @@ build_spm:
     - git submodule update --init --recursive
     - mkdir build
     - cd build
-    - cmake .. -DCMAKE_INSTALL_PREFIX=${PWD}/../install -DBUILD_SHARED_LIBS=ON -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_C_FLAGS="-O0 -g -fPIC --coverage -Wall -fdiagnostics-show-option -fno-inline" -DCMAKE_EXE_LINKER_FLAGS="--coverage" -DPASTIX_INT64=OFF
+    - cmake .. -DCMAKE_INSTALL_PREFIX=${PWD}/../install -DBUILD_SHARED_LIBS=ON -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_C_FLAGS="-O0 -g -fPIC --coverage -Wall -fdiagnostics-show-option -fno-inline" -DCMAKE_EXE_LINKER_FLAGS="--coverage" -DSPM_INT64=OFF -DCMAKE_EXPORT_COMPILE_COMMANDS=ON
     - make -j 4 | tee ../pastix-build.log
     - make install | tee -a ../pastix-build.log
   only:
@@ -34,6 +34,7 @@ test_spm:
       - spm.lcov
       - spm-gcov.log
   script:
+    - git submodule update --init --recursive
     - source install/bin/spm_env.sh
     - (cd build &&
        eval "ctest
@@ -68,9 +69,12 @@ sonar_spm:
     - master@solverstack/spm
 
 # pages:
-#   stage: doc
+#   stage: deploy
 #   dependencies:
 #     - build_spm
+#   artifacts:
+#     paths:
+#     - public
 #   script:
 #     - git submodule update --init --recursive
 #     - mkdir -p build
@@ -78,9 +82,6 @@ sonar_spm:
 #     - cmake .. -DBUILD_DOCUMENTATION=ON -DSPM_INT64=OFF
 #     - make docs
 #     - mv docs/out/html ../public/
-#   artifacts:
-#     paths:
-#       - public
 #   only:
 #     - master@solverstack/spm
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ee61de1085c7a5456ec54ecebc1f66c229da3a17..8c904a5fa0fc36c17c57281b462e92b507c2c991 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -249,9 +249,17 @@ install(FILES
   DESTINATION include )
 
 ### Build pkg-config and environment file
-include(GenSPMPkgConfig)
-spm_generate_pkgconfig_file()
-spm_generate_env_file()
+include(GenPkgConfig)
+list(APPEND SPM_PKGCONFIG_LIBS_PRIVATE
+  ${LAPACKE_LIBRARIES_DEP}
+  ${CBLAS_LIBRARIES_DEP}
+  )
+generate_pkgconfig_files(
+  "${CMAKE_CURRENT_SOURCE_DIR}/tools/spm.pc.in"
+  "${CMAKE_CURRENT_SOURCE_DIR}/tools/spmf.pc.in"
+  PROJECTNAME SPM )
+
+generate_env_file( PROJECTNAME SPM )
 
 ### Add documented files to the global property
 add_documented_files(
diff --git a/cmake_modules/GenSPMPkgConfig.cmake b/cmake_modules/GenSPMPkgConfig.cmake
deleted file mode 100644
index 935b432f3b3df4aa1eb9fd34509c4c162e815cfb..0000000000000000000000000000000000000000
--- a/cmake_modules/GenSPMPkgConfig.cmake
+++ /dev/null
@@ -1,125 +0,0 @@
-###
-#
-# @copyright (c) 2009-2014 The University of Tennessee and The University
-#                          of Tennessee Research Foundation.
-#                          All rights reserved.
-# @copyright (c) 2012-2017 Inria. All rights reserved.
-# @copyright (c) 2012-2014 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria, Univ. Bordeaux. All rights reserved.
-#
-###
-#
-#  @file GenPkgConfig.cmake
-#
-#  @project MORSE
-#  MORSE is a software package provided by:
-#     Inria Bordeaux - Sud-Ouest,
-#     Univ. of Tennessee,
-#     King Abdullah Univesity of Science and Technology
-#     Univ. of California Berkeley,
-#     Univ. of Colorado Denver.
-#
-#  @version 0.9.1
-#  @author Cedric Castagnede
-#  @author Emmanuel Agullo
-#  @author Mathieu Faverge
-#  @author Florent Pruvost
-#  @date 10-11-2014
-#
-###
-
-###
-#
-# CONVERT_LIBSTYLE_TO_PKGCONFIG: convert a libraries list to follow the pkg-config style
-#                                used in CLEAN_LIB_LIST
-#
-###
-# macro(CONVERT_LIBSTYLE_TO_PKGCONFIG _liblist)
-#     set(${_liblist}_CPY "${${_liblist}}")
-#     set(${_liblist} "")
-#     foreach(_dep ${${_liblist}_CPY})
-#         if (${_dep} MATCHES "^/")
-#             get_filename_component(dep_libname ${_dep} NAME)
-#             get_filename_component(dep_libdir  ${_dep} DIRECTORY)
-#             string(REPLACE "lib"    "" dep_libname "${dep_libname}")
-#             string(REPLACE ".so"    "" dep_libname "${dep_libname}")
-#             string(REPLACE ".a"     "" dep_libname "${dep_libname}")
-#             string(REPLACE ".dylib" "" dep_libname "${dep_libname}")
-#             string(REPLACE ".dll"   "" dep_libname "${dep_libname}")
-#             list(APPEND ${_liblist} -L${dep_libdir} -l${dep_libname})
-#         elseif(NOT ${_dep} MATCHES "^-")
-#             list(APPEND ${_liblist} "-l${_dep}")
-#         else()
-#             list(APPEND ${_liblist} ${_dep})
-#         endif()
-#     endforeach()
-# endmacro(CONVERT_LIBSTYLE_TO_PKGCONFIG)
-
-###
-#
-# CLEAN_LIB_LIST: clean libraries lists to follow the pkg-config style
-#                 used in GENERATE_PKGCONFIG_FILE
-#
-###
-#macro(CLEAN_LIB_LIST _package)
-#    list(REMOVE_DUPLICATES ${_package}_PKGCONFIG_LIBS)
-#    list(REMOVE_DUPLICATES ${_package}_PKGCONFIG_LIBS_PRIVATE)
-#    list(REMOVE_DUPLICATES ${_package}_PKGCONFIG_REQUIRED)
-#    list(REMOVE_DUPLICATES ${_package}_PKGCONFIG_REQUIRED_PRIVATE)
-#    convert_libstyle_to_pkgconfig(${_package}_PKGCONFIG_LIBS)
-#    convert_libstyle_to_pkgconfig(${_package}_PKGCONFIG_LIBS_PRIVATE)
-#    string(REPLACE ";" " " ${_package}_PKGCONFIG_LIBS "${${_package}_PKGCONFIG_LIBS}")
-#    string(REPLACE ";" " " ${_package}_PKGCONFIG_LIBS_PRIVATE "${${_package}_PKGCONFIG_LIBS_PRIVATE}")
-#    string(REPLACE ";" " " ${_package}_PKGCONFIG_REQUIRED "${${_package}_PKGCONFIG_REQUIRED}")
-#    string(REPLACE ";" " " ${_package}_PKGCONFIG_REQUIRED_PRIVATE "${${_package}_PKGCONFIG_REQUIRED_PRIVATE}")
-#endmacro(CLEAN_LIB_LIST)
-
-###
-#
-# GENERATE_PKGCONFIG_FILE: generate files spm.pc
-#
-###
-macro(spm_generate_pkgconfig_file)
-
-  set(SPM_PKGCONFIG_LIBS "-lspm")
-  set(SPM_PKGCONFIG_LIBS_PRIVATE "-lm")
-  set(SPM_PKGCONFIG_REQUIRED "")
-  set(SPM_PKGCONFIG_REQUIRED_PRIVATE "")
-
-  #clean_lib_list(SPM)
-
-  set(_output_spm_file "${CMAKE_BINARY_DIR}/spm.pc")
-  configure_file(
-    "${CMAKE_CURRENT_SOURCE_DIR}/tools/spm.pc.in"
-    "${_output_spm_file}"
-    @ONLY
-    )
-  install(
-    FILES ${_output_spm_file}
-    DESTINATION lib/pkgconfig
-    )
-
-endmacro(spm_generate_pkgconfig_file)
-
-###
-#
-# generate_env_file: generate files pastix.pc
-#
-###
-macro(spm_generate_env_file)
-
-    # Create .sh file
-    # ---------------
-    configure_file(
-      "${CMAKE_CURRENT_SOURCE_DIR}/tools/spm_env.sh.in"
-      "${CMAKE_BINARY_DIR}/bin/spm_env.sh" @ONLY)
-
-    # installation
-    # ------------
-    install(FILES "${CMAKE_BINARY_DIR}/bin/spm_env.sh"
-      DESTINATION bin)
-
-endmacro(spm_generate_env_file)
-
-##
-## @end file GenPkgConfig.cmake
-##
diff --git a/cmake_modules/morse_cmake b/cmake_modules/morse_cmake
index 47b8241b579347ecc7213ccb7eba0b387f036c98..d44efd0e4af75cf4be6b9483198c595b03ec083a 160000
--- a/cmake_modules/morse_cmake
+++ b/cmake_modules/morse_cmake
@@ -1 +1 @@
-Subproject commit 47b8241b579347ecc7213ccb7eba0b387f036c98
+Subproject commit d44efd0e4af75cf4be6b9483198c595b03ec083a
diff --git a/include/cblas.h b/include/cblas.h
index 14756f24ea4c49cf4cff4a070ad355e09c8baa54..fa70cffa23ecddbc76b2ce0f694dffc6d8f15720 100644
--- a/include/cblas.h
+++ b/include/cblas.h
@@ -16,11 +16,11 @@ extern "C" {            /* Assume C declarations for C++ */
     #define CBLAS_INDEX int
 #endif
 
-typedef enum {CblasRowMajor=101, CblasColMajor=102} CBLAS_LAYOUT;
-typedef enum {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
-typedef enum {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
-typedef enum {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
-typedef enum {CblasLeft=141, CblasRight=142} CBLAS_SIDE;
+typedef enum CBLAS_LAYOUT {CblasRowMajor=101, CblasColMajor=102} CBLAS_LAYOUT;
+typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
+typedef enum CBLAS_UPLO {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
+typedef enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
+typedef enum CBLAS_SIDE {CblasLeft=141, CblasRight=142} CBLAS_SIDE;
 
 typedef CBLAS_LAYOUT CBLAS_ORDER; /* this for backward compatibility with CBLAS_ORDER */
 
diff --git a/include/spm.h b/include/spm.h
index 6fb956e2ae241cf076c3e3c527956cbed4605f18..9d4ac2f90958503c5f63129762f8891ba4e3de3c 100644
--- a/include/spm.h
+++ b/include/spm.h
@@ -81,8 +81,9 @@ typedef struct spmatrix_s {
  * @name SPM basic subroutines
  * @{
  */
-void spmInit( spmatrix_t *spm );
-void spmExit( spmatrix_t *spm );
+void spmInit ( spmatrix_t *spm );
+void spmAlloc( spmatrix_t *spm );
+void spmExit ( spmatrix_t *spm );
 
 spmatrix_t *spmCopy( const spmatrix_t *spm );
 void        spmBase( spmatrix_t *spm, int baseval );
@@ -110,10 +111,10 @@ void   spmScalVector( spm_coeftype_t flt, double alpha, spm_int_t n, void *x, sp
  * @name SPM subroutines to check format
  * @{
  */
-int         spmSort( spmatrix_t *spm );
-spm_int_t   spmMergeDuplicate( spmatrix_t *spm );
-spm_int_t   spmSymmetrize( spmatrix_t *spm );
-spmatrix_t *spmCheckAndCorrect( spmatrix_t *spm );
+int       spmSort( spmatrix_t *spm );
+spm_int_t spmMergeDuplicate( spmatrix_t *spm );
+spm_int_t spmSymmetrize( spmatrix_t *spm );
+int       spmCheckAndCorrect( const spmatrix_t *spm_in, spmatrix_t *spm_out );
 
 /**
  * @}
diff --git a/src/drivers/skitf.f b/src/drivers/skitf.f
index 5c5973459219657a2d8fcec4fc0306eaf972155c..cf609894d89de78f67652c6339f09805d7d9c951 100644
--- a/src/drivers/skitf.f
+++ b/src/drivers/skitf.f
@@ -29,16 +29,16 @@ c         for the rows: perm(i) is the destination of row i in the
 c         permuted matrix -- also the destination of column i in case
 c         permutation is symmetric (job .le. 2)
 c
-c qperm	= same thing for the columns. This should be provided only
+c qperm = same thing for the columns. This should be provided only
 c         if job=3 or job=4, i.e., only in the case of a nonsymmetric
 c         permutation of rows and columns. Otherwise qperm is a dummy
 c
-c job	= integer indicating the work to be done:
+c job = integer indicating the work to be done:
 c * job = 1,2 permutation is symmetric  Ao :== P * A * transp(P)
-c               job = 1	permute a, ja, ia into ao, jao, iao
+c               job = 1 permute a, ja, ia into ao, jao, iao
 c               job = 2 permute matrix ignoring real values.
 c * job = 3,4 permutation is non-symmetric  Ao :== P * A * Q
-c               job = 3	permute a, ja, ia into ao, jao, iao
+c               job = 3 permute a, ja, ia into ao, jao, iao
 c               job = 4 permute matrix ignoring real values.
 c
 c on return:
@@ -168,7 +168,7 @@ c perm  = integer array of length n containing the permutation arrays
 c         for the rows: perm(i) is the destination of row i in the
 c         permuted matrix -- also the destination of column i.
 c
-c rperm	= reverse permutation. defined by rperm(iperm(i))=i for all i
+c rperm = reverse permutation. defined by rperm(iperm(i))=i for all i
 c
 c job   = job indicator. if (job .ne.1) values are not copied (i.e.,
 c         only pattern is copied).
@@ -351,17 +351,17 @@ c           ao, jao, iao in csr format
 c
 c a, ja, ia = input matrix in csr format
 c
-c iperm	= integer array of length nrow containing the reverse permutation
+c iperm = integer array of length nrow containing the reverse permutation
 c         array for the rows. row number iperm(j) in permuted matrix PA
 c         used to be row number j in unpermuted matrix.
 c         ---> a(i,j) in the permuted matrix was a(iperm(i),j)
 c         in the inout matrix.
 c
-c job	= integer indicating the work to be done:
+c job = integer indicating the work to be done:
 c               job .ne. 1 : get structure only of output matrix,,
 c               i.e., ignore real values. (in which case arrays a
 c               and ao are not used nor accessed).
-c               job = 1	get complete data structure of output matrix.
+c               job = 1 get complete data structure of output matrix.
 c               (i.e., including arrays ao and iao).
 c------------
 c on return:
@@ -809,128 +809,15 @@ c
 c---------end-of-readmt-------------------------------------------------
 c-----------------------------------------------------------------------
       end
-c-----------------------------------------------------------------------
-      subroutine roscal(nrow,job,nrm,a,ja,ia,diag,b,jb,ib,ierr)
-      real*8 a(*), b(*), diag(nrow)
-      integer nrow,job,nrm,ja(*),jb(*),ia(nrow+1),ib(nrow+1),ierr
-c-----------------------------------------------------------------------
-c scales the rows of A such that their norms are one on return
-c 3 choices of norms: 1-norm, 2-norm, max-norm.
-c-----------------------------------------------------------------------
-c on entry:
-c ---------
-c nrow	= integer. The row dimension of A
-c
-c job   = integer. job indicator. Job=0 means get array b only
-c         job = 1 means get b, and the integer arrays ib, jb.
-c
-c nrm   = integer. norm indicator. nrm = 1, means 1-norm, nrm =2
-c                  means the 2-nrm, nrm = 0 means max norm
-c
-c a,
-c ja,
-c ia   = Matrix A in compressed sparse row format.
-c
-c on return:
-c----------
-c
-c diag = diagonal matrix stored as a vector containing the matrix
-c        by which the rows have been scaled, i.e., on return
-c        we have B = Diag*A.
-c
-c b,
-c jb,
-c ib	= resulting matrix B in compressed sparse row sparse format.
-c
-c ierr  = error message. ierr=0     : Normal return
-c                        ierr=i > 0 : Row number i is a zero row.
-c Notes:
-c-------
-c 1)        The column dimension of A is not needed.
-c 2)        algorithm in place (B can take the place of A).
-c-----------------------------------------------------------------
-      call rnrms (nrow,nrm,a,ja,ia,diag)
-
-      ierr = 0
-      do 1 j=1, nrow
-         if (diag(j) .eq. 0.0d0) then
-            ierr = j
-            return
-         else
-            diag(j) = 1.0d0/diag(j)
-         endif
- 1    continue
-      call diamua(nrow,job,a,ja,ia,diag,b,jb,ib)
-      return
-c-------end-of-roscal---------------------------------------------------
-c-----------------------------------------------------------------------
-      end
-c-----------------------------------------------------------------------
-      subroutine coscal(nrow,job,nrm,a,ja,ia,diag,b,jb,ib,ierr)
-c-----------------------------------------------------------------------
-      real*8 a(*),b(*),diag(nrow)
-      integer nrow,job,ja(*),jb(*),ia(nrow+1),ib(nrow+1),ierr
-c-----------------------------------------------------------------------
-c scales the columns of A such that their norms are one on return
-c result matrix written on b, or overwritten on A.
-c 3 choices of norms: 1-norm, 2-norm, max-norm. in place.
-c-----------------------------------------------------------------------
-c on entry:
-c ---------
-c nrow	= integer. The row dimension of A
-c
-c job   = integer. job indicator. Job=0 means get array b only
-c         job = 1 means get b, and the integer arrays ib, jb.
-c
-c nrm   = integer. norm indicator. nrm = 1, means 1-norm, nrm =2
-c                  means the 2-nrm, nrm = 0 means max norm
-c
-c a,
-c ja,
-c ia   = Matrix A in compressed sparse row format.
-c
-c on return:
-c----------
-c
-c diag = diagonal matrix stored as a vector containing the matrix
-c        by which the columns have been scaled, i.e., on return
-c        we have B = A * Diag
-c
-c b,
-c jb,
-c ib	= resulting matrix B in compressed sparse row sparse format.
-c
-c ierr  = error message. ierr=0     : Normal return
-c                        ierr=i > 0 : Column number i is a zero row.
-c Notes:
-c-------
-c 1)        The column dimension of A is not needed.
-c 2)       algorithm in place (B can take the place of A).
-c-----------------------------------------------------------------
-      call cnrms (nrow,nrm,a,ja,ia,diag)
-      ierr = 0
-      do 1 j=1, nrow
-         if (diag(j) .eq. 0.0) then
-            ierr = j
-            return
-         else
-            diag(j) = 1.0d0/diag(j)
-         endif
- 1    continue
-      call amudia (nrow,job,a,ja,ia,diag,b,jb,ib)
-      return
-c--------end-of-coscal--------------------------------------------------
-c-----------------------------------------------------------------------
-      end
-      subroutine rnrms   (nrow, nrm, a, ja, ia, diag)
+      subroutine rnrms   (nrow, nrm, a, ia, diag)
       real*8 a(*), diag(nrow), scal
-      integer ja(*), ia(nrow+1)
+      integer ia(nrow+1)
 c-----------------------------------------------------------------------
 c gets the norms of each row of A. (choice of three norms)
 c-----------------------------------------------------------------------
 c on entry:
 c ---------
-c nrow	= integer. The row dimension of A
+c nrow = integer. The row dimension of A
 c
 c nrm   = integer. norm indicator. nrm = 1, means 1-norm, nrm =2
 c                  means the 2-nrm, nrm = 0 means max norm
@@ -981,7 +868,7 @@ c gets the norms of each column of A. (choice of three norms)
 c-----------------------------------------------------------------------
 c on entry:
 c ---------
-c nrow	= integer. The row dimension of A
+c nrow = integer. The row dimension of A
 c
 c nrm   = integer. norm indicator. nrm = 1, means 1-norm, nrm =2
 c                  means the 2-nrm, nrm = 0 means max norm
@@ -1030,7 +917,7 @@ c performs the matrix by matrix product B = Diag * A  (in place)
 c-----------------------------------------------------------------------
 c on entry:
 c ---------
-c nrow	= integer. The row dimension of A
+c nrow = integer. The row dimension of A
 c
 c job   = integer. job indicator. Job=0 means get array b only
 c         job = 1 means get b, and the integer arrays ib, jb.
@@ -1046,7 +933,7 @@ c----------
 c
 c b,
 c jb,
-c ib	= resulting matrix B in compressed sparse row sparse format.
+c ib = resulting matrix B in compressed sparse row sparse format.
 c
 c Notes:
 c-------
@@ -1087,7 +974,7 @@ c performs the matrix by matrix product B = A * Diag  (in place)
 c-----------------------------------------------------------------------
 c on entry:
 c ---------
-c nrow	= integer. The row dimension of A
+c nrow = integer. The row dimension of A
 c
 c job   = integer. job indicator. Job=0 means get array b only
 c         job = 1 means get b, and the integer arrays ib, jb.
@@ -1103,7 +990,7 @@ c----------
 c
 c b,
 c jb,
-c ib	= resulting matrix B in compressed sparse row sparse format.
+c ib = resulting matrix B in compressed sparse row sparse format.
 c
 c Notes:
 c-------
@@ -1144,7 +1031,7 @@ c performs the matrix sum  C = A+B.
 c-----------------------------------------------------------------------
 c on entry:
 c ---------
-c nrow	= integer. The row dimension of A and B
+c nrow = integer. The row dimension of A and B
 c ncol  = integer. The column dimension of A and B.
 c job   = integer. Job indicator. When job = 0, only the structure
 c                  (i.e. the arrays jc, ic) is computed and the
@@ -1156,9 +1043,9 @@ c ia   = Matrix A in compressed sparse row format.
 c
 c b,
 c jb,
-c ib	=  Matrix B in compressed sparse row format.
+c ib =  Matrix B in compressed sparse row format.
 c
-c nzmax	= integer. The  length of the arrays c and jc.
+c nzmax = integer. The  length of the arrays c and jc.
 c         amub will stop if the result matrix C  has a number
 c         of elements that exceeds exceeds nzmax. See ierr.
 c
@@ -1166,9 +1053,9 @@ c on return:
 c----------
 c c,
 c jc,
-c ic	= resulting matrix C in compressed sparse row sparse format.
+c ic = resulting matrix C in compressed sparse row sparse format.
 c
-c ierr	= integer. serving as error message.
+c ierr = integer. serving as error message.
 c         ierr = 0 means normal return,
 c         ierr .gt. 0 means that amub stopped while computing the
 c         i-th row  of C with i=ierr, because the number
@@ -1176,7 +1063,7 @@ c         of elements in C exceeds nzmax.
 c
 c work arrays:
 c------------
-c iw	= integer work array of length equal to the number of
+c iw = integer work array of length equal to the number of
 c         columns in A.
 c
 c-----------------------------------------------------------------------
@@ -1239,8 +1126,8 @@ c this subroutine transposes a matrix stored in a, ja, ia format.
 c ---------------
 c on entry:
 c----------
-c n	= dimension of A.
-c job	= integer to indicate whether to fill the values (job.eq.1) of the
+c n = dimension of A.
+c job = integer to indicate whether to fill the values (job.eq.1) of the
 c         matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1)
 c
 c ipos  = starting position in ao, jao of the transposed matrix.
@@ -1249,19 +1136,19 @@ c         Note: this may be useful if one needs to append the data structure
 c         of the transpose to that of A. In this case use for example
 c                call csrcsc (n,1,n+2,a,ja,ia,a,ja,ia(n+2))
 c         for any other normal usage, enter ipos=1.
-c a	= real array of length nnz (nnz=number of nonzero elements in input
+c a = real array of length nnz (nnz=number of nonzero elements in input
 c         matrix) containing the nonzero elements.
-c ja	= integer array of length nnz containing the column positions
+c ja = integer array of length nnz containing the column positions
 c         of the corresponding elements in a.
-c ia	= integer of size n+1. ia(k) contains the position in a, ja of
+c ia = integer of size n+1. ia(k) contains the position in a, ja of
 c         the beginning of the k-th row.
 c
 c on return:
 c ----------
 c output arguments:
-c ao	= real array of size nzz containing the "a" part of the transpose
-c jao	= integer array of size nnz containing the column indices.
-c iao	= integer array of size n+1 containing the "ia" index array of
+c ao = real array of size nzz containing the "a" part of the transpose
+c jao = integer array of size nnz containing the column indices.
+c iao = integer array of size n+1 containing the "ia" index array of
 c         the transpose.
 c
 c-----------------------------------------------------------------------
@@ -1284,9 +1171,9 @@ c this subroutine transposes a matrix stored in a, ja, ia format.
 c ---------------
 c on entry:
 c----------
-c n	= number of rows of CSR matrix.
+c n = number of rows of CSR matrix.
 c n2    = number of columns of CSC matrix.
-c job	= integer to indicate whether to fill the values (job.eq.1) of the
+c job = integer to indicate whether to fill the values (job.eq.1) of the
 c         matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1)
 c
 c ipos  = starting position in ao, jao of the transposed matrix.
@@ -1295,19 +1182,19 @@ c         Note: this may be useful if one needs to append the data structure
 c         of the transpose to that of A. In this case use for example
 c                call csrcsc2 (n,n,1,n+2,a,ja,ia,a,ja,ia(n+2))
 c         for any other normal usage, enter ipos=1.
-c a	= real array of length nnz (nnz=number of nonzero elements in input
+c a = real array of length nnz (nnz=number of nonzero elements in input
 c         matrix) containing the nonzero elements.
-c ja	= integer array of length nnz containing the column positions
+c ja = integer array of length nnz containing the column positions
 c         of the corresponding elements in a.
-c ia	= integer of size n+1. ia(k) contains the position in a, ja of
+c ia = integer of size n+1. ia(k) contains the position in a, ja of
 c         the beginning of the k-th row.
 c
 c on return:
 c ----------
 c output arguments:
-c ao	= real array of size nzz containing the "a" part of the transpose
-c jao	= integer array of size nnz containing the column indices.
-c iao	= integer array of size n+1 containing the "ia" index array of
+c ao = real array of size nzz containing the "a" part of the transpose
+c jao = integer array of size nnz containing the column indices.
+c iao = integer array of size n+1 containing the "ia" index array of
 c         the transpose.
 c
 c-----------------------------------------------------------------------
@@ -1808,18 +1695,18 @@ c this subroutine performs an in-place permutation of a real vector x
 c according to the permutation array perm(*), i.e., on return,
 c the vector x satisfies,
 c
-c	x(perm(j)) :== x(j), j=1,2,.., n
+c x(perm(j)) :== x(j), j=1,2,.., n
 c
 c-----------------------------------------------------------------------
 c on entry:
 c---------
 c n     = length of vector x.
 c perm  = integer array of length n containing the permutation  array.
-c x	= input vector
+c x = input vector
 c
 c on return:
 c----------
-c x	= vector x permuted according to x(perm(*)) :=  x(*)
+c x = vector x permuted according to x(perm(*)) :=  x(*)
 c
 c----------------------------------------------------------------------c
 c           Y. Saad, Sep. 21 1989                                      c
@@ -1828,7 +1715,7 @@ c local variables
       real*8 tmp, tmp1
 c
       init      = 1
-      tmp	= x(init)
+      tmp = x(init)
       ii        = perm(init)
       perm(init)= -perm(init)
       k         = 0
@@ -1860,8 +1747,8 @@ c
  65   init      = init+1
       if (init .gt. n) goto 101
       if (perm(init) .lt. 0) goto 65
-      tmp	= x(init)
-      ii	= perm(init)
+      tmp = x(init)
+      ii = perm(init)
       perm(init)=-perm(init)
       goto 6
 c
@@ -1883,18 +1770,18 @@ c this subroutine performs an in-place permutation of an integer vector
 c ix according to the permutation array perm(*), i.e., on return,
 c the vector x satisfies,
 c
-c	ix(perm(j)) :== ix(j), j=1,2,.., n
+c ix(perm(j)) :== ix(j), j=1,2,.., n
 c
 c-----------------------------------------------------------------------
 c on entry:
 c---------
 c n     = length of vector x.
 c perm  = integer array of length n containing the permutation  array.
-c ix	= input vector
+c ix = input vector
 c
 c on return:
 c----------
-c ix	= vector x permuted according to ix(perm(*)) :=  ix(*)
+c ix = vector x permuted according to ix(perm(*)) :=  ix(*)
 c
 c----------------------------------------------------------------------c
 c           Y. Saad, Sep. 21 1989                                      c
@@ -1903,7 +1790,7 @@ c local variables
       integer tmp, tmp1
 c
       init      = 1
-      tmp	= ix(init)
+      tmp = ix(init)
       ii        = perm(init)
       perm(init)= -perm(init)
       k         = 0
@@ -1935,8 +1822,8 @@ c
  65   init      = init+1
       if (init .gt. n) goto 101
       if (perm(init) .lt. 0) goto 65
-      tmp	= ix(init)
-      ii	= perm(init)
+      tmp = ix(init)
+      ii = perm(init)
       perm(init)=-perm(init)
       goto 6
 c
@@ -1972,8 +1859,8 @@ c         permuted matrix.
 c         ---> a(i,j) in the original matrix becomes a(perm(i),j)
 c         in the output  matrix.
 c
-c job	= integer indicating the work to be done:
-c               job = 1	permute a, ja, ia into ao, jao, iao
+c job = integer indicating the work to be done:
+c               job = 1 permute a, ja, ia into ao, jao, iao
 c                       (including the copying of real values ao and
 c                       the array iao).
 c               job .ne. 1 :  ignore real values.
@@ -2042,13 +1929,13 @@ c nrow  = row dimension of the matrix
 c
 c a, ja, ia = input matrix in csr format.
 c
-c perm	= integer array of length ncol (number of columns of A
+c perm = integer array of length ncol (number of columns of A
 c         containing the permutation array  the columns:
 c         a(i,j) in the original matrix becomes a(i,perm(j))
 c         in the output matrix.
 c
-c job	= integer indicating the work to be done:
-c               job = 1	permute a, ja, ia into ao, jao, iao
+c job = integer indicating the work to be done:
+c               job = 1 permute a, ja, ia into ao, jao, iao
 c                       (including the copying of real values ao and
 c                       the array iao).
 c               job .ne. 1 :  ignore real values ao and ignore iao.
@@ -2558,6 +2445,7 @@ c
 c     select domain with largest size
 c
       maxsiz = 0
+      nextdom = 0
       do j=1, idom
          if (size(j) .gt. maxsiz) then
             maxsiz = size(j)
@@ -2830,11 +2718,10 @@ c         print *,'mask',riord(j),j,nfirst
          mask(riord(j)) = 0
  12   continue
 c-----------------------------------------------------------------------
- 13   continue
 c
  1    nlev = nlev+1
       levels(nlev) = istart + 1
-      call add_lvst (istart,iend,nlev,riord,ja,ia,mask,maskval)
+      call add_lvst (istart,iend,riord,ja,ia,mask,maskval)
       if (istart .lt. iend) goto 1
  2    ii = ii+1
       if (ii .le. n) then
@@ -2854,7 +2741,7 @@ c
          endif
       endif
 c-----------------------------------------------------------------------
- 3    levels(nlev+1) = iend+1
+      levels(nlev+1) = iend+1
       do j=1, iend
          mask(riord(j)) = maskval
       enddo
@@ -2862,8 +2749,8 @@ c-----------------------------------------------------------------------
       return
       end
 c-----------------------------------------------------------------------
-      subroutine add_lvst(istart,iend,nlev,riord,ja,ia,mask,maskval)
-      integer nlev, nod, riord(*), ja(*), ia(*), mask(*)
+      subroutine add_lvst(istart,iend,riord,ja,ia,mask,maskval)
+      integer nod, riord(*), ja(*), ia(*), mask(*)
 c-------------------------------------------------------------
 c     adds one level set to the previous sets..
 c     span all nodes of previous mask
@@ -3120,7 +3007,7 @@ c
       nstuck = 0
 c-----------------------------------------------------------------------
  100  continue
-      idom = mindom(n,ndom,levst,link)
+      idom = mindom(n,ndom,link)
 c-----------------------------------------------------------------------
 c     begin level-set loop
 c-----------------------------------------------------------------------
@@ -3141,7 +3028,7 @@ c
       do 2 kk=ia(i), ia(i+1)-1
          j = ja(kk)
          if (marker(j) .eq. 0) then
-            call add_lk(j,nod,idom,ndom,lkend,levst,link,nodes,marker)
+            call add_lk(j,nod,idom,lkend,levst,link,nodes,marker)
          endif
  2    continue
 c
@@ -3192,11 +3079,11 @@ c
 c
 c     if no neighboring domain select smallest one
 c
-               if (idom .eq. 0) idom = mindom(n,ndom,levst,link)
+               if (idom .eq. 0) idom = mindom(n,ndom,link)
 c
 c     add ii to sudomain idom at end of linked list
 c
-            call add_lk(ii,nod,idom,ndom,lkend,levst,link,nodes,marker)
+            call add_lk(ii,nod,idom,lkend,levst,link,nodes,marker)
             goto 3
          else
             goto 20
@@ -3242,9 +3129,9 @@ c
 c-----------------------------------------------------------------------
       end
 c-----------------------------------------------------------------------
-      function mindom(n, ndom, levst, link)
+      function mindom(n, ndom, link)
       implicit none
-      integer mindom, n, ndom, levst(2*ndom),link(n)
+      integer mindom, n, ndom, link(n)
 c-----------------------------------------------------------------------
 c     returns  the domain with smallest size
 c-----------------------------------------------------------------------
@@ -3252,6 +3139,7 @@ c      locals
 c
       integer i, nsize, isiz
 c
+      mindom = -1
       isiz = n+1
       do 10 i=1, ndom
          nsize = - link(i)
@@ -3264,9 +3152,9 @@ c
       return
       end
 c-----------------------------------------------------------------------
-      subroutine add_lk(new,nod,idom,ndom,lkend,levst,link,nodes,marker)
+      subroutine add_lk(new,nod,idom,lkend,levst,link,nodes,marker)
       implicit none
-      integer new,nod,idom,ndom,lkend,levst(*),link(*),nodes(*),
+      integer new,nod,idom,lkend,levst(*),link(*),nodes(*),
      *     marker(*)
 c-----------------------------------------------------------------------
 c     adds from head --
@@ -3275,7 +3163,6 @@ c     adds one entry (new) to linked list and ipdates everything.
 c     new  = node to be added
 c     nod  = current number of marked nodes
 c     idom = domain to which new is to be added
-c     ndom = total number of domains
 c     lkend= location of end of structure (link and nodes)
 c     levst= pointer array for link, nodes
 c     link = link array
diff --git a/src/spm.c b/src/spm.c
index 49d1daf94083f5c63a5855f1ab853361b9f8602c..c38e807d860b22bd918f14759133d94dd592fa55 100644
--- a/src/spm.c
+++ b/src/spm.c
@@ -190,6 +190,58 @@ spmUpdateComputedFields( spmatrix_t *spm )
     spm->gnnzexp = spm->nnzexp;
 }
 
+/**
+ *******************************************************************************
+ *
+ * @brief Allocate the arrays of an spm structure with PaStiX internal
+ * allocations.
+ *
+ * This function must be called after initialization of the non-computed fields,
+ * and the call to spmUpdateComputedFields(). It allocates the colptr, rowptr,
+ * dof, and values arrays with PaStiX C malloc function. This is highly
+ * recommended to use this function when using PaStiX from Fortran.
+ *
+ *******************************************************************************
+ *
+ * @param[inout] spm
+ *          The sparse matrix fr which the internal arrays needs to be allocated.
+ *
+ *******************************************************************************/
+void
+spmAlloc( spmatrix_t *spm )
+{
+    spm_int_t colsize, rowsize, valsize, dofsize;
+
+    switch(spm->fmttype){
+    case SpmCSC:
+        colsize = spm->n + 1;
+        rowsize = spm->nnz;
+        valsize = spm->nnzexp;
+        dofsize = spm->n + 1;
+        break;
+    case SpmCSR:
+        colsize = spm->nnz;
+        rowsize = spm->n + 1;
+        valsize = spm->nnzexp;
+        dofsize = spm->n + 1;
+        break;
+    case SpmIJV:
+    default:
+        colsize = spm->nnz;
+        rowsize = spm->nnz;
+        valsize = spm->nnzexp;
+        dofsize = spm->n + 1;
+    }
+
+    spm->colptr = (spm_int_t*)malloc( colsize * sizeof(spm_int_t) );
+    spm->rowptr = (spm_int_t*)malloc( rowsize * sizeof(spm_int_t) );
+    if ( spm->dof < 1 ) {
+        spm->dofs = (spm_int_t*)malloc( dofsize * sizeof(spm_int_t) );
+    }
+    valsize = valsize * spm_size_of( spm->flttype );
+    spm->values = malloc(valsize);
+}
+
 /**
  *******************************************************************************
  *
@@ -676,24 +728,28 @@ spmSymmetrize( spmatrix_t *spm )
  *
   *******************************************************************************
  *
- * @param[inout] spm
+ * @param[in] in
  *          The pointer to the sparse matrix structure to check, and correct.
  *
+ * @param[inout] out
+ *          On entry, an allocated structure to hold the corrected spm.
+ *          On exit, holds the pointer to spm corrected.
+ *
  *******************************************************************************
  *
- * @return if no changes have been made, the initial sparse matrix is returned,
- * otherwise a copy of the sparse matrix, fixed to meet the PaStiX requirements,
- * is returned.
+ * @return 0 if no changes have been made to the spm matrix.
+ * @return 1 if corrections have been applied to the in sparse matrix.
  *
  *******************************************************************************/
-spmatrix_t *
-spmCheckAndCorrect( spmatrix_t *spm )
+int
+spmCheckAndCorrect( const spmatrix_t *spm_in,
+                    spmatrix_t       *spm_out )
 {
     spmatrix_t *newspm = NULL;
     spm_int_t count;
 
     /* Let's work on a copy */
-    newspm = spmCopy( spm );
+    newspm = spmCopy( spm_in );
 
     /* PaStiX works on CSC matrices */
     spmConvert( SpmCSC, newspm );
@@ -733,15 +789,18 @@ spmCheckAndCorrect( spmatrix_t *spm )
      * Check if we return the new one, or the original one because no changes
      * have been made
      */
-    if (( spm->fmttype != newspm->fmttype ) ||
-        ( spm->nnzexp  != newspm->nnzexp  ) )
+    if (( spm_in->fmttype != newspm->fmttype ) ||
+        ( spm_in->nnzexp  != newspm->nnzexp  ) )
     {
-        return newspm;
+        memcpy( spm_out, newspm, sizeof(spmatrix_t) );
+        free( newspm );
+        return 1;
     }
     else {
+        memcpy( spm_out, spm_in, sizeof(spmatrix_t) );
         spmExit( newspm );
-        free(newspm);
-        return (spmatrix_t*)spm;
+        free( newspm );
+        return 0;
     }
 }
 
@@ -815,6 +874,7 @@ spmCopy( const spmatrix_t *spm )
         newspm->values = malloc(valsize);
         memcpy( newspm->values, spm->values, valsize );
     }
+
     return newspm;
 }
 
diff --git a/src/spm_io.c b/src/spm_io.c
index 11da33cfec0c5bfcc6ea690f7047f2cd56e1108e..9d536cbe787b86a5d25865fec9e6b0b35e5da258 100644
--- a/src/spm_io.c
+++ b/src/spm_io.c
@@ -508,14 +508,14 @@ spmLoad( spmatrix_t  *spm,
             return SPM_ERR_BADPARAMETER;
         }
 
-        spm->mtxtype = mtxtype;
-        spm->flttype = flttype;
-        spm->fmttype = fmttype;
+        spm->mtxtype = (spm_mtxtype_t)mtxtype;
+        spm->flttype = (spm_coeftype_t)flttype;
+        spm->fmttype = (spm_fmttype_t)fmttype;
         spm->gN      = gN;
         spm->n       = n;
         spm->nnz     = nnz;
         spm->dof     = dof;
-        spm->layout  = layout;
+        spm->layout  = (spm_layout_t)layout;
 
         spmUpdateComputedFields( spm );
 
diff --git a/src/z_spm.c b/src/z_spm.c
index e0d45d07d69f3b8f60789ba76de3bd13642c111d..457e7ba94c7eabe5e9b4fa07a4d84138efe2bfab 100644
--- a/src/z_spm.c
+++ b/src/z_spm.c
@@ -139,8 +139,11 @@ z_spmMergeDuplicate( spmatrix_t *spm )
     spm_int_t n       = spm->n;
     spm_int_t baseval = spm->colptr[0];
     spm_int_t dof2    = spm->dof * spm->dof;
-    spm_int_t idx, i, j, d, size;
+    spm_int_t idx, i, j, size;
     spm_int_t merge = 0;
+#if !defined(PRECISION_p)
+    spm_int_t d;
+#endif
 
     if ( spm->fmttype == SpmCSC ) {
         idx = 0;
@@ -196,7 +199,6 @@ z_spmMergeDuplicate( spmatrix_t *spm )
         }
     }
 
-    (void)d;
     return merge;
 }
 
diff --git a/src/z_spm_expand.c b/src/z_spm_expand.c
index 78a3785cfee037512528ca263bd287905f245737..0fa2d27828061b8e60930e4b61dc60a6e94c4b63 100644
--- a/src/z_spm_expand.c
+++ b/src/z_spm_expand.c
@@ -46,7 +46,10 @@ z_spmCSCExpand(const spmatrix_t *spm)
     spm_int_t        i, j, k, ii, jj, dofi, dofj, col, row, baseval, lda;
     spm_int_t        diag, height;
     spm_int_t       *newcol, *newrow, *oldcol, *oldrow, *dofs;
-    spm_complex64_t *newval, *oldval, *oldval2;
+#if !defined(PRECISION_p)
+    spm_complex64_t *newval = NULL;
+#endif
+    spm_complex64_t *oldval2, *oldval = NULL;
 
     assert( spm->fmttype == SpmCSC );
     assert( spm->flttype == SpmComplex64 );
@@ -200,7 +203,6 @@ z_spmCSCExpand(const spmatrix_t *spm)
 
     assert(spm->loc2glob == NULL);
 
-    (void)newval;
     (void)lda;
     return newspm;
 }
@@ -234,7 +236,8 @@ z_spmCSRExpand(const spmatrix_t *spm)
     spm_int_t        i, j, k, ii, jj, dofi, dofj, col, row, baseval, lda;
     spm_int_t        diag, height;
     spm_int_t       *newcol, *newrow, *oldcol, *oldrow, *dofs;
-    spm_complex64_t *newval, *oldval, *oldval2;
+    spm_complex64_t *newval = NULL;
+    spm_complex64_t *oldval2, *oldval = NULL;
 
     assert( spm->fmttype == SpmCSR );
     assert( spm->flttype == SpmComplex64 );
@@ -419,8 +422,10 @@ z_spmIJVExpand(const spmatrix_t *spm)
     spmatrix_t       *newspm;
     spm_int_t        i, j, k, ii, jj, dofi, dofj, col, row, baseval;
     spm_int_t       *newcol, *newrow, *oldcol, *oldrow, *dofs;
-    spm_complex64_t *newval, *oldval=NULL;
-
+#if !defined(PRECISION_p)
+    spm_complex64_t *newval = NULL;
+#endif
+    spm_complex64_t *oldval = NULL;
     assert( spm->fmttype == SpmIJV );
     assert( spm->flttype == SpmComplex64 );
 
@@ -566,7 +571,6 @@ z_spmIJVExpand(const spmatrix_t *spm)
 
     assert(spm->loc2glob == NULL);
 
-    (void)newval;
     return newspm;
 }
 
diff --git a/src/z_spm_genrhs.c b/src/z_spm_genrhs.c
index e5b4783d7c94350d6d80d011a3e94a0f85448eb3..76f3a248e5403f56b6bd9ff4ed5ed5556ba0514b 100644
--- a/src/z_spm_genrhs.c
+++ b/src/z_spm_genrhs.c
@@ -451,7 +451,7 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
      */
     if ( x0 != NULL ) {
         double normX0;
-        double forw, nr, nx;
+        double forw, nr, nx, nx0;
         int fail;
 
         forward = 0.;
@@ -461,14 +461,15 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
 
         for( i=0; i<nrhs; i++, zx += ldx, zx0 += ldx0 ) {
 
-            nx = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, zx0, ldx0 );
+            nx0 = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, zx0, ldx0 );
+            nx  = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, zx,  ldx  );
 
             cblas_zaxpy( spm->n, CBLAS_SADDR(mzone),
                          zx, 1, zx0, 1);
 
             nr = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, zx0, ldx0 );
 
-            forw = (nr / nx) / eps;
+            forw = (nr / nx0) / eps;
 
             normX0  = ( nx   > normX0  ) ? nx   : normX0;
             normR   = ( nr   > normR   ) ? nr   : normR;
@@ -477,9 +478,11 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
             fail = isnan(nx) || isinf(nx) || isnan(forw) || isinf(forw) || (forw > 1.e2);
             if ( fail ) {
                 fprintf( stdout,
-                         "   || x0_%d ||_oo                                           %e\n"
+                         "   || x_%d ||_oo                                            %e\n"
+                         "   || x0_%d - x_%d ||_oo                                     %e\n"
                          "   || x0_%d - x_%d ||_oo / (||x0_%d||_oo * eps)               %e (%s)\n",
-                         i, nr,
+                         i, nx,
+                         i, i, nr,
                          i, i, i, forw,
                          fail ? "FAILED" : "SUCCESS" );
             }
@@ -488,7 +491,7 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
         }
 
         fprintf( stdout,
-                 "   max(|| x0_i ||_oo)                                      %e\n"
+                 "   max(|| x_i ||_oo)                                       %e\n"
                  "   max(|| x0_i - x_i ||_oo)                                %e\n"
                  "   max(|| x0_i - x_i ||_oo / || x0_i ||_oo)                %e (%s)\n",
                  normX0, normR, forward,
diff --git a/src/z_spm_laplacian.c b/src/z_spm_laplacian.c
index 68fc4361946e9586f1bd505e882602551e49afb2..135e61cd3c2110ac1c4aa9fc4e93a1d148c2ac3e 100644
--- a/src/z_spm_laplacian.c
+++ b/src/z_spm_laplacian.c
@@ -310,10 +310,10 @@ z_spmLaplacian_27points( spmatrix_t   *spm,
 
                 /* Diagonal value */
                 *rowptr = l;
-#if !defined(PRECISION_p)
                 degree = 1;
                 d = 1;
 
+#if !defined(PRECISION_p)
                 if (k > 1) {
                     d++;
                 }
diff --git a/tests/get_options.c b/tests/get_options.c
index 8d8a227b501cb15fed689b2213b5268ca3866b54..683d8da0e45b2aab2752a4823f0274ebdc7c78bf 100644
--- a/tests/get_options.c
+++ b/tests/get_options.c
@@ -105,7 +105,7 @@ spmGetOptions( int argc, char **argv,
         exit(0);
     }
 
-    *driver = -1;
+    *driver = (spm_driver_t)-1;
     do
     {
 #if defined(HAVE_GETOPT_LONG)
diff --git a/tests/spm_convert_tests.c b/tests/spm_convert_tests.c
index 1bc33bb8e7d9fd367432300cac4343c52a7ce890..d644cd10adef6bd896ad3c9cf2353820a0addfe3 100644
--- a/tests/spm_convert_tests.c
+++ b/tests/spm_convert_tests.c
@@ -45,7 +45,7 @@ int spmComp( const spmatrix_t *spm1,
 {
     spm_int_t *colptr1, *colptr2;
     spm_int_t *rowptr1, *rowptr2;
-    int          *valptr1, *valptr2;
+    int       *valptr1, *valptr2;
     spm_int_t  i;
 
     if ( spm1->fmttype != SpmCSC ) {
@@ -108,7 +108,8 @@ int main (int argc, char **argv)
     char *filename;
     spmatrix_t  spm, *spm2;
     spm_driver_t driver;
-    int mtxtype, baseval;
+    spm_mtxtype_t mtxtype;
+    int baseval;
     int ret = SPM_SUCCESS;
     int err = 0;
     FILE *f;
diff --git a/tests/spm_dof_expand_tests.c b/tests/spm_dof_expand_tests.c
index 0a343aa99c7f140b4f11edae9a4aea8db3abf9d7..e27e495cdfc3fe873ed0274d726f880f9863b523 100644
--- a/tests/spm_dof_expand_tests.c
+++ b/tests/spm_dof_expand_tests.c
@@ -54,8 +54,10 @@ int main (int argc, char **argv)
     spmatrix_t    original, *spm;
     spm_driver_t driver;
     char *filename;
-    int spmtype, mtxtype, fmttype, baseval;
-    int i, rc, dofmax = 3;
+    spm_mtxtype_t spmtype;
+    spm_mtxtype_t mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval, i, rc, dofmax = 3;
 
     /**
      * Get options from command line
@@ -95,7 +97,7 @@ int main (int argc, char **argv)
             {
                 spmBase( &original, baseval );
 
-                for( fmttype=0; fmttype<3; fmttype++ )
+                for( fmttype=SpmCSC; fmttype<=SpmIJV; fmttype++ )
                 {
                     spmConvert( fmttype, &original );
                     spm = spmDofExtend( &original, i, dofmax );
diff --git a/tests/spm_dof_matvec_tests.c b/tests/spm_dof_matvec_tests.c
index 8d5e1a9891211ea35bc3d4af9668b4fb5778daf6..41723259ba0f4159a88771cad01df24a69d608db 100644
--- a/tests/spm_dof_matvec_tests.c
+++ b/tests/spm_dof_matvec_tests.c
@@ -23,11 +23,6 @@
 #include <time.h>
 #include "spm_tests.h"
 
-int z_spm_matvec_check( int trans, const spmatrix_t *spm );
-int c_spm_matvec_check( int trans, const spmatrix_t *spm );
-int d_spm_matvec_check( int trans, const spmatrix_t *spm );
-int s_spm_matvec_check( int trans, const spmatrix_t *spm );
-
 #define PRINT_RES(_ret_)                        \
     if(_ret_) {                                 \
         printf("FAILED(%d)\n", _ret_);          \
@@ -46,7 +41,9 @@ int main (int argc, char **argv)
     spmatrix_t    original, *spm;
     spm_driver_t driver;
     char *filename;
-    int spmtype, mtxtype, fmttype, baseval;
+    spm_mtxtype_t spmtype, mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval;
     int ret = SPM_SUCCESS;
     int err = 0;
     int rc, i, dofmax = 3;
diff --git a/tests/spm_dof_norm_tests.c b/tests/spm_dof_norm_tests.c
index 6d4e74f6a228fa6534f6235d63ec8347a2d96021..fe229a3a207308c6c50ac77787ae3712e51ee20f 100644
--- a/tests/spm_dof_norm_tests.c
+++ b/tests/spm_dof_norm_tests.c
@@ -43,9 +43,11 @@ char* mtxnames[] = { "General", "Symmetric", "Hermitian" };
 int main (int argc, char **argv)
 {
     spmatrix_t    original, *spm;
-    spm_driver_t driver;
-    char *filename;
-    int spmtype, mtxtype, fmttype, baseval;
+    spm_driver_t  driver;
+    char         *filename;
+    spm_mtxtype_t spmtype, mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval;
     int rc = SPM_SUCCESS;
     int err = 0;
     int i, dofmax = 4;
@@ -92,7 +94,7 @@ int main (int argc, char **argv)
             {
                 spmBase( &original, baseval );
 
-                for( fmttype=0; fmttype<3; fmttype++ )
+                for( fmttype=SpmCSC; fmttype<=SpmIJV; fmttype++ )
                 {
                     spmConvert( fmttype, &original );
                     spm = spmDofExtend( &original, i, dofmax );
diff --git a/tests/spm_matvec_tests.c b/tests/spm_matvec_tests.c
index 1a7c5e40fb6b72717cd41b2d02d4edd68a6b331a..2ccd492b949d17475d4468deeea9139ed5e92a59 100644
--- a/tests/spm_matvec_tests.c
+++ b/tests/spm_matvec_tests.c
@@ -22,11 +22,6 @@
 #include <time.h>
 #include "spm_tests.h"
 
-int z_spm_matvec_check( int trans, const spmatrix_t *spm );
-int c_spm_matvec_check( int trans, const spmatrix_t *spm );
-int d_spm_matvec_check( int trans, const spmatrix_t *spm );
-int s_spm_matvec_check( int trans, const spmatrix_t *spm );
-
 #define PRINT_RES(_ret_)                        \
     if(_ret_) {                                 \
         printf("FAILED(%d)\n", _ret_);          \
@@ -46,7 +41,11 @@ int main (int argc, char **argv)
     spmatrix_t    spm;
     spm_driver_t driver;
     char *filename;
-    int t,spmtype, mtxtype, fmttype, baseval;
+    spm_trans_t t;
+    spm_mtxtype_t spmtype;
+    spm_mtxtype_t mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval;
     int rc = SPM_SUCCESS;
     int err = 0;
 
diff --git a/tests/spm_norm_tests.c b/tests/spm_norm_tests.c
index d6824b737159545edd007a8c386492753b5bfe97..b2ae5a42b278652b1a223ea980d9339fb49f1092 100644
--- a/tests/spm_norm_tests.c
+++ b/tests/spm_norm_tests.c
@@ -45,7 +45,9 @@ int main (int argc, char **argv)
     spmatrix_t    spm;
     spm_driver_t driver;
     char *filename;
-    int spmtype, mtxtype, fmttype, baseval;
+    spm_mtxtype_t spmtype, mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval;
     int rc = SPM_SUCCESS;
     int err = 0;
 
@@ -70,7 +72,7 @@ int main (int argc, char **argv)
     spmtype = spm.mtxtype;
     printf(" -- SPM Norms Test --\n");
 
-    for( fmttype=0; fmttype<3; fmttype++ ) {
+    for( fmttype=SpmCSC; fmttype<=SpmIJV; fmttype++ ) {
 
         spmConvert( fmttype, &spm );
 
diff --git a/tests/spm_tests.h b/tests/spm_tests.h
index 765de2347d2cb09665554e57245af0d11229aa59..7d79984f5407899dbc2dee1fa8343c398aed8630 100644
--- a/tests/spm_tests.h
+++ b/tests/spm_tests.h
@@ -72,19 +72,19 @@ int  core_sgeadd( spm_trans_t  trans,
                   spm_int_t    LDB );
 
 void z_spm_print_check( char *filename, const spmatrix_t *spm );
-int  z_spm_matvec_check( int trans, const spmatrix_t *spm );
+int  z_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  z_spm_norm_check( const spmatrix_t *spm );
 
 void c_spm_print_check( char *filename, const spmatrix_t *spm );
-int  c_spm_matvec_check( int trans, const spmatrix_t *spm );
+int  c_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  c_spm_norm_check( const spmatrix_t *spm );
 
 void d_spm_print_check( char *filename, const spmatrix_t *spm );
-int  d_spm_matvec_check( int trans, const spmatrix_t *spm );
+int  d_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  d_spm_norm_check( const spmatrix_t *spm );
 
 void s_spm_print_check( char *filename, const spmatrix_t *spm );
-int  s_spm_matvec_check( int trans, const spmatrix_t *spm );
+int  s_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  s_spm_norm_check( const spmatrix_t *spm );
 
 static inline int
diff --git a/tests/z_spm_tests.c b/tests/z_spm_tests.c
index f6487f297f19e46739b0a069835da3fdb02381d5..17f0bb8c3d32cc552a617026bd60310d21fe8448 100644
--- a/tests/z_spm_tests.c
+++ b/tests/z_spm_tests.c
@@ -109,7 +109,7 @@ z_spm_print_check( char *filename, const spmatrix_t *spm )
  *  Check the accuracy of the solution
  */
 int
-z_spm_matvec_check( int trans, const spmatrix_t *spm )
+z_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm )
 {
     unsigned long long int seed = 35469;
     spm_complex64_t *A, *x, *y0, *ys, *yd;
@@ -161,7 +161,8 @@ z_spm_matvec_check( int trans, const spmatrix_t *spm )
     }
 
     /* Compute the dense matrix-vector product */
-    cblas_zgemm( CblasColMajor, trans, CblasNoTrans, spm->gNexp, 1, spm->gNexp,
+    cblas_zgemm( CblasColMajor, (CBLAS_TRANSPOSE)trans, CblasNoTrans,
+                 spm->gNexp, 1, spm->gNexp,
                  CBLAS_SADDR(zalpha), A, spm->gNexp,
                                       x, spm->gNexp,
                  CBLAS_SADDR(zbeta), yd, spm->gNexp );
diff --git a/tools/analysis.sh b/tools/analysis.sh
index 0ec0039d64e1a5c75312033add6de7c08da777a8..837391731bf290ce4d248c04ed8b44650b497c9d 100755
--- a/tools/analysis.sh
+++ b/tools/analysis.sh
@@ -62,6 +62,7 @@ sonar.c.compiler.reportPath=spm-build.log
 sonar.c.coverage.reportPath=spm-coverage.xml
 sonar.c.cppcheck.reportPath=spm-cppcheck.xml
 sonar.c.rats.reportPath=spm-rats.xml
+sonar.c.jsonCompilationDatabase=build/compile_commands.json
 EOF
 
 # run sonar analysis + publish on sonarqube-dev
diff --git a/tools/gen_wrappers.py b/tools/gen_wrappers.py
index e03a3c5bd0fbd8bea0bdbd3ca1cbb273c80444d2..b923b517a97f8debe7f68147762df5c9c5dd30a2 100755
--- a/tools/gen_wrappers.py
+++ b/tools/gen_wrappers.py
@@ -450,14 +450,13 @@ enums_python_coeftype='''
             return -1
 '''
 
-enums_fortran_footer='''
-  integer, parameter :: spm_int_t = SPM_INT_KIND
+enums_fortran_footer='''  integer, parameter :: spm_int_t = SPM_INT_KIND
 
 contains
 
   function spm_getintsize()
     integer :: spm_getintsize
-    spm_getintsize = SPM_INT_KIND
+    spm_getintsize = kind(SPM_INT_KIND)
     return
   end function spm_getintsize
 '''
diff --git a/tools/spm.pc.in b/tools/spm.pc.in
index 56c3d7b31d2291cf3affc63a1afa4c1ec20887e7..42bea7f45f6f22a015cfe73c7380b3d302eeb453 100644
--- a/tools/spm.pc.in
+++ b/tools/spm.pc.in
@@ -7,7 +7,7 @@ Name: SPM
 Description: SParse Matrix package
 Version: @SPM_VERSION_MAJOR@.@SPM_VERSION_MINOR@.@SPM_VERSION_MICRO@
 Cflags: -I${includedir}
-Libs: -L${libdir} @SPM_PKGCONFIG_LIBS@
+Libs: -L${libdir} -lspm @SPM_PKGCONFIG_LIBS@
 Libs.private: @SPM_PKGCONFIG_LIBS_PRIVATE@
 Requires: @SPM_PKGCONFIG_REQUIRED@
 Requires.private: @SPM_PKGCONFIG_REQUIRED_PRIVATE@
diff --git a/tools/spmf.pc.in b/tools/spmf.pc.in
new file mode 100644
index 0000000000000000000000000000000000000000..c368a6d7fb3371d6a3303a931b9bc0ea2b119484
--- /dev/null
+++ b/tools/spmf.pc.in
@@ -0,0 +1,23 @@
+#
+#  @file spmf.pc
+#
+#  @copyright 2016-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+#                       Univ. Bordeaux. All rights reserved.
+#
+#  @version 1.0.0
+#  @author Mathieu Faverge
+#  @date 2017-06-24
+#
+prefix=@CMAKE_INSTALL_PREFIX@
+exec_prefix=${prefix}
+libdir=${exec_prefix}/lib
+includedir=${exec_prefix}/include
+
+Name: SPM - Fortran
+Description: SParse Matrix package - Fortran interface
+Version: @SPM_VERSION_MAJOR@.@SPM_VERSION_MINOR@.@SPM_VERSION_MICRO@
+Cflags: -I${includedir}
+Libs: -L${libdir} -lspmf @SPM_PKGCONFIG_LIBS@
+Libs.private: @SPM_PKGCONFIG_LIBS_PRIVATE@
+Requires: spm @SPM_PKGCONFIG_REQUIRED@
+Requires.private: @SPM_PKGCONFIG_REQUIRED_PRIVATE@
diff --git a/tools/wrappers/wrap_fortran.py b/tools/wrappers/wrap_fortran.py
index e4250f832399f6fab099d1d1086cc0807713a482..69242d335e57953b8fd69f7f14845b5212e8f12f 100644
--- a/tools/wrappers/wrap_fortran.py
+++ b/tools/wrappers/wrap_fortran.py
@@ -30,6 +30,7 @@ iindent=3
 # translation_table of types
 types_dict = {
     "int":            ("integer(kind=c_int)"),
+    "int8_t":         ("integer(kind=c_int8_t)"),
     "spm_coeftype_t": ("integer(c_int)"),
     "spm_dir_t":      ("integer(c_int)"),
     "spm_trans_t":    ("integer(c_int)"),
@@ -87,7 +88,7 @@ def iso_c_wrapper_type(arg, args_list, args_size):
     else:
         is_pointer = False
 
-    if (is_pointer and arg[0].strip() == "void"):
+    if (is_pointer and (arg[0].strip() == "void")):
         f_type = "type(c_ptr), "
     else:
         f_type = types_dict[arg[0]] + ", "
@@ -108,11 +109,6 @@ def iso_c_wrapper_type(arg, args_list, args_size):
 
     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 = "(:,:)"
diff --git a/tools/wrappers/wrap_python.py b/tools/wrappers/wrap_python.py
index 60945858ffae065efb42db29682bb9d4d9952044..44a75e53cd466d93e7fbc5b2edd22801ed34c684 100644
--- a/tools/wrappers/wrap_python.py
+++ b/tools/wrappers/wrap_python.py
@@ -26,6 +26,7 @@ iindent=4
 # translation_table of types
 types_dict = {
     "int":            ("c_int"),
+    "int8_t":         ("c_int8"),
     "spm_coeftype_t": ("c_int"),
     "spm_dir_t":      ("c_int"),
     "spm_trans_t":    ("c_int"),
diff --git a/wrappers/fortran90/CMakeLists.txt b/wrappers/fortran90/CMakeLists.txt
index dd80adc3ae8ba3c1e84ffe8e6242a5335df6732a..a7a81fea05e7f7e8d6a9a10fdd6ed3d47d04dffa 100644
--- a/wrappers/fortran90/CMakeLists.txt
+++ b/wrappers/fortran90/CMakeLists.txt
@@ -21,11 +21,11 @@ add_library( spmf
 if ( SPM_INT64 )
   set_source_files_properties(
     src/spm_enums.F90
-    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=8")
+    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=c_int64_t")
 else()
   set_source_files_properties(
     src/spm_enums.F90
-    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=4")
+    PROPERTIES COMPILE_DEFINITIONS "SPM_INT_KIND=c_int32_t")
 endif()
 
 target_link_libraries( spmf spm )
diff --git a/wrappers/fortran90/examples/spm_driver.f90 b/wrappers/fortran90/examples/spm_driver.f90
index fc7f05e450898a1801c52a72aebcf95cd2121c48..3e856d21a1e79150ffd076e11bdf5dae38f2ae6f 100644
--- a/wrappers/fortran90/examples/spm_driver.f90
+++ b/wrappers/fortran90/examples/spm_driver.f90
@@ -17,7 +17,7 @@ program spm_driver
   implicit none
 
   type(spmatrix_t),           target                       :: spm
-  type(spmatrix_t),           pointer                      :: spm2
+  type(spmatrix_t),           target                       :: spm2
   real(kind=c_double)                                      :: normA
   real(kind=c_double)                                      :: eps = 1.e-15
   integer(c_int)                                           :: info
@@ -30,8 +30,8 @@ program spm_driver
   !   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 spmCheckAndCorrect( spm, spm2, info )
+  if ( info .ne. 0 ) then
      call spmExit( spm )
      spm = spm2
   end if
diff --git a/wrappers/fortran90/examples/spm_user.f90 b/wrappers/fortran90/examples/spm_user.f90
index cdf8df9bb84dd30ed8f0de8d9221f460661fb0a6..833bd4099faaa0cc120d0f453972351fcdbadffc 100644
--- a/wrappers/fortran90/examples/spm_user.f90
+++ b/wrappers/fortran90/examples/spm_user.f90
@@ -15,14 +15,14 @@ program spm_user
   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
+  integer(kind=spm_int_t), dimension(:), pointer :: rowptr
+  integer(kind=spm_int_t), dimension(:), pointer :: colptr
+  real(kind=c_double),  dimension(:), pointer    :: 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
+  type(spmatrix_t),        target                            :: 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
@@ -36,9 +36,23 @@ program spm_user
   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))
+  !
+  ! 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
+
+  call spmUpdateComputedFields( spm )
+  call spmAlloc( spm )
+
+  call c_f_pointer( spm%rowptr, rowptr, [nnz] )
+  call c_f_pointer( spm%colptr, colptr, [nnz] )
+  call c_f_pointer( spm%values, values, [nnz] )
 
   l = 1
   do i=1,dim1
@@ -73,19 +87,19 @@ program spm_user
            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
+              values(l) = - 1.
               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
+              values(l) = - 1.
               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
+              values(l) = -1.
               l = l + 1
            end if
         end do
@@ -96,32 +110,8 @@ program spm_user
      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 spmCheckAndCorrect( spm, spm2, info )
+  if ( info .ne. 0 ) then
      call spmExit( spm )
      spm = spm2
   end if
diff --git a/wrappers/fortran90/src/spm_enums.F90 b/wrappers/fortran90/src/spm_enums.F90
index 5b7f9c93d6cb68a2aa2c0b5f54e95b03a24d5330..c1919678628cb014fa1ade67fccd8269db9a26aa 100644
--- a/wrappers/fortran90/src/spm_enums.F90
+++ b/wrappers/fortran90/src/spm_enums.F90
@@ -129,14 +129,13 @@ module spm_enums
      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
+    spm_getintsize = kind(SPM_INT_KIND)
     return
   end function spm_getintsize
 
diff --git a/wrappers/fortran90/src/spmf.f90 b/wrappers/fortran90/src/spmf.f90
index f3c722a02063b2f313d280036ff7f144651b4115..ac0ce6d69b98cb3e65e68debaccacbe0810f8c62 100644
--- a/wrappers/fortran90/src/spmf.f90
+++ b/wrappers/fortran90/src/spmf.f90
@@ -49,6 +49,16 @@ module spmf
      end subroutine spmInit_c
   end interface
 
+  interface
+     subroutine spmAlloc_c(spm) &
+          bind(c, name='spmAlloc')
+       use iso_c_binding
+       import spmatrix_t
+       implicit none
+       type(c_ptr), value :: spm
+     end subroutine spmAlloc_c
+  end interface
+
   interface
      subroutine spmExit_c(spm) &
           bind(c, name='spmExit')
@@ -234,13 +244,14 @@ module spmf
   end interface
 
   interface
-     function spmCheckAndCorrect_c(spm) &
+     function spmCheckAndCorrect_c(spm_in, spm_out) &
           bind(c, name='spmCheckAndCorrect')
        use iso_c_binding
        import spmatrix_t
        implicit none
-       type(c_ptr)        :: spmCheckAndCorrect_c
-       type(c_ptr), value :: spm
+       integer(kind=c_int)   :: spmCheckAndCorrect_c
+       type(c_ptr),    value :: spm_in
+       type(c_ptr),    value :: spm_out
      end function spmCheckAndCorrect_c
   end interface
 
@@ -416,6 +427,14 @@ contains
     call spmInit_c(c_loc(spm))
   end subroutine spmInit
 
+  subroutine spmAlloc(spm)
+    use iso_c_binding
+    implicit none
+    type(spmatrix_t), intent(inout), target :: spm
+
+    call spmAlloc_c(c_loc(spm))
+  end subroutine spmAlloc
+
   subroutine spmExit(spm)
     use iso_c_binding
     implicit none
@@ -566,13 +585,14 @@ contains
     value = spmSymmetrize_c(c_loc(spm))
   end subroutine spmSymmetrize
 
-  subroutine spmCheckAndCorrect(spm, spmo)
+  subroutine spmCheckAndCorrect(spm_in, spm_out, info)
     use iso_c_binding
     implicit none
-    type(spmatrix_t), intent(inout), target  :: spm
-    type(spmatrix_t), intent(out),   pointer :: spmo
+    type(spmatrix_t),    intent(in),    target :: spm_in
+    type(spmatrix_t),    intent(inout), target :: spm_out
+    integer(kind=c_int), intent(out)           :: info
 
-    call c_f_pointer(spmCheckAndCorrect_c(c_loc(spm)), spmo)
+    info = spmCheckAndCorrect_c(c_loc(spm_in), c_loc(spm_out))
   end subroutine spmCheckAndCorrect
 
   subroutine spmGenRHS(type, nrhs, spm, x, ldx, b, ldb, info)
diff --git a/wrappers/python/spm/__spm__.py b/wrappers/python/spm/__spm__.py
index a4b902fb73541cadc99f7b147ab18854acc6e4d3..97148fc180fda384e3528f2f8a1384f8ed4c93a1 100644
--- a/wrappers/python/spm/__spm__.py
+++ b/wrappers/python/spm/__spm__.py
@@ -46,6 +46,10 @@ def pyspm_spmInit( spm ):
     libspm.spmInit.argtypes = [ POINTER(pyspm_spmatrix_t) ]
     libspm.spmInit( spm )
 
+def pyspm_spmAlloc( spm ):
+    libspm.spmAlloc.argtypes = [ POINTER(pyspm_spmatrix_t) ]
+    libspm.spmAlloc( spm )
+
 def pyspm_spmExit( spm ):
     libspm.spmExit.argtypes = [ POINTER(pyspm_spmatrix_t) ]
     libspm.spmExit( spm )
@@ -119,10 +123,11 @@ def pyspm_spmSymmetrize( spm ):
     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_spmCheckAndCorrect( spm_in, spm_out ):
+    libspm.spmCheckAndCorrect.argtypes = [ POINTER(pyspm_spmatrix_t),
+                                           POINTER(pyspm_spmatrix_t) ]
+    libspm.spmCheckAndCorrect.restype = c_int
+    return libspm.spmCheckAndCorrect( spm_in, spm_out )
 
 def pyspm_spmGenRHS( type, nrhs, spm, x, ldx, b, ldb ):
     libspm.spmGenRHS.argtypes = [ c_int, __spm_int__, POINTER(pyspm_spmatrix_t),
diff --git a/wrappers/python/spm/spm.py b/wrappers/python/spm/spm.py
index 016a4cb05f8fb46847b40f980ddb5b11f5e132d2..79aab068b8192e6a33ab8e2701abac7ab1c87ffd 100644
--- a/wrappers/python/spm/spm.py
+++ b/wrappers/python/spm/spm.py
@@ -39,13 +39,14 @@ class spmatrix():
             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 )
+                                       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:
@@ -150,13 +151,13 @@ class spmatrix():
 
     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  ) ):
+        spm2 = spmatrix()
+        rc = pyspm_spmCheckAndCorrect( self.id_ptr, spm2.id_ptr )
+        if ( rc == 0 ):
             return
+        self = spm2
         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()
@@ -225,3 +226,34 @@ class spmatrix():
                          x.ctypes.data_as(c_void_p), ldx,
                          b.ctypes.data_as(c_void_p), ldb )
         return x, b
+
+    def mult( self, B, C, trans=trans.NoTrans, n=-1, alpha=1.0, beta=0. ):
+
+        m = self.spm_c.n
+
+        B = np.asarray(B, self.dtype)
+        C = np.asarray(C, self.dtype)
+
+        if n == -1:
+            if B.ndim == 1:
+                n = 1
+            else:
+                n = B.shape[1]
+
+        self.__checkVector( m, n, B )
+        self.__checkVector( m, n, C )
+
+        ldb  = B.shape[0]
+        ldc  = B.shape[0]
+
+        if n == 1:
+            pyspm_spmMatVec( trans, alpha, self.id_ptr,
+                             B.ctypes.data_as(c_void_p),
+                             beta,
+                             C.ctypes.data_as(c_void_p) )
+        else:
+            pyspm_spmMatMat( trans, n, alpha, self.id_ptr,
+                             B.ctypes.data_as(c_void_p), ldb,
+                             beta,
+                             C.ctypes.data_as(c_void_p), ldc )
+
diff --git a/wrappers/python/spm_driver.py b/wrappers/python/spm_driver.py
index fa88e594052daa326843ce7c7741aa2288392d64..99317ab7a26ff4b66d428fbbc47a6fd393485d8d 100755
--- a/wrappers/python/spm_driver.py
+++ b/wrappers/python/spm_driver.py
@@ -37,6 +37,5 @@ nrhs = 10
 x0, b = A.genRHS( spm.rhstype.RndX, nrhs )
 
 # Check that A * x = b
-x = x0.copy()
-A.checkAxb( x0, b, x )
+A.checkAxb( None, b, x0 )
 
diff --git a/wrappers/python/spm_scipy.py b/wrappers/python/spm_scipy.py
index 6a574b6a5100de854189502139aac8db74e3fb97..d21672380e90a4db5108111a7262e7bcba3d0675 100755
--- a/wrappers/python/spm_scipy.py
+++ b/wrappers/python/spm_scipy.py
@@ -26,11 +26,13 @@ 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)
+b  = np.zeros(n)
 
 spmA = spm.spmatrix( A )
 
+# Multiply A by x
+spmA.mult( x0, b, trans=spm.trans.NoTrans, alpha=1., beta=0. )
+
 # Check that A * x = b
-x = x0.copy()
-spmA.checkAxb( x0, b, x )
+spmA.checkAxb( None, b, x0 )