Mercurial > hg > octave-lyh

comparison liboctave/DASRT.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 53b4eab68976
comparison
equal deleted inserted replaced
3990:46388d6a4e44 3991:48d2bc4a3729
43 #include "parse.h" 43 #include "parse.h"
44 #include "unwind-prot.h" 44 #include "unwind-prot.h"
45 #include "utils.h" 45 #include "utils.h"
46 #include "variables.h" 46 #include "variables.h"
47 47
48 // For instantiating the Array<Matrix> object.
49 #include "Array.h"
50 #include "Array.cc"
51
52 #include "DASRT.h" 48 #include "DASRT.h"
53 #include "f77-fcn.h" 49 #include "f77-fcn.h"
54 #include "lo-error.h" 50 #include "lo-error.h"
55 51
56 #ifndef F77_FUNC 52 #ifndef F77_FUNC
57 #define F77_FUNC(x, X) F77_FCN (x, X) 53 #define F77_FUNC(x, X) F77_FCN (x, X)
58 #endif 54 #endif
59 55
56 typedef int (*dasrt_fcn_ptr) (const double&, const double*, const double*,
57 double*, int&, double*, int*);
58
59 typedef int (*dasrt_jac_ptr) (const double&, const double*, const double*,
60 double*, const double&, double*, int*);
61
62 typedef int (*dasrt_constr_ptr) (const int&, const double&, const double*,
63 const int&, double*, double*, int*);
64
60 extern "C" 65 extern "C"
61 { 66 int F77_FUNC (ddasrt, DASRT) (dasrt_fcn_ptr, const int&, double&,
62 int F77_FUNC (ddasrt, DASRT) (int (*)(const double&, double*, double*, 67 double*, double*, const double&, int*,
63 double*, int&, double*, int*), 68 double*, double*, int&, double*,
64 const int&, const double&, double*, double*, 69 const int&, int*, const int&, double*,
65 const double&, int*, double*, double*, 70 int*, dasrt_jac_ptr, dasrt_constr_ptr,
66 int&, double*, const int&, int*, 71 const int&, int*);
67 const int&, double*, int*,
68 int (*)(const double&, double*,
69 double*, double*,
70 const double&, double*, int*),
71 int (*)(const int&, const double&, double*,
72 const int&, double*, double*, int*),
73 const int&, int*);
74 }
75
76 template class Array<Matrix>;
77 72
78 static DAEFunc::DAERHSFunc user_fsub; 73 static DAEFunc::DAERHSFunc user_fsub;
79 static DAEFunc::DAEJacFunc user_jsub; 74 static DAEFunc::DAEJacFunc user_jsub;
80 static DAERTFunc::DAERTConstrFunc user_csub; 75 static DAERTFunc::DAERTConstrFunc user_csub;
81 static int nn; 76 static int nn;
82 77
83 static int 78 static int
84 ddasrt_f (const double& t, double *state, double *deriv, double *delta, 79 ddasrt_f (const double& t, const double *state, const double *deriv,
85 int& ires, double *rpar, int *ipar) 80 double *delta, int& ires, double *rpar, int *ipar)
86 { 81 {
87 ColumnVector tmp_state (nn); 82 ColumnVector tmp_state (nn);
88 for (int i = 0; i < nn; i++) 83 for (int i = 0; i < nn; i++)
89 tmp_state(i) = state[i]; 84 tmp_state(i) = state[i];
90 85
109 // double *ework, double *rpar, int *ipar, 104 // double *ework, double *rpar, int *ipar,
110 // const int& ieform, int& ires); 105 // const int& ieform, int& ires);
111 106
112 //static efptr e_fun; 107 //static efptr e_fun;
113 108
114 static int 109 int
115 ddasrt_j (const double& t, double *state, double *deriv, 110 ddasrt_j (const double& time, const double *state, const double *deriv,
116 double *pdwork, const double& cj, double *rpar, int *ipar) 111 double *pd, const double& cj, double *, int *)
117 { 112 {
113 // XXX FIXME XXX -- would be nice to avoid copying the data.
114
118 ColumnVector tmp_state (nn); 115 ColumnVector tmp_state (nn);
116 ColumnVector tmp_deriv (nn);
117
119 for (int i = 0; i < nn; i++) 118 for (int i = 0; i < nn; i++)
120 tmp_state(i) = state[i]; 119 {
121 120 tmp_deriv.elem (i) = deriv [i];
122 ColumnVector tmp_deriv (nn); 121 tmp_state.elem (i) = state [i];
123 for (int i = 0; i < nn; i++) 122 }
124 tmp_deriv(i) = deriv[i]; 123
125 124 Matrix tmp_pd = user_jsub (tmp_state, tmp_deriv, time, cj);
126 // XXX FIXME XXX
127
128 Matrix tmp_dfdxdot (nn, nn);
129 Matrix tmp_dfdx (nn, nn);
130
131 DAEFunc::DAEJac tmp_jac;
132 tmp_jac.dfdxdot = &tmp_dfdxdot;
133 tmp_jac.dfdx = &tmp_dfdx;
134
135 tmp_jac = user_jsub (tmp_state, tmp_deriv, t);
136
137 // Fix up the matrix of partial derivatives for dasrt.
138
139 tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
140 125
141 for (int j = 0; j < nn; j++) 126 for (int j = 0; j < nn; j++)
142 for (int i = 0; i < nn; i++) 127 for (int i = 0; i < nn; i++)
143 pdwork[j*nn+i] = tmp_dfdx.elem (i, j); 128 pd [nn * j + i] = tmp_pd.elem (i, j);
144 129
145 return 0; 130 return 0;
146 } 131 }
147 132
148 static int 133 static int
149 ddasrt_g (const int& neq, const double& t, double *state, const int& ng, 134 ddasrt_g (const int& neq, const double& t, const double *state,
150 double *gout, double *rpar, int *ipar) 135 const int& ng, double *gout, double *rpar, int *ipar)
151 { 136 {
152 int n = neq; 137 int n = neq;
153 138
154 ColumnVector tmp_state (n); 139 ColumnVector tmp_state (n);
155 for (int i = 0; i < n; i++) 140 for (int i = 0; i < n; i++)