path: root/r845.patch

Index: setup.py
===================================================================
--- setup.py	(revision 844)
+++ setup.py	(revision 845)
@@ -309,15 +309,30 @@
             self.with_SUNDIALS=True
             L.debug('SUNDIALS found.')
             sundials_version = None
+            sundials_vector_type_size = None
             
             try:
                 if os.path.exists(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')):
                     with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f:
                         for line in f:
-                            if "SUNDIALS_PACKAGE_VERSION" in line:
+                            if "SUNDIALS_PACKAGE_VERSION" in line or "SUNDIALS_VERSION" in line:
                                 sundials_version = tuple([int(f) for f in line.split()[-1][1:-1].split(".")])
                                 L.debug('SUNDIALS %d.%d found.'%(sundials_version[0], sundials_version[1]))
-            except Exception:
+                                break
+                    with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f:
+                        for line in f:
+                            if "SUNDIALS_INT32_T" in line and line.startswith("#define"):
+                                sundials_vector_type_size = "32"
+                                L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size))
+                                break
+                            if "SUNDIALS_INT64_T" in line and line.startswith("#define"):
+                                sundials_vector_type_size = "64"
+                                L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size))
+                                if self.with_SLU:
+                                    L.warning("It is recommended to set the SUNDIALS_INDEX_TYPE to an 32bit integer when using SUNDIALS together with SuperLU.")
+                                    L.warning("SuperLU may not function properly.")
+                                break
+            except Exception as e:
                 if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6
                     sundials_version = (2,6,0)
                     L.debug('SUNDIALS 2.6 found.')
@@ -326,6 +341,7 @@
                     L.debug('SUNDIALS 2.5 found.')
                 
             self.SUNDIALS_version = sundials_version
+            self.SUNDIALS_vector_size = sundials_vector_type_size
             
         else:    
             L.warning(("Could not find Sundials, check the provided path (--sundials-home={}) "+ 
@@ -375,7 +391,8 @@
         # SUNDIALS
         if self.with_SUNDIALS:
             compile_time_env = {'SUNDIALS_VERSION': self.SUNDIALS_version,
-                                'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu}
+                                'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu,
+                                'SUNDIALS_VECTOR_SIZE': self.SUNDIALS_vector_size}
             #CVode and IDA
             ext_list += cythonize(["assimulo" + os.path.sep + "solvers" + os.path.sep + "sundials.pyx"], 
                                  include_path=[".","assimulo","assimulo" + os.sep + "lib"],
@@ -382,11 +399,19 @@
                                  compile_time_env=compile_time_env, force=True)
             ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs]
             ext_list[-1].library_dirs = [self.libdirs]
-            ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"]
+            
+            if self.SUNDIALS_version >= (3,0,0):
+                ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas", "sundials_sunlinsoldense", "sundials_sunlinsolspgmr", "sundials_sunmatrixdense", "sundials_sunmatrixsparse"]
+            else:
+                ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"]
             if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU
+                if self.SUNDIALS_version >= (3,0,0):
+                    ext_list[-1].libraries.extend(["sundials_sunlinsolsuperlumt"])
+                
                 ext_list[-1].include_dirs.append(self.SLUincdir)
                 ext_list[-1].library_dirs.append(self.SLUlibdir)
                 ext_list[-1].libraries.extend(self.superLUFiles)
+                
         
             #Kinsol
             ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"kinsol.pyx"], 
@@ -395,7 +420,11 @@
             ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs]
             ext_list[-1].library_dirs = [self.libdirs]
             ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"]
-    
+            
+            if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU
+                ext_list[-1].include_dirs.append(self.SLUincdir)
+                ext_list[-1].library_dirs.append(self.SLUlibdir)
+                ext_list[-1].libraries.extend(self.superLUFiles)
         
         for el in ext_list:
             #Debug
Index: assimulo/solvers/__init__.py
===================================================================
--- assimulo/solvers/__init__.py	(revision 844)
+++ assimulo/solvers/__init__.py	(revision 845)
@@ -25,7 +25,7 @@
     from .euler import ExplicitEuler
     from .euler import ImplicitEuler
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .radau5 import Radau5ODE
     from .radau5 import Radau5DAE
@@ -32,43 +32,43 @@
     from .radau5 import _Radau5ODE
     from .radau5 import _Radau5DAE 
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .sundials import IDA
     from .sundials import CVode
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .kinsol import KINSOL
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .runge_kutta import RungeKutta34
     from .runge_kutta import RungeKutta4
     from .runge_kutta import Dopri5
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .rosenbrock import RodasODE
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .odassl import ODASSL
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .odepack import LSODAR
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .radar5 import Radar5ODE
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .dasp3 import DASP3ODE
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
 try:
     from .glimda import GLIMDA
 except ImportError as ie:
-    sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+    sys.stderr.write("Could not find " + str(ie) + "\n")
Index: src/solvers/kinsol.pyx
===================================================================
--- assimulo/solvers/kinsol.pyx	(revision 844)
+++ assimulo/solvers/kinsol.pyx	(revision 845)
@@ -30,8 +30,8 @@
 cimport sundials_includes as SUNDIALS
 
 #Various C includes transfered to namespace
-from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
-from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
 from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
 from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
 
@@ -57,6 +57,8 @@
     
     cdef object pt_fcn, pt_jac, pt_jacv, pt_prec_setup, pt_prec_solve
     cdef object _added_linear_solver
+    cdef SUNDIALS.SUNMatrix sun_matrix
+    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
     
     def __init__(self, problem):
         Algebraic.__init__(self, problem) #Calls the base class
@@ -86,6 +88,7 @@
         self.options["no_min_epsilon"] = False #Specifies wheter the scaled linear residual is bounded from below
         self.options["max_beta_fails"] = 10
         self.options["max_krylov"] = 0
+        self.options["precond"] = PREC_NONE
         
         #Statistics
         self.statistics["nfevals"]    = 0 #Function evaluations
@@ -109,6 +112,13 @@
         if self.kinsol_mem != NULL:
             #Free Memory
             SUNDIALS.KINFree(&self.kinsol_mem)
