Changeset 845


Ignore:
Timestamp:
Dec 22, 2017 10:06:17 AM (8 months ago)
Author:
chria
Message:

Added support for Sundials 3.1. Related to ticket:418

Location:
trunk
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/CHANGELOG

    r844 r845  
    22
    33--- Assimulo-trunk ---
     4    * Added support for Sundials 3.1 (ticket:418)
    45    * Renamed the option stablimit to stablimit (ticket:417)
    56
  • trunk/examples/__init__.py

    r805 r845  
    2525           "cvode_basic_backward","ida_basic_backward","dasp3_basic", "cvode_with_preconditioning",
    2626           "kinsol_basic","kinsol_with_jac", "radau5dae_time_events", "kinsol_ors", "lsodar_bouncing_ball",
    27            "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result"]
     27           "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result", "cvode_with_jac_sparse"]
    2828
    2929
  • trunk/examples/cvode_with_jac.py

    r605 r845  
    4747        yd_0 = y[1]
    4848        yd_1 = -9.82
    49         #print y, yd_0, yd_1
    5049        return N.array([yd_0,yd_1])
    5150   
  • trunk/setup.py

    r843 r845  
    310310            L.debug('SUNDIALS found.')
    311311            sundials_version = None
     312            sundials_vector_type_size = None
    312313           
    313314            try:
     
    315316                    with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f:
    316317                        for line in f:
    317                             if "SUNDIALS_PACKAGE_VERSION" in line:
     318                            if "SUNDIALS_PACKAGE_VERSION" in line or "SUNDIALS_VERSION" in line:
    318319                                sundials_version = tuple([int(f) for f in line.split()[-1][1:-1].split(".")])
    319320                                L.debug('SUNDIALS %d.%d found.'%(sundials_version[0], sundials_version[1]))
    320             except Exception:
     321                                break
     322                    with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f:
     323                        for line in f:
     324                            if "SUNDIALS_INT32_T" in line and line.startswith("#define"):
     325                                sundials_vector_type_size = "32"
     326                                L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size))
     327                                break
     328                            if "SUNDIALS_INT64_T" in line and line.startswith("#define"):
     329                                sundials_vector_type_size = "64"
     330                                L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size))
     331                                if self.with_SLU:
     332                                    L.warning("It is recommended to set the SUNDIALS_INDEX_TYPE to an 32bit integer when using SUNDIALS together with SuperLU.")
     333                                    L.warning("SuperLU may not function properly.")
     334                                break
     335            except Exception as e:
    321336                if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6
    322337                    sundials_version = (2,6,0)
     
    327342               
    328343            self.SUNDIALS_version = sundials_version
     344            self.SUNDIALS_vector_size = sundials_vector_type_size
    329345           
    330346        else:   
     
    376392        if self.with_SUNDIALS:
    377393            compile_time_env = {'SUNDIALS_VERSION': self.SUNDIALS_version,
    378                                 'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu}
     394                                'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu,
     395                                'SUNDIALS_VECTOR_SIZE': self.SUNDIALS_vector_size}
    379396            #CVode and IDA
    380397            ext_list += cythonize(["assimulo" + os.path.sep + "solvers" + os.path.sep + "sundials.pyx"],
     
    383400            ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs]
    384401            ext_list[-1].library_dirs = [self.libdirs]
    385             ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"]
     402           
     403            if self.SUNDIALS_version >= (3,0,0):
     404                ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas", "sundials_sunlinsoldense", "sundials_sunlinsolspgmr", "sundials_sunmatrixdense", "sundials_sunmatrixsparse"]
     405            else:
     406                ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"]
    386407            if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU
     408                if self.SUNDIALS_version >= (3,0,0):
     409                    ext_list[-1].libraries.extend(["sundials_sunlinsolsuperlumt"])
     410               
    387411                ext_list[-1].include_dirs.append(self.SLUincdir)
    388412                ext_list[-1].library_dirs.append(self.SLUlibdir)
    389413                ext_list[-1].libraries.extend(self.superLUFiles)
     414               
    390415       
    391416            #Kinsol
     
    396421            ext_list[-1].library_dirs = [self.libdirs]
    397422            ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"]
    398    
     423           
     424            if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU
     425                ext_list[-1].include_dirs.append(self.SLUincdir)
     426                ext_list[-1].library_dirs.append(self.SLUlibdir)
     427                ext_list[-1].libraries.extend(self.superLUFiles)
    399428       
    400429        for el in ext_list:
  • trunk/src/lib/sundials_callbacks.pxi

    r840 r845  
    9595        return CV_UNREC_RHSFUNC_ERR
    9696
    97 @cython.boundscheck(False)
    98 @cython.wraparound(False)
    99 cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
    100                 void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
    101     """
    102     This method is used to connect the Assimulo.Problem.jac to the Sundials
    103     Sparse Jacobian function.
    104     """
    105     cdef ProblemData pData = <ProblemData>problem_data
    106     #cdef N.ndarray y = nv2arr(yv)
    107     cdef N.ndarray y = pData.work_y
    108     cdef int i
    109     cdef int nnz = Jacobian.NNZ
    110     cdef int ret_nnz
    111     cdef int dim = Jacobian.N
    112     cdef realtype* data = Jacobian.data
    113    
    114     IF SUNDIALS_VERSION >= (2,6,3):
    115         cdef int* rowvals = Jacobian.rowvals[0]
    116         cdef int* colptrs = Jacobian.colptrs[0]
    117     ELSE:
    118         cdef int* rowvals = Jacobian.rowvals
    119         cdef int* colptrs = Jacobian.colptrs
    120    
    121     nv2arr_inplace(yv, y)
    122     """
    123         realtype *data;
    124         int *rowvals;
    125         int *colptrs;
    126     """
    127     try:
    128         if pData.dimSens > 0: #Sensitivity activated
    129             p = realtype2arr(pData.p,pData.dimSens)
    130             if pData.sw != NULL:
    131                 jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
     97
     98IF SUNDIALS_VERSION >= (3,0,0):
     99    @cython.boundscheck(False)
     100    @cython.wraparound(False)
     101    cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
     102                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
     103        """
     104        This method is used to connect the Assimulo.Problem.jac to the Sundials
     105        Sparse Jacobian function.
     106        """
     107        cdef ProblemData pData = <ProblemData>problem_data
     108        cdef SUNMatrixContent_Sparse Jacobian = <SUNMatrixContent_Sparse>Jac.content
     109        cdef N.ndarray y = pData.work_y
     110        cdef int i
     111        cdef sunindextype nnz = Jacobian.NNZ
     112        cdef int ret_nnz
     113        cdef sunindextype dim = Jacobian.N
     114        cdef realtype* data = Jacobian.data
     115        cdef sunindextype* rowvals = Jacobian.rowvals[0]
     116        cdef sunindextype* colptrs = Jacobian.colptrs[0]
     117       
     118        nv2arr_inplace(yv, y)
     119
     120        try:
     121            if pData.dimSens > 0: #Sensitivity activated
     122                p = realtype2arr(pData.p,pData.dimSens)
     123                if pData.sw != NULL:
     124                    jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
     125                else:
     126                    jac=(<object>pData.JAC)(t,y,p=p)
    132127            else:
    133                 jac=(<object>pData.JAC)(t,y,p=p)
    134         else:
    135             if pData.sw != NULL:
    136                 jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
    137             else:
    138                 jac=(<object>pData.JAC)(t,y)
    139            
    140         if not isinstance(jac, sparse.csc.csc_matrix):
    141             jac = sparse.csc.csc_matrix(jac)
    142             raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
    143         ret_nnz = jac.nnz
    144         if ret_nnz > nnz:
    145             raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")   
    146 
    147         for i in range(min(ret_nnz,nnz)):
    148             data[i]    = jac.data[i]
    149             rowvals[i] = jac.indices[i]
    150         for i in range(dim+1):
    151             colptrs[i] = jac.indptr[i]
    152        
    153         return CVDLS_SUCCESS
    154     except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
    155         return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
    156     except:
    157         traceback.print_exc()
    158         return CVDLS_JACFUNC_UNRECVR
    159 
    160 cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
    161                 void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
    162     """
    163     This method is used to connect the Assimulo.Problem.jac to the Sundials
    164     Jacobian function.
    165     """
    166     cdef ProblemData pData = <ProblemData>problem_data
    167     #cdef ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
    168     cdef realtype* col_i=DENSE_COL(Jacobian,0)
    169     #(<ndarray>pData.y).data =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
    170     #cdef N.ndarray y = nv2arr(yv)
    171     cdef N.ndarray y = pData.work_y
    172     cdef int i,j
    173    
    174     nv2arr_inplace(yv, y)
    175 
    176     if pData.dimSens>0: #Sensitivity activated
    177         p = realtype2arr(pData.p,pData.dimSens)
    178         try:
    179             if pData.sw != NULL:
    180                 jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
    181             else:
    182                 jac=(<object>pData.JAC)(t,y,p)
     128                if pData.sw != NULL:
     129                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
     130                else:
     131                    jac=(<object>pData.JAC)(t,y)
    183132               
    184             for i in range(Neq):
    185                 col_i = DENSE_COL(Jacobian, i)
    186                 for j in range(Neq):
    187                     col_i[j] = jac[j,i]
    188 
     133            if not isinstance(jac, sparse.csc.csc_matrix):
     134                jac = sparse.csc.csc_matrix(jac)
     135                raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
     136            ret_nnz = jac.nnz
     137            if ret_nnz > nnz:
     138                raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")   
     139
     140            for i in range(min(ret_nnz,nnz)):
     141                data[i]    = jac.data[i]
     142                rowvals[i] = jac.indices[i]
     143            for i in range(dim+1):
     144                colptrs[i] = jac.indptr[i]
     145           
    189146            return CVDLS_SUCCESS
    190147        except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     
    193150            traceback.print_exc()
    194151            return CVDLS_JACFUNC_UNRECVR
    195     else:
    196         try:
    197             if pData.sw != NULL:
    198                 jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
     152ELSE:
     153    @cython.boundscheck(False)
     154    @cython.wraparound(False)
     155    cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
     156                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
     157        """
     158        This method is used to connect the Assimulo.Problem.jac to the Sundials
     159        Sparse Jacobian function.
     160        """
     161        cdef ProblemData pData = <ProblemData>problem_data
     162        cdef N.ndarray y = pData.work_y
     163        cdef int i
     164        cdef int nnz = Jacobian.NNZ
     165        cdef int ret_nnz
     166        cdef int dim = Jacobian.N
     167        cdef realtype* data = Jacobian.data
     168       
     169        IF SUNDIALS_VERSION >= (2,6,3):
     170            cdef int* rowvals = Jacobian.rowvals[0]
     171            cdef int* colptrs = Jacobian.colptrs[0]
     172        ELSE:
     173            cdef int* rowvals = Jacobian.rowvals
     174            cdef int* colptrs = Jacobian.colptrs
     175       
     176        nv2arr_inplace(yv, y)
     177        """
     178            realtype *data;
     179            int *rowvals;
     180            int *colptrs;
     181        """
     182        try:
     183            if pData.dimSens > 0: #Sensitivity activated
     184                p = realtype2arr(pData.p,pData.dimSens)
     185                if pData.sw != NULL:
     186                    jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
     187                else:
     188                    jac=(<object>pData.JAC)(t,y,p=p)
    199189            else:
    200                 jac=(<object>pData.JAC)(t,y)
    201    
    202             for i in range(Neq):
    203                 col_i = DENSE_COL(Jacobian, i)
    204                 for j in range(Neq):
    205                     col_i[j] = jac[j,i]
    206 
     190                if pData.sw != NULL:
     191                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
     192                else:
     193                    jac=(<object>pData.JAC)(t,y)
     194               
     195            if not isinstance(jac, sparse.csc.csc_matrix):
     196                jac = sparse.csc.csc_matrix(jac)
     197                raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
     198            ret_nnz = jac.nnz
     199            if ret_nnz > nnz:
     200                raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")   
     201
     202            for i in range(min(ret_nnz,nnz)):
     203                data[i]    = jac.data[i]
     204                rowvals[i] = jac.indices[i]
     205            for i in range(dim+1):
     206                colptrs[i] = jac.indptr[i]
     207           
    207208            return CVDLS_SUCCESS
    208209        except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     
    211212            traceback.print_exc()
    212213            return CVDLS_JACFUNC_UNRECVR
     214
     215
     216IF SUNDIALS_VERSION >= (3,0,0):
     217    cdef int cv_jac(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
     218                void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
     219        """
     220        This method is used to connect the Assimulo.Problem.jac to the Sundials
     221        Jacobian function.
     222        """
     223        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
     224        cdef ProblemData pData = <ProblemData>problem_data
     225        cdef realtype* col_i=Jacobian.cols[0]
     226        cdef N.ndarray y = pData.work_y
     227        cdef int i,j, Neq = pData.dim
     228       
     229        nv2arr_inplace(yv, y)
     230
     231        if pData.dimSens>0: #Sensitivity activated
     232            p = realtype2arr(pData.p,pData.dimSens)
     233            try:
     234                if pData.sw != NULL:
     235                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
     236                else:
     237                    jac=(<object>pData.JAC)(t,y,p)
     238                   
     239                for i in range(Neq):
     240                    col_i = Jacobian.cols[i]
     241                    for j in range(Neq):
     242                        col_i[j] = jac[j,i]
     243
     244                return CVDLS_SUCCESS
     245            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     246                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     247            except:
     248                traceback.print_exc()
     249                return CVDLS_JACFUNC_UNRECVR
     250        else:
     251            try:
     252                if pData.sw != NULL:
     253                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
     254                else:
     255                    jac=(<object>pData.JAC)(t,y)
     256       
     257                for i in range(Neq):
     258                    col_i = Jacobian.cols[i]
     259                    for j in range(Neq):
     260                        col_i[j] = jac[j,i]
     261
     262                return CVDLS_SUCCESS
     263            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     264                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     265            except:
     266                traceback.print_exc()
     267                return CVDLS_JACFUNC_UNRECVR
     268ELSE:
     269    cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
     270                    void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
     271        """
     272        This method is used to connect the Assimulo.Problem.jac to the Sundials
     273        Jacobian function.
     274        """
     275        cdef ProblemData pData = <ProblemData>problem_data
     276        cdef realtype* col_i=DENSE_COL(Jacobian,0)
     277        cdef N.ndarray y = pData.work_y
     278        cdef int i,j
     279       
     280        nv2arr_inplace(yv, y)
     281
     282        if pData.dimSens>0: #Sensitivity activated
     283            p = realtype2arr(pData.p,pData.dimSens)
     284            try:
     285                if pData.sw != NULL:
     286                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
     287                else:
     288                    jac=(<object>pData.JAC)(t,y,p)
     289                   
     290                for i in range(Neq):
     291                    col_i = DENSE_COL(Jacobian, i)
     292                    for j in range(Neq):
     293                        col_i[j] = jac[j,i]
     294
     295                return CVDLS_SUCCESS
     296            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     297                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     298            except:
     299                traceback.print_exc()
     300                return CVDLS_JACFUNC_UNRECVR
     301        else:
     302            try:
     303                if pData.sw != NULL:
     304                    jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
     305                else:
     306                    jac=(<object>pData.JAC)(t,y)
     307       
     308                for i in range(Neq):
     309                    col_i = DENSE_COL(Jacobian, i)
     310                    for j in range(Neq):
     311                        col_i[j] = jac[j,i]
     312
     313                return CVDLS_SUCCESS
     314            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     315                return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     316            except:
     317                traceback.print_exc()
     318                return CVDLS_JACFUNC_UNRECVR
    213319       
    214320cdef int cv_jacv(N_Vector vv, N_Vector Jv, realtype t, N_Vector yv, N_Vector fyv,
     
    338444    """
    339445    cdef ProblemData pData = <ProblemData>problem_data
    340     #cdef ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
    341     #(<ndarray>pData.y).data =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
    342     #cdef N.ndarray y = nv2arr(yv)
    343446    cdef N.ndarray y = pData.work_y
    344447    cdef int i
     
    367470    cdef ProblemData pData = <ProblemData>problem_data
    368471    cdef N.ndarray[realtype, ndim=1, mode='c'] res #Used for return from the user function
    369     #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
    370     #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
    371472    cdef N.ndarray y = pData.work_y
    372473    cdef N.ndarray yd = pData.work_yd
    373     # cdef N.ndarray y = nv2arr(yv)
    374     # cdef N.ndarray yd = nv2arr(yvdot)
    375474    cdef realtype* resptr=(<N_VectorContent_Serial>residual.content).data
    376475    cdef int i
     
    415514            traceback.print_exc()
    416515            return IDA_RES_FAIL
    417            
    418 cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
     516
     517IF SUNDIALS_VERSION >= (3,0,0):
     518    cdef int ida_jac(realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, SUNMatrix Jac,
     519                 void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
     520        """
     521        This method is used to connect the Assimulo.Problem.jac to the Sundials
     522        Jacobian function.
     523        """
     524        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
     525        cdef ProblemData pData = <ProblemData>problem_data
     526        cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
     527        cdef realtype* col_i=Jacobian.cols[0]
     528        cdef N.ndarray y = pData.work_y
     529        cdef N.ndarray yd = pData.work_yd
     530        cdef int i,j, Neq = pData.dim
     531       
     532        nv2arr_inplace(yv, y)
     533        nv2arr_inplace(yvdot, yd)
     534       
     535        if pData.dimSens!=0: #SENSITIVITY
     536            p = realtype2arr(pData.p,pData.dimSens)
     537            try:
     538                if pData.sw != NULL:
     539                    jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
     540                else:
     541                    jac=(<object>pData.JAC)(c,t,y,yd,p=p)
     542               
     543                for i in range(Neq):
     544                    col_i = Jacobian.cols[i]
     545                    for j in range(Neq):
     546                        col_i[j] = jac[j,i]
     547                return IDADLS_SUCCESS
     548            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     549                return IDADLS_JACFUNC_RECVR #Recoverable Error
     550            except:
     551                traceback.print_exc()
     552                return IDADLS_JACFUNC_UNRECVR
     553        else:
     554            try:
     555                if pData.sw != NULL:
     556                    jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
     557                else:
     558                    jac=(<object>pData.JAC)(c,t,y,yd)
     559               
     560                for i in range(Neq):
     561                    col_i = Jacobian.cols[i]
     562                    for j in range(Neq):
     563                        col_i[j] = jac[j,i]
     564                return IDADLS_SUCCESS
     565            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     566                return IDADLS_JACFUNC_RECVR #Recoverable Error
     567            except:
     568                traceback.print_exc()
     569                return IDADLS_JACFUNC_UNRECVR
     570ELSE:
     571    cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
    419572                 void* problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
    420     """
    421     This method is used to connect the Assimulo.Problem.jac to the Sundials
    422     Jacobian function.
    423     """
    424     cdef ProblemData pData = <ProblemData>problem_data
    425     cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
    426     cdef realtype* col_i=DENSE_COL(Jacobian,0)
    427     #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
    428     #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
     573        """
     574        This method is used to connect the Assimulo.Problem.jac to the Sundials
     575        Jacobian function.
     576        """
     577        cdef ProblemData pData = <ProblemData>problem_data
     578        cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
     579        cdef realtype* col_i=DENSE_COL(Jacobian,0)
     580        cdef N.ndarray y = pData.work_y
     581        cdef N.ndarray yd = pData.work_yd
     582        cdef int i,j
     583       
     584        nv2arr_inplace(yv, y)
     585        nv2arr_inplace(yvdot, yd)
     586       
     587        if pData.dimSens!=0: #SENSITIVITY
     588            p = realtype2arr(pData.p,pData.dimSens)
     589            try:
     590                if pData.sw != NULL:
     591                    jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
     592                else:
     593                    jac=(<object>pData.JAC)(c,t,y,yd,p=p)
     594               
     595                for i in range(Neq):
     596                    col_i = DENSE_COL(Jacobian, i)
     597                    for j in range(Neq):
     598                        col_i[j] = jac[j,i]
     599                return IDADLS_SUCCESS
     600            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     601                return IDADLS_JACFUNC_RECVR #Recoverable Error
     602            except:
     603                traceback.print_exc()
     604                return IDADLS_JACFUNC_UNRECVR
     605        else:
     606            try:
     607                if pData.sw != NULL:
     608                    jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
     609                else:
     610                    jac=(<object>pData.JAC)(c,t,y,yd)
     611               
     612                for i in range(Neq):
     613                    col_i = DENSE_COL(Jacobian, i)
     614                    for j in range(Neq):
     615                        col_i[j] = jac[j,i]
     616                return IDADLS_SUCCESS
     617            except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
     618                return IDADLS_JACFUNC_RECVR #Recoverable Error
     619            except:
     620                traceback.print_exc()
     621                return IDADLS_JACFUNC_UNRECVR
     622           
     623
     624cdef int ida_root(realtype t, N_Vector yv, N_Vector yvdot, realtype *gout,  void* problem_data):
     625    """
     626    This method is used to connect the Assimulo.Problem.state_events to the Sundials
     627    root function.
     628    """
     629    cdef ProblemData pData = <ProblemData>problem_data
     630    cdef N.ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
    429631    cdef N.ndarray y = pData.work_y
    430632    cdef N.ndarray yd = pData.work_yd
    431     #cdef N.ndarray y = nv2arr(yv)
    432     #cdef N.ndarray yd = nv2arr(yvdot)
    433     cdef int i,j
    434    
    435     nv2arr_inplace(yv, y)
    436     nv2arr_inplace(yvdot, yd)
    437    
    438     if pData.dimSens!=0: #SENSITIVITY
    439         p = realtype2arr(pData.p,pData.dimSens)
    440         try:
    441             if pData.sw != NULL:
    442                 jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p)  # call to the python residual function
    443             else:
    444                 jac=(<object>pData.JAC)(c,t,y,yd,p=p)
    445            
    446             for i in range(Neq):
    447                 col_i = DENSE_COL(Jacobian, i)
    448                 for j in range(Neq):
    449                     col_i[j] = jac[j,i]
    450             return IDADLS_SUCCESS
    451         except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
    452             return IDADLS_JACFUNC_RECVR #Recoverable Error
    453         except:
    454             traceback.print_exc()
    455             return IDADLS_JACFUNC_UNRECVR
    456     else:
    457         try:
    458             if pData.sw != NULL:
    459                 jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw)  # call to the python residual function
    460             else:
    461                 jac=(<object>pData.JAC)(c,t,y,yd)
    462            
    463             for i in range(Neq):
    464                 col_i = DENSE_COL(Jacobian, i)
    465                 for j in range(Neq):
    466                     col_i[j] = jac[j,i]
    467             return IDADLS_SUCCESS
    468         except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
    469             return IDADLS_JACFUNC_RECVR #Recoverable Error
    470         except:
    471             traceback.print_exc()
    472             return IDADLS_JACFUNC_UNRECVR
    473        
    474 
    475 cdef int ida_root(realtype t, N_Vector yv, N_Vector yvdot, realtype *gout,  void* problem_data):
    476     """
    477     This method is used to connect the Assimulo.Problem.state_events to the Sundials
    478     root function.
    479     """
    480     cdef ProblemData pData = <ProblemData>problem_data
    481     cdef N.ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
    482     #(<ndarray>pData.y).data  =  <realtype*>((<N_VectorContent_Serial>yv.content).data)
    483     #(<ndarray>pData.yd).data =  <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
    484     cdef N.ndarray y = pData.work_y
    485     cdef N.ndarray yd = pData.work_yd
    486     #cdef N.ndarray y = nv2arr(yv)
    487     #cdef N.ndarray yd = nv2arr(yvdot)
    488633    cdef int i
    489634   
     
    554699            return SPGMR_PSOLVE_FAIL_UNREC
    555700   
    556 
    557 cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian,
    558                 void *problem_data, N_Vector tmp1, N_Vector tmp2):
    559     """
    560     This method is used to connect the assimulo.Problem.jac to the Sundials
    561     Jacobian function.
    562     """
    563     cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
    564     cdef realtype* col_i=DENSE_COL(Jacobian,0)
    565     cdef N.ndarray x = nv2arr(xv)
    566     cdef int i,j
    567    
    568     try:
    569         jac=(<object>pData.JAC)(x)
    570 
    571         for i in range(Neq):
    572             col_i = DENSE_COL(Jacobian, i)
    573             for j in range(Neq):
    574                 col_i[j] = jac[j,i]
    575 
    576         return KINDLS_SUCCESS
    577     except:
    578         return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
    579        
     701IF SUNDIALS_VERSION >= (3,0,0):
     702    cdef int kin_jac(N_Vector xv, N_Vector fval, SUNMatrix Jac,
     703                    void *problem_data, N_Vector tmp1, N_Vector tmp2):
     704        """
     705        This method is used to connect the assimulo.Problem.jac to the Sundials
     706        Jacobian function.
     707        """
     708        cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
     709        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
     710        cdef realtype* col_i=Jacobian.cols[0]
     711        cdef N.ndarray x = nv2arr(xv)
     712        cdef int i,j, Neq = pData.dim
     713       
     714        try:
     715            jac=(<object>pData.JAC)(x)
     716
     717            for i in range(Neq):
     718                col_i = Jacobian.cols[i]
     719                for j in range(Neq):
     720                    col_i[j] = jac[j,i]
     721
     722            return KINDLS_SUCCESS
     723        except:
     724            return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     725ELSE:
     726    cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian,
     727                    void *problem_data, N_Vector tmp1, N_Vector tmp2):
     728        """
     729        This method is used to connect the assimulo.Problem.jac to the Sundials
     730        Jacobian function.
     731        """
     732        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
     733        cdef realtype* col_i=DENSE_COL(Jacobian,0)
     734        cdef N.ndarray x = nv2arr(xv)
     735        cdef int i,j
     736       
     737        try:
     738            jac=(<object>pData.JAC)(x)
     739
     740            for i in range(Neq):
     741                col_i = DENSE_COL(Jacobian, i)
     742                for j in range(Neq):
     743                    col_i[j] = jac[j,i]
     744
     745            return KINDLS_SUCCESS
     746        except:
     747            return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
     748           
    580749cdef int kin_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
    581750            void *problem_data):
  • trunk/src/lib/sundials_constants.pxi

    r736 r845  
    9595DEF CV_BAD_IS             = -45
    9696
     97
     98DEF CSC_MAT = 0
     99DEF CSR_MAT = 1
    97100
    98101#==========
  • trunk/src/lib/sundials_includes.pxd

    r840 r845  
    3737    ctypedef double realtype
    3838    ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
    39 
    40 #==============================================
    41 # C headers
    42 #==============================================
    43 cdef extern from "string.h":
    44     void *memcpy(void *s1, void *s2, int n)
    45 cdef extern from "stdlib.h":
    46     void *malloc(int size)
    47     void free(void *ptr)
    48    
    49 #==============================================
    50 #External definitions from Sundials headers
    51 #==============================================
    5239
    5340cdef extern from "sundials/sundials_nvector.h":
     
    7764    void N_VDestroy_Serial(N_Vector v)
    7865    void N_VPrint_Serial(N_Vector v)
     66
     67IF SUNDIALS_VERSION >= (3,0,0):
     68    cdef extern from "sundials/sundials_types.h":
     69        IF SUNDIALS_VECTOR_SIZE == "64":
     70            ctypedef long int sunindextype
     71        ELSE:
     72            ctypedef int sunindextype
     73    cdef extern from "sundials/sundials_matrix.h":
     74        ctypedef _generic_SUNMatrix *SUNMatrix;
     75        void SUNMatDestroy(SUNMatrix A);
     76       
     77        cdef struct _generic_SUNMatrix_Ops:
     78            SUNMatrix_ID (*getid)(SUNMatrix);
     79            SUNMatrix    (*clone)(SUNMatrix);
     80            void         (*destroy)(SUNMatrix);
     81            int          (*zero)(SUNMatrix);
     82            int          (*copy)(SUNMatrix, SUNMatrix);
     83            int          (*scaleadd)(realtype, SUNMatrix, SUNMatrix);
     84            int          (*scaleaddi)(realtype, SUNMatrix);
     85            int          (*matvec)(SUNMatrix, N_Vector, N_Vector);
     86            int          (*space)(SUNMatrix, long int*, long int*);
     87
     88        cdef struct _generic_SUNMatrix:
     89            void *content;
     90            _generic_SUNMatrix_Ops *ops;
     91           
     92        cdef enum SUNMatrix_ID:
     93            SUNMATRIX_DENSE,
     94            SUNMATRIX_BAND,
     95            SUNMATRIX_SPARSE,
     96            SUNMATRIX_CUSTOM
     97   
     98    cdef extern from "sundials/sundials_linearsolver.h":
     99        ctypedef _generic_SUNLinearSolver *SUNLinearSolver;
     100        int SUNLinSolFree(SUNLinearSolver S);
     101       
     102        cdef struct _generic_SUNLinearSolver_Ops:
     103            SUNLinearSolver_Type (*gettype)(SUNLinearSolver);
     104            int                  (*setatimes)(SUNLinearSolver, void*, ATimesFn);
     105            int                  (*setpreconditioner)(SUNLinearSolver, void*,
     106                                                    PSetupFn, PSolveFn);
     107            int                  (*setscalingvectors)(SUNLinearSolver,
     108                                                    N_Vector, N_Vector);
     109            int                  (*initialize)(SUNLinearSolver);
     110            int                  (*setup)(SUNLinearSolver, SUNMatrix);
     111            int                  (*solve)(SUNLinearSolver, SUNMatrix, N_Vector,
     112                                        N_Vector, realtype);
     113            int                  (*numiters)(SUNLinearSolver);
     114            realtype             (*resnorm)(SUNLinearSolver);
     115            long int             (*lastflag)(SUNLinearSolver);
     116            int                  (*space)(SUNLinearSolver, long int*, long int*);
     117            N_Vector             (*resid)(SUNLinearSolver);
     118            int                  (*free)(SUNLinearSolver);
     119       
     120        cdef struct _generic_SUNLinearSolver:
     121            void *content;
     122            _generic_SUNLinearSolver_Ops *ops;
     123           
     124        cdef enum SUNLinearSolver_Type:
     125            SUNLINEARSOLVER_DIRECT,
     126            SUNLINEARSOLVER_ITERATIVE,
     127            SUNLINEARSOLVER_CUSTOM
     128   
     129    cdef extern from "sunmatrix/sunmatrix_dense.h":
     130        ctypedef _SUNMatrixContent_Dense *SUNMatrixContent_Dense;
     131        cdef struct _SUNMatrixContent_Dense:
     132            sunindextype M;
     133            sunindextype N;
     134            realtype *data;
     135            sunindextype ldata;
     136            realtype **cols;
     137        SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N);
     138    cdef extern from "sunmatrix/sunmatrix_sparse.h":
     139        ctypedef _SUNMatrixContent_Sparse *SUNMatrixContent_Sparse;
     140        cdef struct _SUNMatrixContent_Sparse:
     141            sunindextype M;
     142            sunindextype N;
     143            sunindextype NNZ;
     144            sunindextype NP;
     145            realtype *data;
     146            int sparsetype;
     147            sunindextype *indexvals;
     148            sunindextype *indexptrs;
     149            sunindextype **rowvals;
     150            sunindextype **colptrs;
     151            sunindextype **colvals;
     152            sunindextype **rowptrs;
     153        SUNMatrix SUNSparseMatrix(sunindextype M, sunindextype N, sunindextype NNZ, int sparsetype);
     154    cdef extern from "sunlinsol/sunlinsol_dense.h":
     155        SUNLinearSolver SUNDenseLinearSolver(N_Vector y, SUNMatrix A);
     156    cdef extern from "sunlinsol/sunlinsol_spgmr.h":
     157        SUNLinearSolver SUNSPGMR(N_Vector y, int pretype, int maxl);
     158       
     159ELSE:
     160    #Dummy defines
     161    ctypedef void *SUNLinearSolver
     162    ctypedef void *SUNMatrix
     163    ctypedef void *SUNMatrixContent_Dense
     164    ctypedef void *SUNMatrixContent_Sparse
     165    ctypedef int sunindextype
     166
    79167
    80168#Struct for handling the Jacobian data
     
    93181    ctypedef _DlsMat* DlsMat
    94182    cdef realtype* DENSE_COL(DlsMat A, int j)
    95 
     183   
    96184IF SUNDIALS_VERSION >= (2,6,3):
    97185    cdef extern from "sundials/sundials_sparse.h":
     
    129217        int *colptrs
    130218    ctypedef _SlsMat* SlsMat
     219   
     220#==============================================
     221# C headers
     222#==============================================
     223cdef extern from "string.h":
     224    void *memcpy(void *s1, void *s2, int n)
     225cdef extern from "stdlib.h":
     226    void *malloc(int size)
     227    void free(void *ptr)
     228   
     229#==============================================
     230#External definitions from Sundials headers
     231#==============================================
     232
     233IF SUNDIALS_WITH_SUPERLU:
     234    cdef inline int with_superlu(): return 1
     235ELSE:
     236    cdef inline int with_superlu(): return 0
     237
    131238
    132239cdef extern from "cvodes/cvodes.h":
     
    226333    int CVodeGetStgrSensNumNonlinSolvIters(void *cvode_mem, long int *nSTGR1niters)
    227334    int CVodeGetStgrSensNumNonlinSolvConvFails(void *cvode_mem, long int *nSTGR1ncfails)
    228    
    229    
    230 cdef extern from "cvodes/cvodes_dense.h":
    231     int CVDense(void *cvode_mem, long int n)
    232     ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
    233                    DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
    234     int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
    235 
    236 cdef extern from "cvodes/cvodes_spgmr.h":
    237     int CVSpgmr(void *cvode_mem, int pretype, int max1)
    238    
    239 IF SUNDIALS_VERSION >= (2,6,0):
    240     cdef extern from "cvodes/cvodes_sparse.h":
     335
     336cdef extern from "cvodes/cvodes_spils.h":
     337    ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
     338            N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
     339
     340IF SUNDIALS_VERSION >= (3,0,0):
     341    cdef extern from "cvodes/cvodes_direct.h":
     342        ctypedef int (*CVDlsDenseJacFn)(realtype t, N_Vector y, N_Vector fy,
     343                       SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
     344        int CVDlsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS, SUNMatrix A);
     345        int CVDlsSetJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
     346    cdef extern from "cvodes/cvodes_spils.h":
     347        int CVSpilsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS);
     348        ctypedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void *user_data);
     349        int CVSpilsSetJacTimes(void *cvode_mem, CVSpilsJacTimesSetupFn jtsetup, CVSpilsJacTimesVecFn jtimes);
     350   
     351   
     352    IF SUNDIALS_WITH_SUPERLU:
     353        cdef extern from "sunlinsol/sunlinsol_superlumt.h":
     354            SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
     355    ELSE:
     356        cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
     357   
     358    cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0   
     359    cdef inline tuple version(): return (3,0,0)
     360ELSE:
     361    cdef extern from "cvodes/cvodes_dense.h":
     362        int CVDense(void *cvode_mem, long int n)
     363        ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
     364                       DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
     365        int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
     366
     367    cdef extern from "cvodes/cvodes_spgmr.h":
     368        int CVSpgmr(void *cvode_mem, int pretype, int max1)
     369   
     370    cdef extern from "cvodes/cvodes_spils.h":
     371        int CVSpilsSetJacTimesVecFn(void *cvode_mem,  CVSpilsJacTimesVecFn jtv)
     372   
     373    IF SUNDIALS_VERSION >= (2,6,0):
     374        cdef extern from "cvodes/cvodes_sparse.h":
     375            ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
     376                                      SlsMat Jac, void *user_data, N_Vector tmp1,
     377                                        N_Vector tmp2, N_Vector tmp3)
     378            int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
     379            int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
     380        cdef inline tuple version(): return (2,6,0)
     381        IF SUNDIALS_WITH_SUPERLU:
     382            cdef extern from "cvodes/cvodes_sparse.h":
     383                int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
     384        ELSE:
     385            cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
     386    ELSE:
     387        cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
    241388        ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
    242389                                  SlsMat Jac, void *user_data, N_Vector tmp1,
    243390                                    N_Vector tmp2, N_Vector tmp3)
    244         int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
    245         int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
    246     #cdef inline char* version(): return "2.6.0"
    247     cdef inline tuple version(): return (2,6,0)
    248     IF SUNDIALS_WITH_SUPERLU:
    249         cdef extern from "cvodes/cvodes_sparse.h":
    250             int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
    251     ELSE:
    252         cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
    253 ELSE:
    254     cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
    255     ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
    256                               SlsMat Jac, void *user_data, N_Vector tmp1,
    257                                 N_Vector tmp2, N_Vector tmp3)
    258     cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
    259     cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
    260     #cdef inline char* version(): return "2.5.0"
    261     cdef inline tuple version(): return (2,5,0)
     391        cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
     392        cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
     393        cdef inline tuple version(): return (2,5,0)
    262394   
    263395cdef extern from "cvodes/cvodes_spils.h":
    264     ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
    265                                     N_Vector y, N_Vector fy,
    266                                     void *user_data, N_Vector tmp)
    267396    ctypedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy,
    268397                                  booleantype jok, booleantype *jcurPtr,
     
    274403                                  realtype gamma, realtype delta,
    275404                                  int lr, void *user_data, N_Vector tmp)
    276     int CVSpilsSetJacTimesVecFn(void *cvode_mem,  CVSpilsJacTimesVecFn jtv)
     405   
    277406    int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
    278407    int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
     
    280409    int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
    281410    int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
    282    
    283411
    284412cdef extern from "idas/idas.h":
     
    381509    #End Sensitivities
    382510    #=================
    383 
    384 cdef extern from "idas/idas_dense.h":
    385     int IDADense(void *ida_mem, long int n)
    386     ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
    387                    N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
    388                    N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
    389     int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
    390    
    391 cdef extern from "idas/idas_dense.h":
    392     int IDASpgmr(void *ida_mem, int max1)
    393     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)
    394     int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
    395    
     511   
     512cdef extern from "idas/idas_spils.h":
     513    ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr,
     514            N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
     515   
     516IF SUNDIALS_VERSION >= (3,0,0):
     517    cdef extern from "idas/idas_direct.h":
     518        ctypedef int (*IDADlsDenseJacFn)(realtype tt, realtype cj, N_Vector yy,
     519                       N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data,
     520                       N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
     521        int IDADlsSetJacFn(void *ida_mem, IDADlsDenseJacFn djac)
     522        int IDADlsSetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A);
     523   
     524    cdef extern from "idas/idas_spils.h":
     525        int IDASpilsSetLinearSolver(void *ida_mem, SUNLinearSolver LS);
     526        ctypedef int (*IDASpilsJacTimesSetupFn)(realtype tt, N_Vector yy,
     527                    N_Vector yp, N_Vector rr, realtype c_j, void *user_data);
     528        int IDASpilsSetJacTimes(void *ida_mem,
     529                IDASpilsJacTimesSetupFn jtsetup, IDASpilsJacTimesVecFn jtimes);
     530               
     531    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
     532ELSE:
     533    cdef extern from "idas/idas_dense.h":
     534        int IDADense(void *ida_mem, long int n)
     535        ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
     536                       N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
     537                       N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
     538        int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
     539   
     540        int IDASpgmr(void *ida_mem, int max1)
     541       
     542    cdef extern from "idas/idas_spils.h":
     543        int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
     544
    396545cdef extern from "idas/idas_spils.h":
    397546    int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
     
    456605    void KINFree(void **kinmem)
    457606
    458 # functions used for supplying jacobian, and receiving info from linear solver
     607
     608IF SUNDIALS_VERSION >= (3,0,0):
     609    cdef extern from "kinsol/kinsol_direct.h":
     610        ctypedef int (*KINDlsDenseJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
     611        int KINDlsSetLinearSolver(void *kinmem, SUNLinearSolver LS, SUNMatrix A);
     612        int KINDlsSetJacFn(void *kinmem, KINDlsDenseJacFn djac)
     613   
     614    cdef extern from "kinsol/kinsol_spils.h":
     615        int KINSpilsSetLinearSolver(void *kinsol_mem, SUNLinearSolver LS);
     616ELSE:
     617    # functions used for supplying jacobian, and receiving info from linear solver
     618    cdef extern from "kinsol/kinsol_direct.h":
     619        # user functions
     620        ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
     621       
     622        # function used to link user functions to KINSOL
     623        int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
     624   
     625    cdef extern from "kinsol/kinsol_dense.h":
     626        int KINDense(void *kinmem, int dim)
     627   
     628    cdef extern from "kinsol/kinsol_spgmr.h":
     629        int KINSpgmr(void *kinmem, int maxl)
     630
    459631cdef extern from "kinsol/kinsol_direct.h":
    460     # user functions
    461     ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
    462    
    463     # function used to link user functions to KINSOL
    464     int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
    465 
    466632    # optional output fcts for linear direct solver
    467633    int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
     
    470636    int KINDlsGetLastFlag(void *kinmem, int *flag)
    471637    char *KINDlsGetReturnFlagName(int flag)
    472 
    473 cdef extern from "kinsol/kinsol_dense.h":
    474     int KINDense(void *kinmem, int dim)
    475    
    476 cdef extern from "kinsol/kinsol_spgmr.h":
    477     int KINSpgmr(void *kinmem, int maxl)
    478638
    479639cdef extern from "kinsol/kinsol_spils.h":
  • trunk/src/solvers/__init__.py

    r838 r845  
    2626    from .euler import ImplicitEuler
    2727except ImportError as ie:
    28     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     28    sys.stderr.write("Could not find " + str(ie) + "\n")
    2929try:
    3030    from .radau5 import Radau5ODE
     
    3333    from .radau5 import _Radau5DAE
    3434except ImportError as ie:
    35     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     35    sys.stderr.write("Could not find " + str(ie) + "\n")
    3636try:
    3737    from .sundials import IDA
    3838    from .sundials import CVode
    3939except ImportError as ie:
    40     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     40    sys.stderr.write("Could not find " + str(ie) + "\n")
    4141try:
    4242    from .kinsol import KINSOL
    4343except ImportError as ie:
    44     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     44    sys.stderr.write("Could not find " + str(ie) + "\n")
    4545try:
    4646    from .runge_kutta import RungeKutta34
     
    4848    from .runge_kutta import Dopri5
    4949except ImportError as ie:
    50     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     50    sys.stderr.write("Could not find " + str(ie) + "\n")
    5151try:
    5252    from .rosenbrock import RodasODE
    5353except ImportError as ie:
    54     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     54    sys.stderr.write("Could not find " + str(ie) + "\n")
    5555try:
    5656    from .odassl import ODASSL
    5757except ImportError as ie:
    58     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     58    sys.stderr.write("Could not find " + str(ie) + "\n")
    5959try:
    6060    from .odepack import LSODAR
    6161except ImportError as ie:
    62     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     62    sys.stderr.write("Could not find " + str(ie) + "\n")
    6363try:
    6464    from .radar5 import Radar5ODE
    6565except ImportError as ie:
    66     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     66    sys.stderr.write("Could not find " + str(ie) + "\n")
    6767try:
    6868    from .dasp3 import DASP3ODE
    6969except ImportError as ie:
    70     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     70    sys.stderr.write("Could not find " + str(ie) + "\n")
    7171try:
    7272    from .glimda import GLIMDA
    7373except ImportError as ie:
    74     sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
     74    sys.stderr.write("Could not find " + str(ie) + "\n")
  • trunk/src/solvers/kinsol.pyx

    r742 r845  
    3131
    3232#Various C includes transfered to namespace
    33 from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
    34 from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
     33from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
     34from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
    3535from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
    3636from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
     
    5858    cdef object pt_fcn, pt_jac, pt_jacv, pt_prec_setup, pt_prec_solve
    5959    cdef object _added_linear_solver
     60    cdef SUNDIALS.SUNMatrix sun_matrix
     61    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
    6062   
    6163    def __init__(self, problem):
     
    8789        self.options["max_beta_fails"] = 10
    8890        self.options["max_krylov"] = 0
     91        self.options["precond"] = PREC_NONE
    8992       
    9093        #Statistics
     
    110113            #Free Memory
    111114            SUNDIALS.KINFree(&self.kinsol_mem)
     115           
     116        IF SUNDIALS_VERSION >= (3,0,0):
     117            if self.sun_matrix != NULL:
     118                SUNDIALS.SUNMatDestroy(self.sun_matrix)
     119               
     120            if self.sun_linearsolver != NULL:
     121                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
    112122       
    113123    def update_variable_scaling(self, value="Automatic"):
     
    171181    cdef initialize_kinsol(self):
    172182        cdef int flag #Used for return
    173        
     183
    174184        self.y_temp  = arr2nv(self.y)
    175185        self.y_scale = arr2nv([1.0]*self.problem_info["dim"])
     
    209219    cpdef add_linear_solver(self):
    210220        if self.options["linear_solver"] == "DENSE":
    211             flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
     221            IF SUNDIALS_VERSION >= (3,0,0):
     222                #Create a dense Sundials matrix
     223                self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
     224                #Create a dense Sundials linear solver
     225                self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.y_temp, self.sun_matrix)
     226                #Attach it to Kinsol
     227                flag = SUNDIALS.KINDlsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver, self.sun_matrix)
     228            ELSE:
     229                flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
    212230            if flag < 0:
    213231                raise KINSOLError(flag)
    214232           
    215233            if self.problem_info["jac_fcn"]:
    216                 flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
     234                IF SUNDIALS_VERSION >= (3,0,0):
     235                    flag = SUNDIALS.KINDlsSetJacFn(self.kinsol_mem, kin_jac);
     236                ELSE:
     237                    flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
    217238                if flag < 0:
    218239                    raise KINSOLError(flag)
    219240        elif self.options["linear_solver"] == "SPGMR":
    220             flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
     241            IF SUNDIALS_VERSION >= (3,0,0):
     242                #Create the linear solver
     243                self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.y_temp, self.options["precond"], self.options["max_krylov"])
     244                #Attach it to Kinsol
     245                flag = SUNDIALS.KINSpilsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver)
     246            ELSE:
     247                #Specify the use of KINSpgmr linear solver.
     248                flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
    221249            if flag < 0:
    222250                raise KINSOLError(flag)
  • trunk/src/solvers/sundials.pyx

    r844 r845  
    3535
    3636#Various C includes transfered to namespace
    37 from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
    38 from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
     37from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
     38from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
    3939from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
    4040from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
     
    7171    #cdef N.ndarray _event_info
    7272    cdef public N.ndarray g_old
     73    cdef SUNDIALS.SUNMatrix sun_matrix
     74    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
    7375   
    7476    def __init__(self, problem):
     
    101103        self.options["pbar"] = [1]*self.problem_info["dimSens"]
    102104        self.options["external_event_detection"] = False #Sundials rootfinding is used for event location as default
     105        self.options["precond"] = PREC_NONE
     106
    103107
    104108        #Solver support
     
    174178            #Free Memory
    175179            SUNDIALS.IDAFree(&self.ida_mem)
     180       
     181        IF SUNDIALS_VERSION >= (3,0,0):
     182            if self.sun_matrix != NULL:
     183                SUNDIALS.SUNMatDestroy(self.sun_matrix)
     184               
     185            if self.sun_linearsolver != NULL:
     186                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
    176187   
    177188    cpdef state_event_info(self):
     
    256267                raise IDAError(flag, self.t)
    257268               
    258             #Specify the use of the internal dense linear algebra functions.
    259             flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
    260             if flag < 0:
    261                 raise IDAError(flag, self.t)
    262                
    263                     #Choose a linear solver if and only if NEWTON is choosen
     269            #Choose a linear solver if and only if NEWTON is choosen
    264270            if self.options["linear_solver"] == 'DENSE':
    265                 #Specify the use of the internal dense linear algebra functions.
    266                 flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
     271                IF SUNDIALS_VERSION >= (3,0,0):
     272                    #Create a dense Sundials matrix
     273                    self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
     274                    #Create a dense Sundials linear solver
     275                    self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix)
     276                    #Attach it to IDA
     277                    flag = SUNDIALS.IDADlsSetLinearSolver(self.ida_mem, self.sun_linearsolver, self.sun_matrix);
     278                ELSE:
     279                    #Specify the use of the internal dense linear algebra functions.
     280                    flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim)
    267281                if flag < 0:
    268282                    raise IDAError(flag, self.t)
    269283                       
    270284            elif self.options["linear_solver"] == 'SPGMR':
    271                 #Specify the use of SPGMR linear solver.
    272                 flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations
     285                IF SUNDIALS_VERSION >= (3,0,0):
     286                    #Create the linear solver
     287                    self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], 0)
     288                    #Attach it to IDAS
     289                    flag = SUNDIALS.IDASpilsSetLinearSolver(self.ida_mem, self.sun_linearsolver)
     290                ELSE:
     291                    #Specify the use of SPGMR linear solver.
     292                    flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations
    273293                if flag < 0:
    274294                    raise IDAError(flag, self.t)
     
    311331            #Specify the jacobian to the solver
    312332            if self.pData.JAC != NULL and self.options["usejac"]:
    313                
    314                 flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac)
     333                IF SUNDIALS_VERSION >= (3,0,0):
     334                    flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, ida_jac)
     335                ELSE:
     336                    flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac)
    315337                if flag < 0:
    316338                    raise IDAError(flag,self.t)
    317339            else:
    318                 flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL)
     340                IF SUNDIALS_VERSION >= (3,0,0):
     341                    flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, NULL)
     342                ELSE:
     343                    flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL)
    319344                if flag < 0:
    320345                    raise IDAError(flag,self.t)
     
    323348            #Specify the jacobian times vector function
    324349            if self.pData.JACV != NULL and self.options["usejac"]:
    325                 flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv);
     350                IF SUNDIALS_VERSION >= (3,0,0):
     351                    flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, SUNDIALS.ida_spils_jtsetup_dummy, ida_jacv);
     352                ELSE:
     353                    flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv);
    326354                if flag < 0:
    327355                    raise IDAError(flag, self.t)
    328356            else:
    329                 flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL);
     357                IF SUNDIALS_VERSION >= (3,0,0):
     358                    flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, NULL, NULL);
     359                ELSE:
     360                    flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL);
    330361                if flag < 0:
    331362                    raise IDAError(flag, self.t)
     
    14161447    #cdef N.ndarray _event_info
    14171448    cdef public N.ndarray g_old
     1449    cdef SUNDIALS.SUNMatrix sun_matrix
     1450    cdef SUNDIALS.SUNLinearSolver sun_linearsolver
    14181451   
    14191452    def __init__(self, problem):
     
    14801513            #Free Memory
    14811514            SUNDIALS.CVodeFree(&self.cvode_mem)
     1515       
     1516        IF SUNDIALS_VERSION >= (3,0,0):
     1517            if self.sun_matrix != NULL:
     1518                SUNDIALS.SUNMatDestroy(self.sun_matrix)
     1519               
     1520            if self.sun_linearsolver != NULL:
     1521                SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
    14821522   
    14831523    cpdef get_local_errors(self):
     
    20202060        """
    20212061        cdef flag
    2022        
     2062
    20232063        #Choose a linear solver if and only if NEWTON is choosen
    20242064        if self.options["linear_solver"] == 'DENSE' and self.options["iter"] == "Newton":
    2025             #Specify the use of the internal dense linear algebra functions.
    2026             flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim)
    2027             if flag < 0:
    2028                 raise CVodeError(flag)
     2065            IF SUNDIALS_VERSION >= (3,0,0):
     2066                #Create a dense Sundials matrix
     2067                self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
     2068                #Create a dense Sundials linear solver
     2069                self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix)
     2070                #Attach it to CVode
     2071                flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix);
     2072                if flag < 0:
     2073                    raise CVodeError(flag)
     2074            ELSE:
     2075                #Specify the use of the internal dense linear algebra functions.
     2076                flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim)
     2077                if flag < 0:
     2078                    raise CVodeError(flag)
    20292079               
    20302080            #Specify the jacobian to the solver
    20312081            if self.pData.JAC != NULL and self.options["usejac"]:
    2032                 flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac)
     2082                IF SUNDIALS_VERSION >= (3,0,0):
     2083                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac);
     2084                ELSE:
     2085                    flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac)
    20332086                if flag < 0:
    20342087                    raise CVodeError(flag)
    20352088            else:
    2036                 flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL)
     2089                IF SUNDIALS_VERSION >= (3,0,0):
     2090                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, NULL);
     2091                ELSE:
     2092                    flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL)
    20372093                if flag < 0:
    20382094                    raise CVodeError(flag)
    20392095                   
    20402096        elif self.options["linear_solver"] == 'SPGMR' and self.options["iter"] == "Newton":
    2041             #Specify the use of CVSPGMR linear solver.
    2042             flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"])
     2097            IF SUNDIALS_VERSION >= (3,0,0):
     2098                #Create the linear solver
     2099                self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], self.options["maxkrylov"])
     2100                #Attach it to CVode
     2101                flag = SUNDIALS.CVSpilsSetLinearSolver(self.cvode_mem, self.sun_linearsolver)
     2102            ELSE:
     2103                #Specify the use of CVSPGMR linear solver.
     2104                flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"])
    20432105            if flag < 0:
    2044                 raise CVodeError(flag)
     2106                raise CVodeError(flag) 
    20452107               
    20462108            if self.pData.PREC_SOLVE != NULL:
     
    20562118            #Specify the jacobian times vector function
    20572119            if self.pData.JACV != NULL and self.options["usejac"]:
    2058                 flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv)
    2059                 if flag < 0:
     2120                IF SUNDIALS_VERSION >= (3,0,0):
     2121                    flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, SUNDIALS.cv_spils_jtsetup_dummy, cv_jacv);
     2122                ELSE:
     2123                    flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv)
     2124                if flag < 0:
    20602125                    raise CVodeError(flag)
    20612126            else:
    2062                 flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL)
     2127                IF SUNDIALS_VERSION >= (3,0,0):
     2128                    flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, NULL, NULL);
     2129                ELSE:
     2130                    flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL)
    20632131                if flag < 0:
    20642132                    raise CVodeError(flag)
     
    20672135            if SUNDIALS.version() < (2,6,0):
    20682136                raise AssimuloException("Not supported with this SUNDIALS version.")
    2069            
     2137            if SUNDIALS.with_superlu() == 0:
     2138                raise AssimuloException("No support for SuperLU was detected, please verify that SuperLU and SUNDIALS has been installed correctly.")
     2139               
    20702140            #Specify the use of CVSPGMR linear solver.
    20712141            if self.problem_info["jac_fcn_nnz"] == -1:
    20722142                raise AssimuloException("Need to specify the number of non zero elements in the Jacobian via the option 'jac_nnz'")
    2073             flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"])
     2143               
     2144            IF SUNDIALS_VERSION >= (3,0,0):
     2145                self.sun_matrix = SUNDIALS.SUNSparseMatrix(self.pData.dim, self.pData.dim, self.problem_info["jac_fcn_nnz"], CSC_MAT)
     2146                self.sun_linearsolver = SUNDIALS.SUNSuperLUMT(self.yTemp, self.sun_matrix, self.options["num_threads"])
     2147                flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix)
     2148            ELSE:
     2149                flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"])
    20742150            if flag < 0:
    20752151                    raise CVodeError(flag)
     
    20772153            #Specify the jacobian to the solver
    20782154            if self.pData.JAC != NULL and self.options["usejac"]:
    2079                 flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse)
     2155                IF SUNDIALS_VERSION >= (3,0,0):
     2156                    flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac_sparse)
     2157                ELSE:
     2158                    flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse)
    20802159                if flag < 0:
    20812160                    raise CVodeError(flag)
     
    29002979            self.statistics["njacvecs"]  += njvevals
    29012980        elif self.options["linear_solver"] == "SPARSE":
    2902             flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals)
     2981            IF SUNDIALS_VERSION >= (3,0,0):
     2982                flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals)
     2983            ELSE:
     2984                flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals)
    29032985            self.statistics["njacs"]   += njevals
    29042986        else:
  • trunk/tests/test_examples.py

    r805 r845  
    2323class Test_Examples:
    2424   
     25   
     26    @testattr(stddist = True)
     27    def test_cvode_with_jac_sparse(self):
     28        try:
     29            cvode_with_jac_sparse.run_example(with_plots=False)
     30        except AssimuloException:
     31            pass #Handle the case when SuperLU is not installed
     32   
    2533    @testattr(stddist = True)
    2634    def test_ida_with_user_defined_handle_result(self):
Note: See TracChangeset for help on using the changeset viewer.