Mercurial > hg > octave-lyh

comparison liboctave/NLEqn.cc @ 3:9a4c07481e61

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-08 01:20:23 by jwe] Initial revision
author jwe
date Sun, 08 Aug 1993 01:21:46 +0000
parents
children edfb6cafe85d
comparison
equal deleted inserted replaced
2:c0190df9885d 3:9a4c07481e61
1 // NLEqn.cc -*- C++ -*-
2 /*
3
4 Copyright (C) 1992, 1993 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 #ifdef __GNUG__
25 #pragma implementation
26 #endif
27
28 #include <iostream.h>
29 #include <float.h>
30 #include "NLEqn.h"
31 #include "f77-uscore.h"
32
33 extern "C"
34 {
35 int F77_FCN (hybrd1) (int (*)(), const int*, double*, double*,
36 const double*, int*, double*, const int*);
37
38 int F77_FCN (hybrj1) (int (*)(), const int*, double*, double*,
39 double*, const int*, const double*, int*,
40 double*, const int*);
41 }
42
43 static nonlinear_fcn user_fun;
44 static jacobian_fcn user_jac;
45
46 // error handling
47
48 void
49 NLEqn::error (const char* msg)
50 {
51 cerr << "Fatal NLEqn error. " << msg << "\n";
52 exit(1);
53 }
54
55 // Constructors
56
57 NLEqn::NLEqn (void) : NLFunc (), x (), n (0) {}
58
59 NLEqn::NLEqn (const Vector& xvec, const NLFunc f)
60 : NLFunc (f), x (xvec), n (x.capacity ()) {}
61
62 NLEqn::NLEqn (const NLEqn& a) : NLFunc (a.fun, a.jac), x (a.x), n (a.n) {}
63
64 void
65 NLEqn::resize (int nn)
66 {
67 if (n != nn)
68 {
69 n = nn;
70 x.resize (n);
71 }
72 }
73
74 int
75 NLEqn::size (void) const
76 {
77 return n;
78 }
79
80 // Assignment
81
82 NLEqn&
83 NLEqn::operator = (const NLEqn& a)
84 {
85 fun = a.fun;
86 jac = a.jac;
87 x = a.n;
88
89 return *this;
90 }
91
92 Vector
93 NLEqn::states (void) const
94 {
95 return x;
96 }
97
98 void
99 NLEqn::set_states (const Vector& xvec)
100 {
101 if (xvec.capacity () != n)
102 error ("dimension error");
103
104 x = xvec;
105 }
106
107 // Other operations
108
109 Vector
110 NLEqn::solve (const Vector& xvec)
111 {
112 set_states (xvec);
113 int info;
114 return solve (info);
115 }
116
117 Vector
118 NLEqn::solve (const Vector& xvec, int& info)
119 {
120 set_states (xvec);
121 return solve (info);
122 }
123
124 Vector
125 NLEqn::solve (void)
126 {
127 int info;
128 return solve (info);
129 }
130
131 int
132 hybrd1_fcn (int *n, double *x, double *fvec, int *iflag)
133 {
134 int nn = *n;
135 Vector tmp_f (nn);
136 Vector tmp_x (nn);
137
138 for (int i = 0; i < nn; i++)
139 tmp_x.elem (i) = x[i];
140
141 tmp_f = (*user_fun) (tmp_x);
142
143 for (i = 0; i < nn; i++)
144 fvec[i] = tmp_f.elem (i);
145
146 return 0;
147 }
148
149 int
150 hybrj1_fcn (int *n, double *x, double *fvec, double *fjac,
151 int *ldfjac, int *iflag)
152 {
153 int nn = *n;
154 Vector tmp_x (nn);
155
156 for (int i = 0; i < nn; i++)
157 tmp_x.elem (i) = x[i];
158
159 int flag = *iflag;
160 if (flag == 1)
161 {
162 Vector tmp_f (nn);
163
164 tmp_f = (*user_fun) (tmp_x);
165
166 for (i = 0; i < nn; i++)
167 fvec[i] = tmp_f.elem (i);
168 }
169 else
170 {
171 Matrix tmp_fj (nn, nn);
172
173 tmp_fj = (*user_jac) (tmp_x);
174
175 int ld = *ldfjac;
176 for (int j = 0; j < nn; j++)
177 for (i = 0; i < nn; i++)
178 fjac[j*ld+i] = tmp_fj.elem (i, j);
179 }
180
181 return 0;
182 }
183
184 Vector
185 NLEqn::solve (int& info)
186 {
187 int tmp_info = 0;
188
189 if (n == 0)
190 error ("Equation set not initialized");
191
192 double tol = sqrt (DBL_EPSILON);
193
194 double *fvec = new double [n];
195 double *px = new double [n];
196 for (int i = 0; i < n; i++)
197 px[i] = x.elem (i);
198
199 user_fun = fun;
200 user_jac = jac;
201
202 if (jac == NULL)
203 {
204 int lwa = (n*(3*n+13))/2;
205 double *wa = new double [lwa];
206
207 F77_FCN (hybrd1) (hybrd1_fcn, &n, px, fvec, &tol, &tmp_info, wa, &lwa);
208
209 delete [] wa;
210 }
211 else
212 {
213 int lwa = (n*(n+13))/2;
214 double *wa = new double [lwa];
215 double *fjac = new double [n*n];
216
217 F77_FCN (hybrj1) (hybrj1_fcn, &n, px, fvec, fjac, &n, &tol,
218 &tmp_info, wa, &lwa);
219
220 delete [] wa;
221 delete [] fjac;
222 }
223
224 info = tmp_info;
225
226 Vector retval (n);
227
228 for (i = 0; i < n; i++)
229 retval.elem (i) = px[i];
230
231 return retval;
232 }
233
234 /*
235 ;;; Local Variables: ***
236 ;;; mode: C++ ***
237 ;;; page-delimiter: "^/\\*" ***
238 ;;; End: ***
239 */