Mercurial > octave

comparison liboctave/numeric/LSODE.cc @ 31249:de6fc38c78c6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Jacobian types offered by dlsode.f accessible by lsode (bug #31626). * liboctave/numeric/LSODE-opts.in: Add options "jacobian type", "lower jacobian subdiagonals", and "upper jacobian subdiagonals". * liboctave/numeric/LSODE.cc (file scope, lsode_j, LSODE::do_integrate (double)): Handle new configurable Jacobian types. * build-aux/mk-opts.pl: Don't implicitly convert to integer in condition.
author Olaf Till <olaf.till@uni-jena.de>
date Fri, 12 Nov 2010 08:53:05 +0100
parents 51a3d3a69193
children a40c0b7aa376
comparison
equal deleted inserted replaced
31248:8b75954a4670 31249:de6fc38c78c6
53 } 53 }
54 54
55 static ODEFunc::ODERHSFunc user_fcn; 55 static ODEFunc::ODERHSFunc user_fcn;
56 static ODEFunc::ODEJacFunc user_jac; 56 static ODEFunc::ODEJacFunc user_jac;
57 static ColumnVector *tmp_x; 57 static ColumnVector *tmp_x;
58 static bool user_jac_ignore_ml_mu;
58 59
59 static F77_INT 60 static F77_INT
60 lsode_f (const F77_INT& neq, const double& time, double *, double *deriv, 61 lsode_f (const F77_INT& neq, const double& time, double *, double *deriv,
61 F77_INT& ierr) 62 F77_INT& ierr)
62 { 63 {
78 79
79 return 0; 80 return 0;
80 } 81 }
81 82
82 static F77_INT 83 static F77_INT
83 lsode_j (const F77_INT& neq, const double& time, double *, const F77_INT&, 84 lsode_j (const F77_INT& neq, const double& time, double *,
84 const F77_INT&, double *pd, const F77_INT& nrowpd) 85 const F77_INT& ml, const F77_INT& mu,
86 double *pd, const F77_INT& nrowpd)
85 { 87 {
86 Matrix tmp_jac (neq, neq); 88 Matrix tmp_jac (neq, neq);
87 89
88 // NOTE: this won't work if LSODE passes copies of the state vector. 90 // NOTE: this won't work if LSODE passes copies of the state vector.
89 // In that case we have to create a temporary vector object 91 // In that case we have to create a temporary vector object
90 // and copy. 92 // and copy.
91 93
92 tmp_jac = (*user_jac) (*tmp_x, time); 94 tmp_jac = (*user_jac) (*tmp_x, time);
93 95
94 for (F77_INT j = 0; j < neq; j++) 96 if (user_jac_ignore_ml_mu)
95 for (F77_INT i = 0; i < neq; i++) 97 for (F77_INT j = 0; j < neq; j++)
96 pd[nrowpd * j + i] = tmp_jac (i, j); 98 for (F77_INT i = 0; i < neq; i++)
97 99 pd[nrowpd * j + i] = tmp_jac (i, j);
100 else
101 // upper left ends of subdiagonals in tmp_jac
102 for (F77_INT i = 0, j = mu; i <= ml; j == 0 ? i++ : j--)
103 // walk down the subdiagonal in tmp_jac
104 for (F77_INT k = i, l = j; k < neq && l < neq; k++, l++)
105 pd[nrowpd * l + k + mu - l] = tmp_jac (k, l);
106
98 return 0; 107 return 0;
99 } 108 }
100 109
101 ColumnVector 110 ColumnVector
102 LSODE::do_integrate (double tout) 111 LSODE::do_integrate (double tout)
118 127
119 nn = n; 128 nn = n;
120 129
121 octave_idx_type max_maxord = 0; 130 octave_idx_type max_maxord = 0;
122 131
132 user_jac_ignore_ml_mu = true;
133
134 m_iwork(0) = lower_jacobian_subdiagonals (); // 'ML' in dlsode.f
135
136 m_iwork(1) = upper_jacobian_subdiagonals (); // 'MU' in dlsode.f
137
123 if (integration_method () == "stiff") 138 if (integration_method () == "stiff")
124 { 139 {
125 max_maxord = 5; 140 max_maxord = 5;
126 141
127 if (m_jac) 142 if (m_jac)
128 m_method_flag = 21; 143 {
144 if (jacobian_type () == "banded")
145 {
146 m_method_flag = 24;
147 user_jac_ignore_ml_mu = false;
148 }
149 else
150 m_method_flag = 21;
151 }
129 else 152 else
130 m_method_flag = 22; 153 {
154 if (jacobian_type () == "full")
155 m_method_flag = 22;
156 else if (jacobian_type () == "banded")
157 m_method_flag = 25;
158 else if (jacobian_type () == "diagonal")
159 m_method_flag = 23;
160 else
161 {
162 // should be prevented by lsode_options
163 (*current_liboctave_error_handler)
164 ("lsode: internal error, wrong jacobian type");
165 m_integration_error = true;
166 return retval;
167 }
168 }
131 169
132 m_liw = 20 + n; 170 m_liw = 20 + n;
133 m_lrw = 22 + n * (9 + n); 171 m_lrw = 22 + n * (9 + n);
134 } 172 }
135 else 173 else