+            
+        IF SUNDIALS_VERSION >= (3,0,0):
+            if self.sun_matrix != NULL:
+                SUNDIALS.SUNMatDestroy(self.sun_matrix)
+                
+            if self.sun_linearsolver != NULL:
+                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
         
     def update_variable_scaling(self, value="Automatic"):
         """
@@ -170,7 +180,7 @@
             
     cdef initialize_kinsol(self):
         cdef int flag #Used for return
-        
+
         self.y_temp  = arr2nv(self.y)
         self.y_scale = arr2nv([1.0]*self.problem_info["dim"])
         self.f_scale = arr2nv([1.0]*self.problem_info["dim"])
@@ -208,16 +218,34 @@
             
     cpdef add_linear_solver(self):
         if self.options["linear_solver"] == "DENSE":
-            flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
+            IF SUNDIALS_VERSION >= (3,0,0):
+                #Create a dense Sundials matrix
+                self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
+                #Create a dense Sundials linear solver
+                self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.y_temp, self.sun_matrix)
+                #Attach it to Kinsol
+                flag = SUNDIALS.KINDlsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver, self.sun_matrix)
+            ELSE:
+                flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
             if flag < 0:
                 raise KINSOLError(flag)
             
             if self.problem_info["jac_fcn"]:
-                flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.KINDlsSetJacFn(self.kinsol_mem, kin_jac);
+                ELSE:
+                    flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
                 if flag < 0:
                     raise KINSOLError(flag)
         elif self.options["linear_solver"] == "SPGMR":
-            flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
+            IF SUNDIALS_VERSION >= (3,0,0):
+                #Create the linear solver
+                self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.y_temp, self.options["precond"], self.options["max_krylov"])
+                #Attach it to Kinsol
+                flag = SUNDIALS.KINSpilsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver)
+            ELSE:
+                #Specify the use of KINSpgmr linear solver.
+                flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
             if flag < 0:
                 raise KINSOLError(flag)
             
Index: src/solvers/sundials.pyx
===================================================================
--- assimulo/solvers/sundials.pyx	(revision 844)
+++ assimulo/solvers/sundials.pyx	(revision 845)
@@ -34,8 +34,8 @@
 cimport sundials_includes as SUNDIALS
 
 #Various C includes transfered to namespace
-from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
-from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
 from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
 from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
 
@@ -70,6 +70,8 @@
     cdef public N.ndarray yS0
     #cdef N.ndarray _event_info
     cdef public N.ndarray g_old
+    cdef SUNDIALS.SUNMatrix sun_matrix
+    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
     
     def __init__(self, problem):
         Implicit_ODE.__init__(self, problem) #Calls the base class
@@ -100,7 +102,9 @@
         self.options["dqrhomax"] = 0.0
         self.options["pbar"] = [1]*self.problem_info["dimSens"]
         self.options["external_event_detection"] = False #Sundials rootfinding is used for event location as default 
+        self.options["precond"] = PREC_NONE
 
+
         #Solver support
         self.supports["report_continuously"] = True
         self.supports["interpolated_output"] = True
@@ -173,6 +177,13 @@
         if self.ida_mem != NULL: 
             #Free Memory
             SUNDIALS.IDAFree(&self.ida_mem)
+        
+        IF SUNDIALS_VERSION >= (3,0,0):
+            if self.sun_matrix != NULL:
+                SUNDIALS.SUNMatDestroy(self.sun_matrix)
+                
+            if self.sun_linearsolver != NULL:
+                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
     
     cpdef state_event_info(self):
         """
@@ -255,21 +266,30 @@
             if flag < 0:
                 raise IDAError(flag, self.t)
                 
-            #Specify the use of the internal dense linear algebra functions.
-            flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
-            if flag < 0:
-                raise IDAError(flag, self.t)
-                
-                    #Choose a linear solver if and only if NEWTON is choosen
+            #Choose a linear solver if and only if NEWTON is choosen
             if self.options["linear_solver"] == 'DENSE':
-                #Specify the use of the internal dense linear algebra functions.
-                flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    #Create a dense Sundials matrix
+                    self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
+                    #Create a dense Sundials linear solver
+                    self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix)
+                    #Attach it to IDA
+                    flag = SUNDIALS.IDADlsSetLinearSolver(self.ida_mem, self.sun_linearsolver, self.sun_matrix);
+                ELSE:
+                    #Specify the use of the internal dense linear algebra functions.
+                    flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
                 if flag < 0:
                     raise IDAError(flag, self.t)
                         
             elif self.options["linear_solver"] == 'SPGMR':
-                #Specify the use of SPGMR linear solver.
-                flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    #Create the linear solver
+                    self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], 0)
+                    #Attach it to IDAS
+                    flag = SUNDIALS.IDASpilsSetLinearSolver(self.ida_mem, self.sun_linearsolver)
+                ELSE:
+                    #Specify the use of SPGMR linear solver.
+                    flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations
                 if flag < 0: 
                     raise IDAError(flag, self.t)
                 
@@ -310,12 +330,17 @@
         if self.options["linear_solver"] == 'DENSE':
             #Specify the jacobian to the solver
             if self.pData.JAC != NULL and self.options["usejac"]:
-                
-                flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, ida_jac)
+                ELSE:
+                    flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac)
                 if flag < 0:
                     raise IDAError(flag,self.t)
             else:
-                flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, NULL)
+                ELSE:
+                    flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL)
                 if flag < 0:
                     raise IDAError(flag,self.t)
                     
@@ -322,11 +347,17 @@
         elif self.options["linear_solver"] == 'SPGMR':
             #Specify the jacobian times vector function
             if self.pData.JACV != NULL and self.options["usejac"]:
-                flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv);
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, SUNDIALS.ida_spils_jtsetup_dummy, ida_jacv);
+                ELSE:
+                    flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv);
                 if flag < 0:
                     raise IDAError(flag, self.t)
             else:
-                flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL);
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, NULL, NULL);
+                ELSE:
+                    flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL);
                 if flag < 0:
                     raise IDAError(flag, self.t)
         else:
@@ -1415,6 +1446,8 @@
     cdef public N.ndarray yS0
     #cdef N.ndarray _event_info
     cdef public N.ndarray g_old
+    cdef SUNDIALS.SUNMatrix sun_matrix
+    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
     
     def __init__(self, problem):
         Explicit_ODE.__init__(self, problem) #Calls the base class
@@ -1479,6 +1512,13 @@
         if self.cvode_mem != NULL:
             #Free Memory
             SUNDIALS.CVodeFree(&self.cvode_mem)
+        
+        IF SUNDIALS_VERSION >= (3,0,0):
+            if self.sun_matrix != NULL:
+                SUNDIALS.SUNMatDestroy(self.sun_matrix)
+                
+            if self.sun_linearsolver != NULL:
+                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
     
     cpdef get_local_errors(self):
         """
@@ -2019,29 +2059,51 @@
         Updates the simulation options.
         """
         cdef flag
-        
+
         #Choose a linear solver if and only if NEWTON is choosen
         if self.options["linear_solver"] == 'DENSE' and self.options["iter"] == "Newton":
-            #Specify the use of the internal dense linear algebra functions.
-            flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim)
-            if flag < 0:
-                raise CVodeError(flag)
+            IF SUNDIALS_VERSION >= (3,0,0):
+                #Create a dense Sundials matrix
+                self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
+                #Create a dense Sundials linear solver
+                self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix)
+                #Attach it to CVode
+                flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix);
+                if flag < 0:
+                    raise CVodeError(flag)
+            ELSE:
+                #Specify the use of the internal dense linear algebra functions.
+                flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim)
+                if flag < 0:
+                    raise CVodeError(flag)
                 
             #Specify the jacobian to the solver
             if self.pData.JAC != NULL and self.options["usejac"]:
