Line data Source code
1 : module test_int_support
2 : use num_def
3 : use num_lib
4 : use utils_lib, only: mesa_error
5 :
6 : implicit none
7 :
8 : integer, parameter :: ipar_sparse_format = 1
9 : ! =0 means compressed row format; else, compressed column format.
10 : integer, parameter :: i_nfcn = 2
11 : integer, parameter :: i_njac = 3
12 :
13 : contains
14 :
15 29 : subroutine do_test_stiff_int( &
16 : which_solver, which_decsol, numerical_jacobian, &
17 : fcn, jac, sjac, solout, iout_input, &
18 : fcn_blk_dble, jac_blk_dble, &
19 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
20 : n, ndisc, mljac, mujac, matrix_type_spec, &
21 29 : mas, imas, mlmas, mumas, m1, m2, t, &
22 : rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
23 : use test_support, only: show_results, show_statistics
24 : use mtx_lib
25 : use mtx_def
26 : integer, intent(in) :: which_solver, which_decsol
27 : interface
28 : include 'num_fcn.dek'
29 : include 'num_fcn_blk_dble.dek'
30 : include 'num_jac.dek'
31 : include 'num_jac_blk_dble.dek'
32 : include 'num_sjac.dek'
33 : include 'num_solout.dek'
34 : include 'num_mas.dek'
35 : end interface
36 : integer :: caller_id, nvar_blk_dble, nz_blk_dble
37 : real(dp), dimension(:), pointer :: lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
38 : integer, intent(in) :: imas, mlmas, mumas, m1, m2, iout_input
39 : integer, intent(in) :: n, ndisc, mljac, mujac, matrix_type_spec, lrpar, lipar, itol
40 : logical, intent(in) :: numerical_jacobian, quiet
41 : real(dp), intent(inout) :: t(0:ndisc + 1), rtol(*), atol(*), h0
42 : real(dp), pointer :: y(:) ! (n)
43 : integer, intent(inout) :: nstep
44 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
45 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
46 : integer, intent(out) :: ierr
47 :
48 : integer :: i, lout, iout, idid, ijac, max_cols_exptrap, max_steps, nrdens, &
49 : nzmax, lrd, lid, isparse, liwork, lwork
50 29 : real(dp) :: h, eps, max_step_size, y_min, y_max
51 :
52 29 : integer, pointer :: iwork(:) !(liwork)
53 29 : real(dp), pointer :: work(:) !(lwork)
54 29 : integer, pointer :: ipar_decsol(:) !(lid)
55 29 : real(dp), pointer :: rpar_decsol(:) !(lrd)
56 :
57 29 : iout = iout_input
58 29 : if (quiet) iout = 0
59 29 : max_steps = 500000
60 29 : max_step_size = 0
61 29 : isparse = 0
62 29 : lout = 6
63 :
64 29 : y_min = -1d199
65 29 : y_max = 1d199
66 :
67 29 : if (numerical_jacobian) then
68 11 : ijac = 0
69 : else
70 18 : ijac = 1
71 : end if
72 :
73 116 : ipar = 0
74 58 : rpar = 0
75 :
76 29 : nrdens = n
77 29 : max_cols_exptrap = 0 ! use default
78 :
79 29 : lid = 0; lrd = 0
80 29 : if (which_decsol == lapack) then
81 29 : nzmax = 0
82 : call lapack_work_sizes(n, lrd, lid)
83 : else
84 0 : if (mljac == n) then
85 0 : nzmax = n*n
86 : else
87 0 : nzmax = n*(mljac + mujac + 1)
88 : end if
89 0 : write (*, *) 'test_int_support: bad which_decsol', which_decsol
90 0 : call mesa_error(__FILE__, __LINE__) ! test_int_support
91 : end if
92 :
93 : call isolve_work_sizes(n, nzmax, imas, mljac, mujac, mlmas, mumas, liwork, lwork)
94 :
95 29 : allocate (iwork(liwork), work(lwork), ipar_decsol(lid), rpar_decsol(lrd), stat=ierr)
96 29 : if (ierr /= 0) then
97 0 : write (*, *) 'allocate ierr', ierr
98 0 : call mesa_error(__FILE__, __LINE__) ! test_int_support
99 : end if
100 :
101 7766 : iwork = 0
102 94991 : work = 0
103 :
104 29 : iwork(9) = m1
105 29 : iwork(10) = m2
106 :
107 29 : nstep = 0
108 29 : eps = rtol(1)
109 66 : do i = 0, ndisc
110 37 : ierr = 0
111 37 : h = h0
112 40 : select case (which_solver)
113 : case (ros2_solver)
114 3 : if (i == 0 .and. .not. quiet) write (*, *) 'ros2'
115 : case (rose2_solver)
116 3 : if (i == 0 .and. .not. quiet) write (*, *) 'rose2'
117 : case (ros3p_solver)
118 3 : if (i == 0 .and. .not. quiet) write (*, *) 'ros3p'
119 : case (ros3pl_solver)
120 7 : if (i == 0 .and. .not. quiet) write (*, *) 'ros3pl'
121 : case (rodas3_solver)
122 7 : if (i == 0 .and. .not. quiet) write (*, *) 'rodas3'
123 : case (rodas4_solver)
124 7 : if (i == 0 .and. .not. quiet) write (*, *) 'rodas4'
125 : case (rodasp_solver)
126 7 : if (i == 0 .and. .not. quiet) write (*, *) 'rodasp'
127 : case default
128 0 : write (*, *) 'unknown value for which_solver'
129 37 : call mesa_error(__FILE__, __LINE__)
130 : end select
131 :
132 37 : if (which_decsol == lapack) then
133 37 : if (i == 0 .and. .not. quiet) write (*, *) 'lapack_decsol'
134 37 : call do_isolve(lapack_decsol, null_decsols, null_decsolblk)
135 : else
136 0 : write (*, *) 'unknown value for which_decsol', which_decsol
137 0 : call mesa_error(__FILE__, __LINE__)
138 : end if
139 :
140 37 : if (idid /= 1) ierr = -1
141 37 : if (ierr /= 0) then
142 0 : write (*, *) 'solver returned ierr /= 0', idid
143 0 : call mesa_error(__FILE__, __LINE__)
144 : end if
145 66 : nstep = nstep + iwork(16) ! nsteps
146 : end do
147 :
148 58 : deallocate (iwork, work, ipar_decsol, rpar_decsol)
149 :
150 : contains
151 :
152 37 : subroutine do_isolve(decsol, decsols, decsolblk)
153 : interface
154 : include "mtx_decsol.dek"
155 : include "mtx_decsols.dek"
156 : include "mtx_decsolblk.dek"
157 : end interface
158 : integer :: j
159 : include 'formats'
160 : call isolve( &
161 37 : which_solver, n, fcn, t(i), y, t(i + 1), &
162 : h, max_step_size, max_steps, &
163 : rtol, atol, itol, y_min, y_max, &
164 : jac, ijac, sjac, nzmax, isparse, mljac, mujac, &
165 : mas, imas, mlmas, mumas, &
166 : solout, iout, &
167 : decsol, decsols, decsolblk, &
168 : lrd, rpar_decsol, lid, ipar_decsol, &
169 : caller_id, nvar_blk_dble, nz_blk_dble, &
170 : lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
171 : fcn_blk_dble, jac_blk_dble, &
172 : work, lwork, iwork, liwork, &
173 : lrpar, rpar, lipar, ipar, &
174 74 : lout, idid)
175 : return
176 : do j = 1, n
177 : write (*, 2) 'y(j)', j, y(j)
178 : end do
179 29 : end subroutine do_isolve
180 :
181 : end subroutine do_test_stiff_int
182 :
183 : end module test_int_support
|
