# Has anyone used DVODE and is it good? **URL:** \ **Category:** Uncategorized\ **Created:** [November 16, 2021, 9:08am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281 "2021-11-16T09:08:06Z")\ **Posts on this page:** 20\ **Page:** 2
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [January 9, 2022, 2:50am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/21 "2022-01-09T02:50:12Z")
Note that the argument K of DVINDY, which is the order of the derivative, can be zero. In fact, subroutine DVODE makes calls with K = 0 to DVINDY\_CORE. Such calls interpolate the function values. The user may make calls with other values of K if derivative values are wanted. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 9, 2022, 4:09am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/22 "2022-01-09T04:09:57Z")
Thank you very much @ivanpribec @mecej4 ! I really appreciate your insights! I have tried ITASK = 2 and other values for DVODE, it seems no matter what number ITASK is, the current DVODE cannot do the interpolation I mentioned. I guess like you suggested, DVINDY may be more for some dense output, and calculate some derivatives. For now, I will just use DVODE and do integrals from 0 to 2, 2 to 4, 4 to 6, … 22 to 24, instead of do 2 to 24 then do interpolation at 2,4,6,8,…,22,24. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [January 10, 2022, 2:43am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/23 "2022-01-10T02:43:03Z")
Here’s an example of using `dvindy` for dense output. > **Example** > > ```auto > program droplet > > implicit none > > integer, parameter :: dp = kind(1.0d0) > real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp) > > integer, parameter :: neq = 3 > > integer :: itol, itask, istate, iopt, mf, iflag, n, ngrid > real(dp) :: rtol, s, send, ds > real(dp), dimension(neq) :: atol, y > > real(dp), allocatable :: ygrid(:,:), sgrid(:) > > real(dp), allocatable :: rwork(:) > integer, allocatable :: iwork(:) > integer :: lrw, liw > > real(dp) :: rpar(2) > integer :: idum(1) > > external :: dvode, dvindy > > rpar(1) = 3.3888_dp > rpar(2) = 0.8414_dp > > s = 0.0_dp > send = 2.0_dp > > ngrid = 41 > allocate(sgrid(ngrid),ygrid(neq,ngrid)) > sgrid(1) = 0.0_dp > ygrid(:,1) = [0.0_dp, 0.0_dp, 0.0_dp] > > ds = (send - s)/real(ngrid - 1,dp) > do n = 2, ngrid > sgrid(n) = (n - 1)*ds > end do > n = 2 > > ds = 0.00001_dp > s = ds > > y(1) = ds > y(2) = cos(y(1))*ds > y(3) = sin(y(1))*ds > > itol = 2 > rtol = 1.0e-7_dp > atol(1) = 1.0e-6_dp > atol(2) = 1.0e-6_dp > atol(3) = 1.0e-6_dp > > mf = 22 > lrw = len_rwork(neq,mf) > liw = len_iwork(neq,mf) > allocate(rwork(lrw), iwork(liw)) > > istate = 1 > itask = 2 > iopt = 0 > > integrate: do > > call dvode(young_laplace, neq, y, s, send, & > itol, rtol, atol, itask, istate, iopt, & > rwork, lrw, iwork, liw, & > young_laplace_jacobian, & > mf, rpar, idum) > > ! Automatic stepping > write(*,'(4(ES12.4))') s, y > > if (istate < 0) then > write(*,*) ' Failed at s = ',s,' with istate = ',istate > error stop > end if > > ! Interpolation > interpolate: do while(s >= sgrid(n)) > call dvindy (sgrid(n), 0, rwork(21), neq, ygrid(:,n), iflag) > n = n + 1 > if (n > ngrid) exit interpolate > end do interpolate > > if (s >= send) then > exit integrate > end if > > end do integrate > > write(*,*) > write(*,*) > > do n = 1, ngrid > write(*,'(4(ES12.4))') sgrid(n), ygrid(:,n) > end do > > contains > > subroutine young_laplace(neq, t, y, ydot, rpar, ipar) > integer, intent(in) :: neq > real(dp), intent(in) :: t, y(neq) > real(dp), intent(out) :: ydot(neq) > real(dp), intent(in) :: rpar(*) > integer, intent(in) :: ipar(*) > > associate(pL => rpar(1), drho => rpar(2)) > ydot(1) = pL - drho*y(3) - sin(y(1))/y(2) > ydot(2) = cos(y(1)) > ydot(3) = sin(y(1)) > end associate > > end subroutine > > subroutine young_laplace_jacobian(neq, t, y, ml, mu, pd, nrpd, rpar, ipar) > integer, intent(in) :: neq, nrpd, ml, mu > real(dp), intent(out) :: pd(nrpd,neq) > real(dp), intent(in) :: t, y(neq), rpar(*) > integer, intent(in) :: ipar(*) > > ! dummy routine > associate(pL => rpar(1), drho => rpar(2)) > pd(1,1) = -cos(y(1))/y(2) > pd(2,1) = -sin(y(1)) > pd(3,1) = cos(y(1)) > > pd(1,2) = sin(y(1))/(y(2))**2 > pd(2,2) = 0.0_dp > pd(3,2) = 0.0_dp > > pd(1,3) = -drho > pd(2,3) = 0.0_dp > pd(3,3) = 0.0_dp > end associate > > end subroutine > > function len_rwork(neq, mf, mlu) result(lrw) > integer, intent(in) :: neq, mf > integer, intent(in), optional :: mlu(2) > ! Upper and lower bandwidth in case of banded Jacobian > > integer :: lrw > > ! RWORK = Real work array of length at least: > ! 20 + 16*NEQ for MF = 10, > ! 22 + 9*NEQ + 2*NEQ**2 for MF = 21 or 22, > ! 22 + 11*NEQ + (3*ML + 2*MU)*NEQ for MF = 24 or 25. > > select case(mf) > case(10) > lrw = 20 + 17*neq > case(21,22) > lrw = 22 + 9*neq + 2*neq**2 > case(24,25) > if (present(mlu)) then > lrw = 22 + 11*neq + (3*mlu(1) + 2*mlu(2))*neq > else > write(*,*) "len_rwork: Missing array mlu containing bandwidths [ml, mu]" > error stop > end if > case default > ! Wrong value for mf > lrw = -1 > end select > > end function > > function len_iwork(neq, mf) result(liw) > integer, intent(in) :: neq, mf > > integer :: liw > > ! IWORK = Integer work array of length at least: > ! 30 for MF = 10, > ! 30 + NEQ for MF = 21, 22, 24, or 25. > ! If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower > ! and upper half-bandwidths ML,MU. > > select case (mf) > case (10) > liw = 30 > case (21, 22, 24, 25) > liw = 30 + neq > case default > ! Wrong value for mf > liw = -1 > end select > > end function > > end program > > ``` The system solved are the Young-Laplace equations for a pendant drop in a cylindrical coordinate system: \frac{d\phi}{ds} = \tilde{p}\_L - \Delta\tilde{\rho}\,z - \frac{\sin{\phi}}{r} \\ \frac{dr}{ds} = \cos{\phi}\\ \frac{dz}{ds} = \sin{\phi} An equivalent system of equations (but with a slightly different set of parameters) was first solved in 1883 by Francis Bashforth and John Couch Adams in [An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid - GDZ](http://resolver.sub.uni-goettingen.de/purl?PPN518083195). Notably, for this occasion Adams derived the family of explicit multi-step methods, now called Adams-Bashforth methods. In VODE, when using `mf = 10` the solver uses an implicit method from Adams. Here is plot of the interpolated (r,z) pairs super-imposed on the automatic step control: ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/a/a607bbe0e49608c50d37720ec6784eac8b955e68.png) The resulting profile traces half of a pendant drop. The interpolated points are spaced equidistantly in terms of the arc length s. By fitting the theoretical drop shape to an experimental shape profile, the surface tension of a liquid can be determined. The method is called pendant drop [tensiometry](https://en.wikipedia.org/wiki/Tensiometer_(surface_tension)). ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 10, 2022, 6:50am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/24 "2022-01-10T06:50:53Z")
That is awesome @ivanpribec . Thank you my man! Cool, I see. So you set ITASK = 2, so all the automated s and y are saved. Then you interpolate y at corresponding s. Brilliant. I think in this way as you did, we actually can achieve what a normal interpolation ODE would do for DVODE. Save all the time (in your case s) and the y array along the way between the initial and final time. Then do interpolation at corresponding s between the initial and final time! Thanks a again! Very impressive! PS. I used your droplet.f90 combined with the dvodeoriginal.f, I use intel fortran. In the debug mode, it shows an error ```auto droplet.f90(67): error #8284: If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic. [RTOL] ``` ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/1/1a1aa1999ad83df1eb21200f5b56f6cb7854a414.png) But in release mode it compile just fine and run perfectly. But I checked, it does looks like your rtol is a scalar, do you know what is wrong? Perhaps it is a Intel compiler bug or somthing? ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [January 10, 2022, 12:18pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/25 "2022-01-10T12:18:34Z")
> [@CRquantum](#): > > PS. I used your droplet.f90 combined with the dvodeoriginal.f, I use intel fortran. > In the debug mode, it shows an error I was relying on an external interface. The issue originates here: ```fortran SUBROUTINE DVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF, 2 RPAR, IPAR) EXTERNAL F, JAC DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK, RPAR INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, 1 MF, IPAR DIMENSION Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW), 1 RPAR(*), IPAR(*) ``` since `RTOL(*)` is an asummed-size array and not a scalar like in my calling code. Easiest way to get rid of the error is to change `rtol` from a scalar to a one-element array in the calling code: ```auto real(dp) :: rtol(1) ! Array rtol ``` ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 11, 2022, 1:07am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/26 "2022-01-11T01:07:03Z")
Thank you @ivanpribec ! That is great! My bad, did not look carefully enough that RTOL is assumed-size, LOL. It seems for the Jacobian matrix, it is by default 0 in DVODE, so there is no need to set 0 elements to PD, pd(2,2), pd(2,3), pd(3,2), pd(3,3) can be commented. ```auto associate(pL => rpar(1), drho => rpar(2)) pd(1,1) = -cos(y(1))/y(2) pd(2,1) = -sin(y(1)) pd(3,1) = cos(y(1)) pd(1,2) = sin(y(1))/(y(2))**2 !pd(2,2) = 0.0_dp !pd(3,2) = 0.0_dp pd(1,3) = -drho !pd(2,3) = 0.0_dp !pd(3,3) = 0.0_dp end associate ``` By the way, for stiff problem, providing a Jacobian matrix .vs. without Jacobian matrix, in your experience, will that cause a huge performance difference? I found that with Jacobian matrix, there is basically no particular performance gain, but it does improve the DVODE’s stability and ability of reaching convergence. Thanks! PS. Besides DVINDY, every time I can learn something new from your code 🙂 I learned function/subroutine can contain function/subroutine inside them before. I apply it in combining DVODE and FLINT ODE solver, such that the function FEX and the Jacobian JEX can be defined and used in FLINT ODE class. With interpolation the DVODE in FLINT can be, ```auto subroutine dvode_interpolate(me, T0, TOUT, Y0, Y, K, T_array, Y_array, Params) class(PM_ODE_Sys), intent(inout) :: me !< Differential Equation object real(kind = r8), intent(in), dimension(:), optional :: Params real(kind =r8) :: T0, T, TOUT, RTOL real(kind=r8) :: Y0(:), Y(:) integer :: NEQ, ITASK, ISTATE, IFLAG real(kind=r8) :: T_array(:), Y_array(:,:) integer :: K integer(kind=i8) :: i, ngrid NEQ = me%n Y = Y0 T = T0 ITASK = 2 ISTATE = 1 ngrid = size(T_array,kind=i8) i = 1 integrate : do while (T <= TOUT) CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX) IF (ISTATE<0) THEN write(6,*) 'DVODE Failed at ',T,' with istate = ',istate error stop ELSE ! do interpolation interpolate: do while(T >= T_array(i)) call DVINDY (T_array(i), K, Y_array(:,i),IFLAG) i = i + 1 if (i > ngrid) exit interpolate end do interpolate END IF enddo integrate return contains SUBROUTINE FEX(NEQ,T,Y,YDOT) INTEGER, INTENT (IN) :: NEQ DOUBLE PRECISION, INTENT (IN) :: T DOUBLE PRECISION, INTENT (IN) :: Y(NEQ) DOUBLE PRECISION, INTENT (OUT) :: YDOT(NEQ) YDOT = me%F(T,Y,Params) RETURN END SUBROUTINE FEX SUBROUTINE JEX(NEQ,T,Y,ML,MU,PD,NRPD) INTEGER, INTENT (IN) :: NEQ, ML, MU, NRPD DOUBLE PRECISION, INTENT (IN) :: T DOUBLE PRECISION, INTENT (IN) :: Y(NEQ) DOUBLE PRECISION, INTENT (OUT) :: PD(NRPD,NEQ) PD = me%J(NEQ,T,Y,ML,MU,NRPD) RETURN END SUBROUTINE JEX end subroutine dvode_interpolate ``` where me%F is the differential equations and me% is the Jacobian matrix defined in FLINT ODE class me. This time I learned associate, ```auto associate(pL => rpar(1), drho => rpar(2)) ``` Cool! ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 17, 2022, 6:04am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/27 "2022-01-17T06:04:57Z")
> [@ivanpribec](#): > > Easiest way to get rid of the error is to change `rtol` from a scalar to a one-element array in the calling code: > > ```auto > real(dp) :: rtol(1) ! Array rtol > > ``` Thank you @ivanpribec ! I just learned from Dr. Fortran @sblionel the […] trick, > [@How to fix the errors in compiling ODEPACK in debug mode?](https://fortran-lang.discourse.group/t/how-to-fix-the-errors-in-compiling-odepack-in-debug-mode/2597/7): > > As I wrote earlier, Fortran does not allow you to pass a scalar to an array. And, again, the reason it works in release mode is that a release configuration doesn’t enable the cross-source-file generated interface check. If the value is an input, and not written to, you can pass [int(RWORK(LWM))] , making it an array constructor. Any changes made in the called routine will be discarded on return. The below code also work, just add [] around atol and rtol to make array interface when call dvode (I see this trick in FLINT ode solver by @prince_mahajan too), ```auto call dvode(young_laplace, neq, y, s, send, & itol, [rtol], [atol], itask, istate, iopt, & rwork, lrw, iwork, liw, & young_laplace_jacobian, & mf, rpar, idum) ``` By the way, have you by any chance used DVODPK? You know like LSODPK, etc, those end with PK or K, means using some thing called Preconditioned Krylov methods? LLNL’s ODEPACK webpage mentioned DVODPK as a variant for DVODE, [https://computing.llnl.gov/projects/odepack](https://computing.llnl.gov/projects/odepack) It looks like CVODE in SUNDIALS is basically VODE + VODPK written in C and added some parallelization stuff. Do you know if those Preconditioned Krylov methods efficient? Thanks a lot! ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [January 17, 2022, 10:13am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/28 "2022-01-17T10:13:58Z")
**PLEASE NOTE:** VODE has been widely used for decades, so there are several diverging versions and descendants that a new user may try, without being aware of the other versions. Here are some comments, mostly from my recollections, and to be verified before being acted upon. In particular, if someone plans to use this solver for large problems (high order coupled differential equations in many variables), they should spend some time reading reports and papers before adopting a particular version… The version from Netlib uses the venerable Linpack solvers DGEFA/DGESV (for dense matrix representation) and DGBFA/DGBSV (for band matrix representation). Consider replacing these by the following Lapack routines: DGETRF/DGETRS and DGBTRF/DGBTRS, especially if you are using a library such as Intel’s MKL. The F90 version from Byrne and Thompson, of which @CRQuantum has a Gitlab backup, allows the use of a modified version of MA28 from the HSL library. MA28 is itself quite old, and there is a successor, MA48, available (licensed, but free for academic use). Look at Sundials from LLNL, which is a modern successor to VODE. It is written in C, but has hooks for calling from Fortran, and allows for parallelization. A large number of other modern sparse solvers are available, including Pardiso, SuperLU and SuiteSparse, which are direct solvers. Other iterative solvers based on Krylov-Arnoldi are available, as well. Thus, a prospective user of VODE should devote some time to planning the work before committing to one combination, trying it out on a small set of small problems, and then scaling it up for testing on large systems. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 17, 2022, 7:52pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/29 "2022-01-17T19:52:05Z")
Thank you very much @mecej4 ! You comments are very helpful! > [@mecej4](#): > > The version from Netlib uses the venerable Linpack solvers DGEFA/DGESV (for dense matrix representation) and DGBFA/DGBSV (for band matrix representation). Consider replacing these by the following Lapack routines: DGETRF/DGETRS and DGBTRF/DGBTRS, especially if you are using a library such as Intel’s MKL. > > The F90 version from Byrne and Thompson, of which @CRQuantum has a Gitlab backup, allows the use of a modified version of MA28 from the HSL library. MA28 is itself quite old, and there is a successor, MA48, available (licensed, but free for academic use). RIght. Actually the F90 version from Byrne and Thompson natively support LAPACK with DGETRF and DGBTRF, as well as MA48. This may be the best DVODE version. There are some changing point in Byrne and Thompson code, just make some changes according to the instructions in the code, e.g., ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/1/185d51d389b4357d07a036d361e852c7f0d9f2fe.png) ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/7/70001073b2e53fbab323d059c90ff0615921e870.png) I have tested, the lapack or intel msl lapack works fine 🙂 The MA48, seems need to be download from HSL, [https://www.hsl.rl.ac.uk/catalogue/hsl\_ma48.html](https://www.hsl.rl.ac.uk/catalogue/hsl_ma48.html) DOVDE seems also have some preconditioned Krylov method variants, such as krysi and vodpk, [https://computing.llnl.gov/projects/odepack](https://computing.llnl.gov/projects/odepack) But I have not figured out how to use those preconditioned Krylov variants, LOL. Are those preconditioned Krylov method good? ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 18, 2022, 10:23am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/30 "2022-01-18T10:23:03Z")
> [@ivanpribec](#): > > A similar mechanism is used also in [RKC (Runge-Kutta-Chebyshev)](http://netlib.org/ode/rkc.f) with which I’m more familiar. The `rkc` subroutine has the interface: > > ```auto > subroutine rkc(neqn,f,y,t,tend,rtol,atol,info,work,idid) > > ``` I have used the rkc solver you mentioned. It seems very fast if setting atol = rtol as its illustration code did! Thanks! ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 18, 2022, 5:18pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/31 "2022-01-18T17:18:02Z")
Thank you @mecej4 ! > [@mecej4](#): > > The F90 version from Byrne and Thompson, of which @CRQuantum has a Gitlab backup, allows the use of a modified version of MA28 from the HSL library. MA28 is itself quite old, and there is a successor, MA48, available (licensed, but free for academic use). About MA48, I just got the f90 version of HSL MA48 files which is the one the F90 version DVODE from Byrne and Thompson can use. I follow the instructions in the F90 DVODE code, but finally I have a problem of getting the Jacobian routines for MA48, as indicated at the end of the F90 DVODE code below, ```auto ! ***** MA48 build change point. Insert the MA48 Jacobian related ! routines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48 here. ! filename = jacobian_for_ma48.f90. Insert the following line ! after the first line of this file: ! USE hsl_ma48_double ``` I do see DVODE have MA28 version of the Jacobian routines such as DVNLSS28, DVJACS28, DVSOLS28. But there is no subroutines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48. **Do you or anyone by any chance know what are those Jacobian subroutines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48?** Thank you very much in advance! PS. In the F90 DVODE, there may be a typo, around line 1882, perhaps it should be ZD11\_TYPE, I did not seem to find ZD01\_TYPE in the MA48 files. ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/a/a2e9d07f6812513b1c05a22bda6522746a2bd6e0.png) ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [January 18, 2022, 6:37pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/32 "2022-01-18T18:37:18Z")
Look near line 12950 of DVODE\_m.f90, where it says “If you have licensed access to the HSL library, an alternate version of DVODE\_F90 based on the successor to MA28, MA48, is available from the authors.” ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 18, 2022, 7:26pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/33 "2022-01-18T19:26:29Z")
Thank you @mecej4 ! I actually have emailed S. Thompson before, but he has not replied yet. Perhaps he was retired? Byrne’s email does not work. May I ask, do you by any chance have a MA48 version of DVODE? You have it, right? LOL. The dvode f90 seems actually good, and with ma48 should be great for sparse Jacobian system 🙂 ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [February 22, 2022, 4:56am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/34 "2022-02-22T04:56:35Z")
Thanks @ivanpribec ! I briefly tested some ODEs by using DVODE and LSODA in the ODEPACK. It seems LSODA is about 2 times faster than DVODE (DVODE is set in stiff ode mode since it cannot automatically switch between stiff and non-stiff) and LSODAR can automatically switch between stiff and non-stiff problem which is also very good. At least for the ODEs I tested, the performance of LSODA is very similar with RKC solver you mentioned before. For DVODE which is basically VODE, it seems in terms of speed, it is not among the most efficient, I found some benchmark [Choosing a Solver — Odes 2.6.2 documentation (scikits-odes.readthedocs.io)](https://scikits-odes.readthedocs.io/en/stable/solvers.html) Although it is in python, I believe the underline code is Fortran or C/C++. So the benchmark is roughly kind of accurate. ---
**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\ **Post date:** [February 22, 2022, 5:27am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/35 "2022-02-22T05:27:30Z")
Some of these Fortran subroutines are wrapped in the R package deSolve. [CRAN - Package deSolve](https://cran.r-project.org/web/packages/deSolve/index.html) ---
**Author:** ![Pap](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pap/32/1305_2.png) [@Pap](https://fortran-lang.discourse.group/u/Pap)\ **Post date:** [February 24, 2022, 8:52am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/36 "2022-02-24T08:52:40Z")
> [@CRquantum](#): > > I found that with Jacobian matrix, there is basically no particular performance gain, but it does improve the DVODE’s stability and ability of reaching convergence. Of course it does… A common programmer’s mistake when using any kind of numerical method that asks for a Jacobian is to omit it, and let the algorithm calculate the Jacobian “automatically” (=numerically). Never do that. Numerical differentiation may seem an easy task but it is _not_; it is actually one of the hardest numerical methods and it is prone to erroneous results. In well-behaving “nice” functions, numerical differentiation _usually_ works well but do not rely on that. Even when it works well, it does so with a specific numerical integration method, which is not necessarily the one implemented in DVODE or whatever else. For example, numerical differentiation via splines may give accurate results with a “not-a-knot” constraint and not so accurate with the popular “natural” spline. Providing a numerical Jacobian to an algorithm that uses it extensively means looking for trouble. It means you give potentially inaccurate information, and that not once or twice, but literally in every single step of the algorithm. You are at the mercy of the specific numerical integration used and the functions involved in the problem at hand: Inaccuracies in the Jacobian (even slight ones) may or may propagate in the next steps of the DVODE (or whenever else algorithm), then you will wonder why it converges slowly or doesn’t converge at all. Whenever you can provide the Jacobian analytically, _do so_. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [February 24, 2022, 9:12am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/37 "2022-02-24T09:12:18Z")
Thank you @Pap ! Nice, now I see the importance of providing Jacobian! Hmmm, talking about Jacobian, if both DVODE and LSODA are given Jacobian, in your opinion which one is faster (if both are solving say, non-stiff problem)? Do you think as @mecej4 pointed, it will be useful to use the MA48 version of dvode? > [@Has anyone used DVODE and is it good?](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/31): > > Thank you @mecej4 ! About MA48, I just got the f90 version of HSL MA48 files which is the one the F90 version DVODE from Byrne and Thompson can use. I follow the instructions in the F90 DVODE code, but finally I have a problem of getting the Jacobian routines for MA48, as indicated at the end of the F90 DVODE code below, ! \*\*\*\*\*MA48 build change point. Insert the MA48 Jacobian related ! routines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48 here. ! filename = jacobian\_for\_ma48.f90. Insert the … Or, does LSODA have something like MA48 already? Thanks! I was looking for a MA48 version of DVODE before but cannot find such a version. I also emailed S. Thompson but perhaps he no longer use his email anymore so no reply 😅 But since LSODA is very decent I guess I just mainly use LOSDA. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [February 24, 2022, 9:21am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/38 "2022-02-24T09:21:26Z")
One more very naïve question @Pap , hmmm, it seems there are some ode solver does not require jacobian, does those solver reliable? Like RKC, [RKC: An explicit solver for parabolic PDEs - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S0377042797002197) which can be found as rkc.f in [ode](https://www.netlib.org/ode/) or the modern ones like FLINT by @prince_mahajan which mainly based on DOP853, [GitHub - princemahajan/FLINT: Fortran Library for numerical INTegration of differential equations](https://github.com/princemahajan/FLINT) ---
**Author:** ![Pap](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pap/32/1305_2.png) [@Pap](https://fortran-lang.discourse.group/u/Pap)\ **Post date:** [February 24, 2022, 10:06am UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/39 "2022-02-24T10:06:15Z")
For non-stiff problems, any decent adaptive step size algorithm should do, and I do not expect a significant performance difference between DVODE and LSODA. However… are you _sure_ your problem is not stiff? In many cases it is easy to tell if a problem is stiff (just a glimpse in the differential equation may smell stiffness ahead). However, it is _not_ so easy to tell if a problem is not stiff. You will be surprised how many problems seem innocent and stiffness-free while they aren’t. Experience in the topic definitely helps in this case. Using an algorithm for stiff problems to solve a non-stiff problem is a bad idea, but if the algorithm is flexible enough, able to switch methods when necessary, then it is worth trying that one first. See what kind of a solution you get; plot the results; do you see any bottlenecks? Any small regions where many integration points are taken? (like the nearly vertical segment in the second figure above). If you don’t see such things, it is relatively safe to assume the problem is not stiff, and you can use DVODE or RKF45 or any other decent method for similar problems; they should perform equally well. And I say “relatively safe” because even the term _stiffness_ doesn’t have a unique, universally adopted definition. Alternatively, use DVODE and see what results you get. If you see nearly horizontal segments with many integration points taken there (like the first figure above), then you are using the wrong method. I must add something here, though. You mention “performance” and “speed” very often, way more often than you mention “convergence”. Sure, it is always nice to have speed, and tons of books dedicate huge sections in algorithm improvements for better performance. But they do so _after_ they assured convergence. if performance is essential (and in many cases is), I would prefer a safe algorithm and will look for ways to parallelize the code, instead of a fast but potentially unsafe algorithm. Solving a “real world” numerical problem is like a gunfight in old wild West: Your first and foremost priority is to make sure you shoot your opponent _then_ to do that fast. You want to shoot him fast, sure, but you want to make sure your bullet goes to the right place first. A slight delay to aim well is better than a fast “intuitive” shoot. Listen to what an old gunslinger says 😆: Fast, inattentive gunslingers won’t see many summers. ---
**Author:** ![oscardssmith](https://avatars.discourse-cdn.com/v4/letter/o/b9e5f3/32.png) [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)\ **Post date:** [February 24, 2022, 2:24pm UTC](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/40 "2022-02-24T14:24:00Z")
Note that modern solvers that use dual number based differentiation (rather than finite difference) don’t require the user to pass in jacobians (and can sometimes be faster by more aggressively combining the code to calculate the primal and jacobian values). [Previous page](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281.md?page=1) [Next page](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281.md?page=3)