Mercurial > hg > octave-nkf

comparison liboctave/DASSL.cc @ 3991:48d2bc4a3729

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-07-16 17:46:50 by jwe]
author jwe
date Tue, 16 Jul 2002 17:46:51 +0000
parents 46388d6a4e44
children f23bc69132cc
comparison
equal deleted inserted replaced
3990:46388d6a4e44 3991:48d2bc4a3729
33 33
34 #include "DASSL.h" 34 #include "DASSL.h"
35 #include "f77-fcn.h" 35 #include "f77-fcn.h"
36 #include "lo-error.h" 36 #include "lo-error.h"
37 37
38 typedef int (*dassl_fcn_ptr) (const double&, double*, double*, 38 typedef int (*dassl_fcn_ptr) (const double&, const double*, const double*,
39 double*, int&, double*, int*); 39 double*, int&, double*, int*);
40 40
41 typedef int (*dassl_jac_ptr) (const double&, double*, double*, 41 typedef int (*dassl_jac_ptr) (const double&, const double*, const double*,
42 double*, const double&, double*, int*); 42 double*, const double&, double*, int*);
43 43
44 extern "C" 44 extern "C"
45 int F77_FUNC (ddassl, DDASSL) (dassl_fcn_ptr, const int&, double&, 45 int F77_FUNC (ddassl, DDASSL) (dassl_fcn_ptr, const int&, double&,
46 double*, double*, double&, const int*, 46 double*, double*, double&, const int*,
47 const double&, const double&, int&, 47 const double&, const double&, int&,
48 double*, const int&, int*, const int&, 48 double*, const int&, int*, const int&,
49 const double*, const int*, 49 const double*, const int*,
50 dassl_jac_ptr); 50 dassl_jac_ptr);
51 51
52 static DAEFunc::DAERHSFunc user_fun; 52 static DAEFunc::DAERHSFunc user_fun;
53 static DAEFunc::DAEJacFunc user_jac; 53 static DAEFunc::DAEJacFunc user_jac;
54 static int nn; 54 static int nn;
55 55
130 { 130 {
131 stop_time_set = 0; 131 stop_time_set = 0;
132 } 132 }
133 133
134 int 134 int
135 ddassl_f (const double& time, double *state, double *deriv, 135 ddassl_f (const double& time, const double *state, const double *deriv,
136 double *delta, int& ires, double *, int *) 136 double *delta, int& ires, double *, int *)
137 { 137 {
138 // XXX FIXME XXX -- would be nice to avoid copying the data.
139
138 ColumnVector tmp_deriv (nn); 140 ColumnVector tmp_deriv (nn);
139 ColumnVector tmp_state (nn); 141 ColumnVector tmp_state (nn);
140 ColumnVector tmp_delta (nn); 142 ColumnVector tmp_delta (nn);
141 143
142 for (int i = 0; i < nn; i++) 144 for (int i = 0; i < nn; i++)
160 162
161 return 0; 163 return 0;
162 } 164 }
163 165
164 int 166 int
165 ddassl_j (const double& time, double *, double *, double *pd, const 167 ddassl_j (const double& time, const double *state, const double *deriv,
166 double& cj, double *, int *) 168 double *pd, const double& cj, double *, int *)
167 { 169 {
170 // XXX FIXME XXX -- would be nice to avoid copying the data.
171
168 ColumnVector tmp_state (nn); 172 ColumnVector tmp_state (nn);
169 ColumnVector tmp_deriv (nn); 173 ColumnVector tmp_deriv (nn);
170 174
171 // XXX FIXME XXX 175 for (int i = 0; i < nn; i++)
172 176 {
173 Matrix tmp_dfdxdot (nn, nn); 177 tmp_deriv.elem (i) = deriv [i];
174 Matrix tmp_dfdx (nn, nn); 178 tmp_state.elem (i) = state [i];
175 179 }
176 DAEFunc::DAEJac tmp_jac; 180
177 tmp_jac.dfdxdot = &tmp_dfdxdot; 181 Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
178 tmp_jac.dfdx = &tmp_dfdx;
179
180 tmp_jac = user_jac (tmp_state, tmp_deriv, time);
181
182 // Fix up the matrix of partial derivatives for dassl.
183
184 tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
185 182
186 for (int j = 0; j < nn; j++) 183 for (int j = 0; j < nn; j++)
187 for (int i = 0; i < nn; i++) 184 for (int i = 0; i < nn; i++)
188 pd [nn * j + i] = tmp_dfdx.elem (i, j); 185 pd [nn * j + i] = tmp_pd.elem (i, j);
189 186
190 return 0; 187 return 0;
191 } 188 }
192 189
193 ColumnVector 190 ColumnVector