-                flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac);
+                ELSE:
+                    flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac)
                 if flag < 0:
                     raise CVodeError(flag)
             else:
-                flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, NULL);
+                ELSE:
+                    flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL)
                 if flag < 0:
                     raise CVodeError(flag)
                     
         elif self.options["linear_solver"] == 'SPGMR' and self.options["iter"] == "Newton":
-            #Specify the use of CVSPGMR linear solver.
-            flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"])
+            IF SUNDIALS_VERSION >= (3,0,0):
+                #Create the linear solver
+                self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], self.options["maxkrylov"])
+                #Attach it to CVode
+                flag = SUNDIALS.CVSpilsSetLinearSolver(self.cvode_mem, self.sun_linearsolver)
+            ELSE:
+                #Specify the use of CVSPGMR linear solver.
+                flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"])
             if flag < 0:
-                raise CVodeError(flag)
+                raise CVodeError(flag) 
                 
             if self.pData.PREC_SOLVE != NULL:
                 if self.pData.PREC_SETUP != NULL: 
@@ -2055,11 +2117,17 @@
                   
             #Specify the jacobian times vector function
             if self.pData.JACV != NULL and self.options["usejac"]:
-                flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv)
-                if flag < 0:
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, SUNDIALS.cv_spils_jtsetup_dummy, cv_jacv);
+                ELSE:
+                    flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv)
+                if flag < 0: 
                     raise CVodeError(flag)
             else:
-                flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, NULL, NULL);
+                ELSE:
+                    flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL)
                 if flag < 0:
                     raise CVodeError(flag)
         elif self.options["linear_solver"] == 'SPARSE' and self.options["iter"] == "Newton":
@@ -2066,17 +2134,28 @@
             
             if SUNDIALS.version() < (2,6,0): 
                 raise AssimuloException("Not supported with this SUNDIALS version.")
-            
+            if SUNDIALS.with_superlu() == 0:
+                raise AssimuloException("No support for SuperLU was detected, please verify that SuperLU and SUNDIALS has been installed correctly.")
+                
             #Specify the use of CVSPGMR linear solver.
             if self.problem_info["jac_fcn_nnz"] == -1:
                 raise AssimuloException("Need to specify the number of non zero elements in the Jacobian via the option 'jac_nnz'")
-            flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"])
+                
+            IF SUNDIALS_VERSION >= (3,0,0):
+                self.sun_matrix = SUNDIALS.SUNSparseMatrix(self.pData.dim, self.pData.dim, self.problem_info["jac_fcn_nnz"], CSC_MAT)
+                self.sun_linearsolver = SUNDIALS.SUNSuperLUMT(self.yTemp, self.sun_matrix, self.options["num_threads"])
+                flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix)
+            ELSE:
+                flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"])
             if flag < 0:
                     raise CVodeError(flag)
             
             #Specify the jacobian to the solver
             if self.pData.JAC != NULL and self.options["usejac"]:
-                flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse)
+                IF SUNDIALS_VERSION >= (3,0,0):
+                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac_sparse)
+                ELSE:
+                    flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse)
                 if flag < 0:
                     raise CVodeError(flag)
             else:
@@ -2899,7 +2978,10 @@
             flag = SUNDIALS.CVSpilsGetNumRhsEvals(self.cvode_mem, &nfevalsLS) #Number of rhs due to jac*vector
             self.statistics["njacvecs"]  += njvevals
         elif self.options["linear_solver"] == "SPARSE":
-            flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals)
+            IF SUNDIALS_VERSION >= (3,0,0):
+                flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals)
+            ELSE:
+                flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals)
             self.statistics["njacs"]   += njevals
         else:
             flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals) #Number of jac evals
Index: assimulo/lib/sundials_includes.pxd
===================================================================
--- assimulo/lib/sundials_includes.pxd	(revision 844)
+++ assimulo/lib/sundials_includes.pxd	(revision 845)
@@ -37,19 +37,6 @@
     ctypedef double realtype
     ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
 
-#==============================================
-# C headers
-#==============================================
-cdef extern from "string.h":
-    void *memcpy(void *s1, void *s2, int n)
-cdef extern from "stdlib.h":
-    void *malloc(int size)
-    void free(void *ptr)
-    
-#==============================================
-#External definitions from Sundials headers
-#==============================================
-
 cdef extern from "sundials/sundials_nvector.h":
     ctypedef _generic_N_Vector* N_Vector
     
