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