@@ -77,6 +64,107 @@
     void N_VDestroy_Serial(N_Vector v)
     void N_VPrint_Serial(N_Vector v)
 
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef extern from "sundials/sundials_types.h":
+        IF SUNDIALS_VECTOR_SIZE == "64":
+            ctypedef long int sunindextype
+        ELSE:
+            ctypedef int sunindextype
+    cdef extern from "sundials/sundials_matrix.h":
+        ctypedef _generic_SUNMatrix *SUNMatrix;
+        void SUNMatDestroy(SUNMatrix A);
+        
+        cdef struct _generic_SUNMatrix_Ops:
+            SUNMatrix_ID (*getid)(SUNMatrix);
+            SUNMatrix    (*clone)(SUNMatrix);
+            void         (*destroy)(SUNMatrix);
+            int          (*zero)(SUNMatrix);
+            int          (*copy)(SUNMatrix, SUNMatrix);
+            int          (*scaleadd)(realtype, SUNMatrix, SUNMatrix);
+            int          (*scaleaddi)(realtype, SUNMatrix);
+            int          (*matvec)(SUNMatrix, N_Vector, N_Vector);
+            int          (*space)(SUNMatrix, long int*, long int*);
+
+        cdef struct _generic_SUNMatrix:
+            void *content;
+            _generic_SUNMatrix_Ops *ops;
+            
+        cdef enum SUNMatrix_ID:
+            SUNMATRIX_DENSE, 
+            SUNMATRIX_BAND, 
+            SUNMATRIX_SPARSE, 
+            SUNMATRIX_CUSTOM
+    
+    cdef extern from "sundials/sundials_linearsolver.h":
+        ctypedef _generic_SUNLinearSolver *SUNLinearSolver;
+        int SUNLinSolFree(SUNLinearSolver S);
+        
+        cdef struct _generic_SUNLinearSolver_Ops:
+            SUNLinearSolver_Type (*gettype)(SUNLinearSolver);
+            int                  (*setatimes)(SUNLinearSolver, void*, ATimesFn);
+            int                  (*setpreconditioner)(SUNLinearSolver, void*, 
+                                                    PSetupFn, PSolveFn);
+            int                  (*setscalingvectors)(SUNLinearSolver,
+                                                    N_Vector, N_Vector);
+            int                  (*initialize)(SUNLinearSolver);
+            int                  (*setup)(SUNLinearSolver, SUNMatrix);
+            int                  (*solve)(SUNLinearSolver, SUNMatrix, N_Vector, 
+                                        N_Vector, realtype);
+            int                  (*numiters)(SUNLinearSolver);
+            realtype             (*resnorm)(SUNLinearSolver);
+            long int             (*lastflag)(SUNLinearSolver);
+            int                  (*space)(SUNLinearSolver, long int*, long int*);
+            N_Vector             (*resid)(SUNLinearSolver);
+            int                  (*free)(SUNLinearSolver);
+        
+        cdef struct _generic_SUNLinearSolver:
+            void *content;
+            _generic_SUNLinearSolver_Ops *ops;
+            
+        cdef enum SUNLinearSolver_Type:
+            SUNLINEARSOLVER_DIRECT,
+            SUNLINEARSOLVER_ITERATIVE,
+            SUNLINEARSOLVER_CUSTOM
+    
+    cdef extern from "sunmatrix/sunmatrix_dense.h":
+        ctypedef _SUNMatrixContent_Dense *SUNMatrixContent_Dense;
+        cdef struct _SUNMatrixContent_Dense:
+            sunindextype M;
+            sunindextype N;
+            realtype *data;
+            sunindextype ldata;
+            realtype **cols;
+        SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N);
+    cdef extern from "sunmatrix/sunmatrix_sparse.h":
+        ctypedef _SUNMatrixContent_Sparse *SUNMatrixContent_Sparse;
+        cdef struct _SUNMatrixContent_Sparse:
+            sunindextype M;
+            sunindextype N;
+            sunindextype NNZ;
+            sunindextype NP;
+            realtype *data;
+            int sparsetype;
+            sunindextype *indexvals;
+            sunindextype *indexptrs;
+            sunindextype **rowvals;
+            sunindextype **colptrs;
+            sunindextype **colvals;
+            sunindextype **rowptrs;
+        SUNMatrix SUNSparseMatrix(sunindextype M, sunindextype N, sunindextype NNZ, int sparsetype);
+    cdef extern from "sunlinsol/sunlinsol_dense.h":
+        SUNLinearSolver SUNDenseLinearSolver(N_Vector y, SUNMatrix A);
+    cdef extern from "sunlinsol/sunlinsol_spgmr.h":
+        SUNLinearSolver SUNSPGMR(N_Vector y, int pretype, int maxl);
+        
+ELSE: 
+    #Dummy defines
+    ctypedef void *SUNLinearSolver
+    ctypedef void *SUNMatrix
+    ctypedef void *SUNMatrixContent_Dense
+    ctypedef void *SUNMatrixContent_Sparse
+    ctypedef int sunindextype
+
+
 #Struct for handling the Jacobian data
 cdef extern from "sundials/sundials_direct.h":
     cdef struct _DlsMat:
@@ -92,7 +180,7 @@
         realtype **cols
     ctypedef _DlsMat* DlsMat
     cdef realtype* DENSE_COL(DlsMat A, int j)
-
+    
 IF SUNDIALS_VERSION >= (2,6,3):
     cdef extern from "sundials/sundials_sparse.h":
         cdef struct _SlsMat:
@@ -128,7 +216,26 @@
         int *rowvals
         int *colptrs
     ctypedef _SlsMat* SlsMat
+    
+#==============================================
+# C headers
+#==============================================
+cdef extern from "string.h":
+    void *memcpy(void *s1, void *s2, int n)
+cdef extern from "stdlib.h":
+    void *malloc(int size)
+    void free(void *ptr)
+    
+#==============================================
+#External definitions from Sundials headers
+#==============================================
 
+IF SUNDIALS_WITH_SUPERLU:
+    cdef inline int with_superlu(): return 1
+ELSE:
+    cdef inline int with_superlu(): return 0
+
+
 cdef extern from "cvodes/cvodes.h":
     void* CVodeCreate(int lmm, int iter)
     ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
@@ -225,45 +332,67 @@
     int CVodeGetSensNonlinSolvStats(void *cvode_mem, long int *nSniters, long int *nSncfails)
     int CVodeGetStgrSensNumNonlinSolvIters(void *cvode_mem, long int *nSTGR1niters)
     int CVodeGetStgrSensNumNonlinSolvConvFails(void *cvode_mem, long int *nSTGR1ncfails)
+
+cdef extern from "cvodes/cvodes_spils.h":
+    ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
+            N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
+
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef extern from "cvodes/cvodes_direct.h":
+        ctypedef int (*CVDlsDenseJacFn)(realtype t, N_Vector y, N_Vector fy, 
+                       SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+        int CVDlsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS, SUNMatrix A);
+        int CVDlsSetJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
+    cdef extern from "cvodes/cvodes_spils.h":
+        int CVSpilsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS);
+        ctypedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void *user_data);
+        int CVSpilsSetJacTimes(void *cvode_mem, CVSpilsJacTimesSetupFn jtsetup, CVSpilsJacTimesVecFn jtimes);
     
     
-cdef extern from "cvodes/cvodes_dense.h":
-    int CVDense(void *cvode_mem, long int n)
-    ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy, 
-                   DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-    int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
+    IF SUNDIALS_WITH_SUPERLU:
+        cdef extern from "sunlinsol/sunlinsol_superlumt.h":
+            SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
+    ELSE:
+        cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
+    
+    cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0    
+    cdef inline tuple version(): return (3,0,0)
+ELSE:
+    cdef extern from "cvodes/cvodes_dense.h":
+        int CVDense(void *cvode_mem, long int n)
+        ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy, 
+                       DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+        int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
 
-cdef extern from "cvodes/cvodes_spgmr.h":
-    int CVSpgmr(void *cvode_mem, int pretype, int max1)
+    cdef extern from "cvodes/cvodes_spgmr.h":
+        int CVSpgmr(void *cvode_mem, int pretype, int max1)
     
-IF SUNDIALS_VERSION >= (2,6,0):
-    cdef extern from "cvodes/cvodes_sparse.h":
+    cdef extern from "cvodes/cvodes_spils.h":
+        int CVSpilsSetJacTimesVecFn(void *cvode_mem,  CVSpilsJacTimesVecFn jtv)
+    
+    IF SUNDIALS_VERSION >= (2,6,0):
+        cdef extern from "cvodes/cvodes_sparse.h":
+            ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
+                                      SlsMat Jac, void *user_data, N_Vector tmp1,
+                                        N_Vector tmp2, N_Vector tmp3)
+            int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
+            int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
+        cdef inline tuple version(): return (2,6,0)
+        IF SUNDIALS_WITH_SUPERLU:
+            cdef extern from "cvodes/cvodes_sparse.h":
+                int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
+        ELSE:
+            cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
+    ELSE:
+        cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
         ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
                                   SlsMat Jac, void *user_data, N_Vector tmp1,
                                     N_Vector tmp2, N_Vector tmp3)
-        int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
-        int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
-    #cdef inline char* version(): return "2.6.0"
-    cdef inline tuple version(): return (2,6,0)
-    IF SUNDIALS_WITH_SUPERLU:
-        cdef extern from "cvodes/cvodes_sparse.h":
-            int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
-    ELSE:
-        cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-ELSE:
-    cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-    ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
-                              SlsMat Jac, void *user_data, N_Vector tmp1,
-                                N_Vector tmp2, N_Vector tmp3)
-    cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
-    cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
-    #cdef inline char* version(): return "2.5.0"
-    cdef inline tuple version(): return (2,5,0)
+        cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
+        cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
+        cdef inline tuple version(): return (2,5,0)
     
 cdef extern from "cvodes/cvodes_spils.h":
-    ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-				    N_Vector y, N_Vector fy,
-				    void *user_data, N_Vector tmp)
     ctypedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy,
 				  booleantype jok, booleantype *jcurPtr,
 				  realtype gamma, void *user_data,
@@ -273,13 +402,12 @@
 				  N_Vector r, N_Vector z,
 				  realtype gamma, realtype delta,
 				  int lr, void *user_data, N_Vector tmp)
-    int CVSpilsSetJacTimesVecFn(void *cvode_mem,  CVSpilsJacTimesVecFn jtv)
+    
     int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
     int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
     int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
     int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
     int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
-    
 
 cdef extern from "idas/idas.h":
     ctypedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
@@ -380,19 +508,40 @@
     
     #End Sensitivities
     #=================
-
-cdef extern from "idas/idas_dense.h":
-    int IDADense(void *ida_mem, long int n)
-    ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy, 
-                   N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data, 
-                   N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-    int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
     
-cdef extern from "idas/idas_dense.h":
-    int IDASpgmr(void *ida_mem, int max1)
-    ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-    int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
+cdef extern from "idas/idas_spils.h":
+    ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, 
+            N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
     
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef extern from "idas/idas_direct.h":
+        ctypedef int (*IDADlsDenseJacFn)(realtype tt, realtype cj, N_Vector yy, 
+                       N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data, 
+                       N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+        int IDADlsSetJacFn(void *ida_mem, IDADlsDenseJacFn djac)
+        int IDADlsSetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A);
+    
+    cdef extern from "idas/idas_spils.h":
+        int IDASpilsSetLinearSolver(void *ida_mem, SUNLinearSolver LS);
+        ctypedef int (*IDASpilsJacTimesSetupFn)(realtype tt, N_Vector yy,
+                    N_Vector yp, N_Vector rr, realtype c_j, void *user_data);
+        int IDASpilsSetJacTimes(void *ida_mem,
+                IDASpilsJacTimesSetupFn jtsetup, IDASpilsJacTimesVecFn jtimes);
+                
+    cdef inline int ida_spils_jtsetup_dummy(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, realtype c_j, void *user_data): return 0
+ELSE:
+    cdef extern from "idas/idas_dense.h":
+        int IDADense(void *ida_mem, long int n)
+        ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy, 
+                       N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data, 
+                       N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+        int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
+    
+        int IDASpgmr(void *ida_mem, int max1)
+        
+    cdef extern from "idas/idas_spils.h":
+        int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
+
 cdef extern from "idas/idas_spils.h":
     int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
     int IDASpilsGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
@@ -455,14 +604,31 @@
     # fuction used to deallocate memory used by KINSOL
     void KINFree(void **kinmem)
 
-# functions used for supplying jacobian, and receiving info from linear solver
-cdef extern from "kinsol/kinsol_direct.h":
-    # user functions
-    ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
+
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef extern from "kinsol/kinsol_direct.h":
+        ctypedef int (*KINDlsDenseJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
+        int KINDlsSetLinearSolver(void *kinmem, SUNLinearSolver LS, SUNMatrix A);
+        int KINDlsSetJacFn(void *kinmem, KINDlsDenseJacFn djac)
     
-    # function used to link user functions to KINSOL
-    int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
+    cdef extern from "kinsol/kinsol_spils.h":
+        int KINSpilsSetLinearSolver(void *kinsol_mem, SUNLinearSolver LS);
+ELSE:
+    # functions used for supplying jacobian, and receiving info from linear solver
+    cdef extern from "kinsol/kinsol_direct.h":
+        # user functions
+        ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
+        
+        # function used to link user functions to KINSOL
+        int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
+    
+    cdef extern from "kinsol/kinsol_dense.h":
+        int KINDense(void *kinmem, int dim)
+    
+    cdef extern from "kinsol/kinsol_spgmr.h":
+        int KINSpgmr(void *kinmem, int maxl)
 
+cdef extern from "kinsol/kinsol_direct.h":
     # optional output fcts for linear direct solver
     int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
     int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB)
@@ -470,12 +636,6 @@
     int KINDlsGetLastFlag(void *kinmem, int *flag)
     char *KINDlsGetReturnFlagName(int flag)
 
-cdef extern from "kinsol/kinsol_dense.h":
-    int KINDense(void *kinmem, int dim)
-    
-cdef extern from "kinsol/kinsol_spgmr.h":
-    int KINSpgmr(void *kinmem, int maxl)
-
 cdef extern from "kinsol/kinsol_spils.h":
     ctypedef int (*KINSpilsJacTimesVecFn)(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
                 void *problem_data)
Index: src/lib/sundials_constants.pxi
===================================================================
--- assimulo/lib/sundials_constants.pxi	(revision 844)
+++ assimulo/lib/sundials_constants.pxi	(revision 845)
@@ -95,6 +95,9 @@
 DEF CV_BAD_IS             = -45
 
 
+DEF CSC_MAT = 0
+DEF CSR_MAT = 1
+
 #==========
 # IDA
 #==========
Index: src/lib/sundials_callbacks.pxi
===================================================================
--- assimulo/lib/sundials_callbacks.pxi	(revision 844)
+++ assimulo/lib/sundials_callbacks.pxi	(revision 845)
@@ -94,98 +94,55 @@
         traceback.print_exc()
         return CV_UNREC_RHSFUNC_ERR 
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
-                void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-    """
-    This method is used to connect the Assimulo.Problem.jac to the Sundials
-    Sparse Jacobian function.
-    """
-    cdef ProblemData pData = <ProblemData>problem_data
-    #cdef N.ndarray y = nv2arr(yv)
-    cdef N.ndarray y = pData.work_y
-    cdef int i
-    cdef int nnz = Jacobian.NNZ
-    cdef int ret_nnz
-    cdef int dim = Jacobian.N
-    cdef realtype* data = Jacobian.data
-    
-    IF SUNDIALS_VERSION >= (2,6,3):
-        cdef int* rowvals = Jacobian.rowvals[0]
-        cdef int* colptrs = Jacobian.colptrs[0]
-    ELSE:
-        cdef int* rowvals = Jacobian.rowvals
-        cdef int* colptrs = Jacobian.colptrs
-    
-    nv2arr_inplace(yv, y)
-    """
-        realtype *data;
-        int *rowvals;
-        int *colptrs;
-    """
-    try:
-        if pData.dimSens > 0: #Sensitivity activated
-            p = realtype2arr(pData.p,pData.dimSens)
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-            else:
-                jac=(<object>pData.JAC)(t,y,p=p)
-        else:
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-            else:
-                jac=(<object>pData.JAC)(t,y)
-            
-        if not isinstance(jac, sparse.csc.csc_matrix):
-            jac = sparse.csc.csc_matrix(jac)
-            raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-        ret_nnz = jac.nnz
-        if ret_nnz > nnz:
-            raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")    
 
-        for i in range(min(ret_nnz,nnz)):
-            data[i]    = jac.data[i]
-            rowvals[i] = jac.indices[i]
-        for i in range(dim+1):
-            colptrs[i] = jac.indptr[i]
+IF SUNDIALS_VERSION >= (3,0,0):
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
+                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Sparse Jacobian function.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef SUNMatrixContent_Sparse Jacobian = <SUNMatrixContent_Sparse>Jac.content
+        cdef N.ndarray y = pData.work_y
+        cdef int i
+        cdef sunindextype nnz = Jacobian.NNZ
+        cdef int ret_nnz
+        cdef sunindextype dim = Jacobian.N
+        cdef realtype* data = Jacobian.data
+        cdef sunindextype* rowvals = Jacobian.rowvals[0]
+        cdef sunindextype* colptrs = Jacobian.colptrs[0]
         
-        return CVDLS_SUCCESS
-    except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-        return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-    except:
-        traceback.print_exc()
-        return CVDLS_JACFUNC_UNRECVR
+        nv2arr_inplace(yv, y)
 
-cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian, 
-                void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-    """
-    This method is used to connect the Assimulo.Problem.jac to the Sundials
-    Jacobian function.
-    """
-    cdef ProblemData pData = <ProblemData>problem_data
-    #cdef ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-    cdef realtype* col_i=DENSE_COL(Jacobian,0)
-    #(<ndarray>pData.y).data =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
-    #cdef N.ndarray y = nv2arr(yv)
-    cdef N.ndarray y = pData.work_y
-    cdef int i,j
-    
-    nv2arr_inplace(yv, y)
-
-    if pData.dimSens>0: #Sensitivity activated
-        p = realtype2arr(pData.p,pData.dimSens)
         try:
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
+            if pData.dimSens > 0: #Sensitivity activated
+                p = realtype2arr(pData.p,pData.dimSens)
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y,p=p)
             else:
-                jac=(<object>pData.JAC)(t,y,p)
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y)
                 
-            for i in range(Neq):
-                col_i = DENSE_COL(Jacobian, i)
-                for j in range(Neq):
-                    col_i[j] = jac[j,i]
+            if not isinstance(jac, sparse.csc.csc_matrix):
+                jac = sparse.csc.csc_matrix(jac)
+                raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
+            ret_nnz = jac.nnz
+            if ret_nnz > nnz:
+                raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")    
 
+            for i in range(min(ret_nnz,nnz)):
+                data[i]    = jac.data[i]
+                rowvals[i] = jac.indices[i]
+            for i in range(dim+1):
+                colptrs[i] = jac.indptr[i]
+            
             return CVDLS_SUCCESS
         except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
             return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
@@ -192,18 +149,62 @@
         except:
             traceback.print_exc()
             return CVDLS_JACFUNC_UNRECVR
-    else:
+ELSE:
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
+                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Sparse Jacobian function.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef N.ndarray y = pData.work_y
+        cdef int i
+        cdef int nnz = Jacobian.NNZ
+        cdef int ret_nnz
+        cdef int dim = Jacobian.N
+        cdef realtype* data = Jacobian.data
+        
+        IF SUNDIALS_VERSION >= (2,6,3):
+            cdef int* rowvals = Jacobian.rowvals[0]
+            cdef int* colptrs = Jacobian.colptrs[0]
+        ELSE:
+            cdef int* rowvals = Jacobian.rowvals
+            cdef int* colptrs = Jacobian.colptrs
+        
+        nv2arr_inplace(yv, y)
+        """
+            realtype *data;
+            int *rowvals;
+            int *colptrs;
+        """
         try:
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
+            if pData.dimSens > 0: #Sensitivity activated
+                p = realtype2arr(pData.p,pData.dimSens)
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y,p=p)
             else:
-                jac=(<object>pData.JAC)(t,y)
-    
-            for i in range(Neq):
-                col_i = DENSE_COL(Jacobian, i)
-                for j in range(Neq):
-                    col_i[j] = jac[j,i]
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y)
+                
+            if not isinstance(jac, sparse.csc.csc_matrix):
+                jac = sparse.csc.csc_matrix(jac)
+                raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
+            ret_nnz = jac.nnz
+            if ret_nnz > nnz:
+                raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")    
 
+            for i in range(min(ret_nnz,nnz)):
+                data[i]    = jac.data[i]
+                rowvals[i] = jac.indices[i]
+            for i in range(dim+1):
+                colptrs[i] = jac.indptr[i]
+            
             return CVDLS_SUCCESS
         except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
             return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
@@ -210,7 +211,112 @@
         except:
             traceback.print_exc()
             return CVDLS_JACFUNC_UNRECVR
+
+
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef int cv_jac(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac, 
+                void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef realtype* col_i=Jacobian.cols[0]
+        cdef N.ndarray y = pData.work_y
+        cdef int i,j, Neq = pData.dim
         
+        nv2arr_inplace(yv, y)
+
+        if pData.dimSens>0: #Sensitivity activated
+            p = realtype2arr(pData.p,pData.dimSens)
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
+                else:
+                    jac=(<object>pData.JAC)(t,y,p)
+                    
+                for i in range(Neq):
+                    col_i = Jacobian.cols[i]
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+
+                return CVDLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+            except:
+                traceback.print_exc()
+                return CVDLS_JACFUNC_UNRECVR
+        else:
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y)
+        
+                for i in range(Neq):
+                    col_i = Jacobian.cols[i]
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+
+                return CVDLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+            except:
+                traceback.print_exc()
+                return CVDLS_JACFUNC_UNRECVR
+ELSE:
+    cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian, 
+                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef realtype* col_i=DENSE_COL(Jacobian,0)
+        cdef N.ndarray y = pData.work_y
+        cdef int i,j
+        
+        nv2arr_inplace(yv, y)
+
+        if pData.dimSens>0: #Sensitivity activated
+            p = realtype2arr(pData.p,pData.dimSens)
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
+                else:
+                    jac=(<object>pData.JAC)(t,y,p)
+                    
+                for i in range(Neq):
+                    col_i = DENSE_COL(Jacobian, i)
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+
+                return CVDLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+            except:
+                traceback.print_exc()
+                return CVDLS_JACFUNC_UNRECVR
+        else:
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
+                else:
+                    jac=(<object>pData.JAC)(t,y)
+        
+                for i in range(Neq):
+                    col_i = DENSE_COL(Jacobian, i)
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+
+                return CVDLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+            except:
+                traceback.print_exc()
+                return CVDLS_JACFUNC_UNRECVR
+        
 cdef int cv_jacv(N_Vector vv, N_Vector Jv, realtype t, N_Vector yv, N_Vector fyv,
 				    void *problem_data, N_Vector tmp):
     """
@@ -337,9 +443,6 @@
     Root-finding function.
     """
     cdef ProblemData pData = <ProblemData>problem_data
-    #cdef ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-    #(<ndarray>pData.y).data =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
-    #cdef N.ndarray y = nv2arr(yv)
     cdef N.ndarray y = pData.work_y
     cdef int i
     
@@ -366,12 +469,8 @@
     """
     cdef ProblemData pData = <ProblemData>problem_data
     cdef N.ndarray[realtype, ndim=1, mode='c'] res #Used for return from the user function
-    #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
-    #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
     cdef N.ndarray y = pData.work_y
     cdef N.ndarray yd = pData.work_yd
-    # cdef N.ndarray y = nv2arr(yv)
-    # cdef N.ndarray yd = nv2arr(yvdot)
     cdef realtype* resptr=(<N_VectorContent_Serial>residual.content).data
     cdef int i
     
@@ -414,63 +513,113 @@
         except:
             traceback.print_exc()
             return IDA_RES_FAIL
-            
-cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
+
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef int ida_jac(realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, SUNMatrix Jac,
+                 void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
+        cdef realtype* col_i=Jacobian.cols[0]
+        cdef N.ndarray y = pData.work_y
+        cdef N.ndarray yd = pData.work_yd
+        cdef int i,j, Neq = pData.dim
+        
+        nv2arr_inplace(yv, y)
+        nv2arr_inplace(yvdot, yd)
+        
+        if pData.dimSens!=0: #SENSITIVITY 
+            p = realtype2arr(pData.p,pData.dimSens)
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
+                else:
+                    jac=(<object>pData.JAC)(c,t,y,yd,p=p)
+                
+                for i in range(Neq):
+                    col_i = Jacobian.cols[i]
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+                return IDADLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return IDADLS_JACFUNC_RECVR #Recoverable Error
+            except:
+                traceback.print_exc()
+                return IDADLS_JACFUNC_UNRECVR
+        else:
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
+                else:
+                    jac=(<object>pData.JAC)(c,t,y,yd)
+                
+                for i in range(Neq):
+                    col_i = Jacobian.cols[i]
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+                return IDADLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return IDADLS_JACFUNC_RECVR #Recoverable Error
+            except:
+                traceback.print_exc()
+                return IDADLS_JACFUNC_UNRECVR
+ELSE:
+    cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
                  void* problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-    """
-    This method is used to connect the Assimulo.Problem.jac to the Sundials
-    Jacobian function.
-    """
-    cdef ProblemData pData = <ProblemData>problem_data
-    cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-    cdef realtype* col_i=DENSE_COL(Jacobian,0)
-    #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
-    #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
-    cdef N.ndarray y = pData.work_y
-    cdef N.ndarray yd = pData.work_yd
-    #cdef N.ndarray y = nv2arr(yv)
-    #cdef N.ndarray yd = nv2arr(yvdot)
-    cdef int i,j
-    
-    nv2arr_inplace(yv, y)
-    nv2arr_inplace(yvdot, yd)
-    
-    if pData.dimSens!=0: #SENSITIVITY 
-        p = realtype2arr(pData.p,pData.dimSens)
-        try:
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
-            else:
-                jac=(<object>pData.JAC)(c,t,y,yd,p=p)
+        """
+        This method is used to connect the Assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
+        cdef realtype* col_i=DENSE_COL(Jacobian,0)
+        cdef N.ndarray y = pData.work_y
+        cdef N.ndarray yd = pData.work_yd
+        cdef int i,j
+        
+        nv2arr_inplace(yv, y)
+        nv2arr_inplace(yvdot, yd)
+        
+        if pData.dimSens!=0: #SENSITIVITY 
+            p = realtype2arr(pData.p,pData.dimSens)
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
+                else:
+                    jac=(<object>pData.JAC)(c,t,y,yd,p=p)
+                
+                for i in range(Neq):
+                    col_i = DENSE_COL(Jacobian, i)
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+                return IDADLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return IDADLS_JACFUNC_RECVR #Recoverable Error
+            except:
+                traceback.print_exc()
+                return IDADLS_JACFUNC_UNRECVR
+        else:
+            try:
+                if pData.sw != NULL:
+                    jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
+                else:
+                    jac=(<object>pData.JAC)(c,t,y,yd)
+                
+                for i in range(Neq):
+                    col_i = DENSE_COL(Jacobian, i)
+                    for j in range(Neq):
+                        col_i[j] = jac[j,i]
+                return IDADLS_SUCCESS
+            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
+                return IDADLS_JACFUNC_RECVR #Recoverable Error
+            except:
+                traceback.print_exc()
+                return IDADLS_JACFUNC_UNRECVR
             
-            for i in range(Neq):
-                col_i = DENSE_COL(Jacobian, i)
-                for j in range(Neq):
-                    col_i[j] = jac[j,i]
-            return IDADLS_SUCCESS
-        except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-            return IDADLS_JACFUNC_RECVR #Recoverable Error
-        except:
-            traceback.print_exc()
-            return IDADLS_JACFUNC_UNRECVR
-    else:
-        try:
-            if pData.sw != NULL:
-                jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
-            else:
-                jac=(<object>pData.JAC)(c,t,y,yd)
-            
-            for i in range(Neq):
-                col_i = DENSE_COL(Jacobian, i)
-                for j in range(Neq):
-                    col_i[j] = jac[j,i]
-            return IDADLS_SUCCESS
-        except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-            return IDADLS_JACFUNC_RECVR #Recoverable Error
-        except:
-            traceback.print_exc()
-            return IDADLS_JACFUNC_UNRECVR
-        
 
 cdef int ida_root(realtype t, N_Vector yv, N_Vector yvdot, realtype *gout,  void* problem_data):
     """
@@ -479,12 +628,8 @@
     """
     cdef ProblemData pData = <ProblemData>problem_data
     cdef N.ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-    #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
-    #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
     cdef N.ndarray y = pData.work_y
     cdef N.ndarray yd = pData.work_yd
-    #cdef N.ndarray y = nv2arr(yv)
-    #cdef N.ndarray yd = nv2arr(yvdot)
     cdef int i
     
     nv2arr_inplace(yv, y)
@@ -553,30 +698,54 @@
             traceback.print_exc()
             return SPGMR_PSOLVE_FAIL_UNREC
     
+IF SUNDIALS_VERSION >= (3,0,0):
+    cdef int kin_jac(N_Vector xv, N_Vector fval, SUNMatrix Jac, 
+                    void *problem_data, N_Vector tmp1, N_Vector tmp2):
+        """
+        This method is used to connect the assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
+        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
+        cdef realtype* col_i=Jacobian.cols[0]
+        cdef N.ndarray x = nv2arr(xv)
+        cdef int i,j, Neq = pData.dim
+        
+        try:
+            jac=(<object>pData.JAC)(x)
 
-cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, 
-                void *problem_data, N_Vector tmp1, N_Vector tmp2):
-    """
-    This method is used to connect the assimulo.Problem.jac to the Sundials
-    Jacobian function.
-    """
-    cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-    cdef realtype* col_i=DENSE_COL(Jacobian,0)
-    cdef N.ndarray x = nv2arr(xv)
-    cdef int i,j
-    
-    try:
-        jac=(<object>pData.JAC)(x)
+            for i in range(Neq):
+                col_i = Jacobian.cols[i]
+                for j in range(Neq):
+                    col_i[j] = jac[j,i]
 
-        for i in range(Neq):
-            col_i = DENSE_COL(Jacobian, i)
-            for j in range(Neq):
-                col_i[j] = jac[j,i]
+            return KINDLS_SUCCESS
+        except:
+            return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+ELSE:
+    cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, 
+                    void *problem_data, N_Vector tmp1, N_Vector tmp2):
+        """
+        This method is used to connect the assimulo.Problem.jac to the Sundials
+        Jacobian function.
+        """
+        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
+        cdef realtype* col_i=DENSE_COL(Jacobian,0)
+        cdef N.ndarray x = nv2arr(xv)
+        cdef int i,j
+        
+        try:
+            jac=(<object>pData.JAC)(x)
 
-        return KINDLS_SUCCESS
-    except:
-        return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-        
+            for i in range(Neq):
+                col_i = DENSE_COL(Jacobian, i)
+                for j in range(Neq):
+                    col_i[j] = jac[j,i]
+
+            return KINDLS_SUCCESS
+        except:
+            return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
+            
 cdef int kin_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
             void *problem_data):
     cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
Index: examples/cvode_with_jac.py
===================================================================
--- assimulo/examples/cvode_with_jac.py	(revision 844)
+++ assimulo/examples/cvode_with_jac.py	(revision 845)
@@ -46,7 +46,6 @@
     def f(t,y):
         yd_0 = y[1]
         yd_1 = -9.82
-        #print y, yd_0, yd_1
         return N.array([yd_0,yd_1])
     
     #Defines the Jacobian
Index: examples/__init__.py
===================================================================
--- assimulo/examples/__init__.py	(revision 844)
+++ assimulo/examples/__init__.py	(revision 845)
@@ -24,6 +24,6 @@
            "mech_system_pendulum", "euler_vanderpol", "cvode_with_parameters_modified",
            "cvode_basic_backward","ida_basic_backward","dasp3_basic", "cvode_with_preconditioning",
            "kinsol_basic","kinsol_with_jac", "radau5dae_time_events", "kinsol_ors", "lsodar_bouncing_ball",
-           "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result"]
+           "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result", "cvode_with_jac_sparse"]
 
 
Index: examples/cvode_with_jac_sparse.py
===================================================================
--- assimulo/examples/cvode_with_jac_sparse.py	(nonexistent)
+++ assimulo/examples/cvode_with_jac_sparse.py	(revision 845)
@@ -0,0 +1,102 @@
+#!/usr/bin/env python 
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2010 Modelon AB
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as N
+import scipy.sparse as SP
+import pylab as P
+import nose
+from assimulo.solvers import CVode
+from assimulo.problem import Explicit_Problem
+
+
+def run_example(with_plots=True):
+    r"""
+    Example for demonstrating the use of a user supplied Jacobian (sparse).
+    Note that this will only work if Assimulo has been configured with
+    Sundials + SuperLU. Based on the SUNDIALS example cvRoberts_sps.c
+    
+    ODE:
+    
+    .. math::
+       
+       \dot y_1 &= -0.04y_1 + 1e4 y_2 y_3 \\
+       \dot y_2 &= - \dot y_1 - \dot y_3 \\
+       \dot y_3 &= 3e7 y_2^2
+       
+    
+    on return:
+    
+       - :dfn:`exp_mod`    problem instance
+    
+       - :dfn:`exp_sim`    solver instance
+       
+    """
+    
+    #Defines the rhs
+    def f(t,y):
+        yd_0 = -0.04*y[0] + 1e4*y[1]*y[2]
+        yd_2 = 3e7*y[1]*y[1]
+        yd_1 = -yd_0 - yd_2
+        return N.array([yd_0,yd_1,yd_2])
+    
+    #Defines the Jacobian
+    def jac(t,y):
+        
+        colptrs = [0,3,6,9]
+        rowvals = [0, 1, 2, 0, 1, 2, 0, 1, 2]
+        data = [-0.04, 0.04, 0.0, 1e4*y[2], -1e4*y[2]-6e7*y[1], 6e7*y[1], 1e4*y[1], -1e4*y[1], 0.0]
+
+        J = SP.csc_matrix((data, rowvals, colptrs))
+        return J
+    
+    #Defines an Assimulo explicit problem
+    y0 = [1.0,0.0,0.0] #Initial conditions
+
+    exp_mod = Explicit_Problem(f,y0, name = 'Example using analytic (sparse) Jacobian')
+    
+    exp_mod.jac = jac #Sets the Jacobian
+    exp_mod.jac_nnz = 9
+   
+    
+    exp_sim = CVode(exp_mod) #Create a CVode solver
+    
+    #Set the parameters
+    exp_sim.iter = 'Newton' #Default 'FixedPoint'
+    exp_sim.discr = 'BDF' #Default 'Adams'
+    exp_sim.atol = [1e-8,1e-14,1e-6] #Default 1e-6
+    exp_sim.rtol = 1e-4 #Default 1e-6
+    exp_sim.linear_solver = "sparse"
+    
+    #Simulate
+    t, y = exp_sim.simulate(0.4) #Simulate 0.4 seconds
+    
+    #Basic tests
+    nose.tools.assert_almost_equal(y[-1][0],0.9851,3)
+        
+    #Plot
+    if with_plots:
+        P.plot(t,y[:,1],linestyle="dashed",marker="o") #Plot the solution
+        P.xlabel('Time')
+        P.ylabel('State')
+        P.title(exp_mod.name)
+        P.show()
+        
+    return exp_mod, exp_sim
+
+
+if __name__=='__main__':
+    mod,sim = run_example()
Index: tests/test_examples.py
===================================================================
--- assimulo/tests/test_examples.py	(revision 844)
+++ assimulo/tests/test_examples.py	(revision 845)
@@ -22,7 +22,15 @@
 
 class Test_Examples:
     
+    
     @testattr(stddist = True)
+    def test_cvode_with_jac_sparse(self):
+        try:
+            cvode_with_jac_sparse.run_example(with_plots=False)
+        except AssimuloException:
+            pass #Handle the case when SuperLU is not installed
+    
+    @testattr(stddist = True)
     def test_ida_with_user_defined_handle_result(self):
         ida_with_user_defined_handle_result.run_example(with_plots=False)