← Back to branch summary

~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to liboctave/LSODE.cc

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
Import upstream version 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
liboctave/LSODE.cc
 
1
/*
 
2
 
 
3
Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2002,
 
4
              2003, 2004, 2005, 2006, 2007 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 3 of the License, or (at your
 
11
option) any 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, see
 
20
<http://www.gnu.org/licenses/>.
 
21
 
 
22
*/
 
23
 
 
24
#ifdef HAVE_CONFIG_H
 
25
#include <config.h>
 
26
#endif
 
27
 
 
28
#include <cfloat>
 
29
 
 
30
#include <sstream>
 
31
 
 
32
#include "LSODE.h"
 
33
#include "f77-fcn.h"
 
34
#include "lo-error.h"
 
35
#include "lo-math.h"
 
36
#include "quit.h"
 
37
 
 
38
typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, const double&, double*,
 
39
                              double*, octave_idx_type&);
 
40
 
 
41
typedef octave_idx_type (*lsode_jac_ptr) (const octave_idx_type&, const double&, double*,
 
42
                              const octave_idx_type&, const octave_idx_type&, double*, const
 
43
                              octave_idx_type&);
 
44
 
 
45
extern "C"
 
46
{
 
47
  F77_RET_T
 
48
  F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, double&,
 
49
                             double&, octave_idx_type&, double&, const double*, octave_idx_type&,
 
50
                             octave_idx_type&, octave_idx_type&, double*, octave_idx_type&, octave_idx_type*, octave_idx_type&,
 
51
                             lsode_jac_ptr, octave_idx_type&);
 
52
}
 
53
 
 
54
static ODEFunc::ODERHSFunc user_fun;
 
55
static ODEFunc::ODEJacFunc user_jac;
 
56
static ColumnVector *tmp_x;
 
57
 
 
58
static octave_idx_type
 
59
lsode_f (const octave_idx_type& neq, const double& time, double *,
 
60
         double *deriv, octave_idx_type& ierr) 
 
61
{
 
62
  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
 
63
 
 
64
  ColumnVector tmp_deriv;
 
65
 
 
66
  // NOTE: this won't work if LSODE passes copies of the state vector.
 
67
  //       In that case we have to create a temporary vector object
 
68
  //       and copy.
 
69
 
 
70
  tmp_deriv = (*user_fun) (*tmp_x, time);
 
71
 
 
72
  if (tmp_deriv.length () == 0)
 
73
    ierr = -1;
 
74
  else
 
75
    {
 
76
      for (octave_idx_type i = 0; i < neq; i++)
 
77
        deriv [i] = tmp_deriv.elem (i);
 
78
    }
 
79
 
 
80
  END_INTERRUPT_WITH_EXCEPTIONS;
 
81
 
 
82
  return 0;
 
83
}
 
84
 
 
85
static octave_idx_type
 
86
lsode_j (const octave_idx_type& neq, const double& time, double *,
 
87
         const octave_idx_type&, const octave_idx_type&, double *pd, const octave_idx_type& nrowpd)
 
88
{
 
89
  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
 
90
 
 
91
  Matrix tmp_jac (neq, neq);
 
92
 
 
93
  // NOTE: this won't work if LSODE passes copies of the state vector.
 
94
  //       In that case we have to create a temporary vector object
 
95
  //       and copy.
 
96
 
 
97
  tmp_jac = (*user_jac) (*tmp_x, time);
 
98
 
 
99
  for (octave_idx_type j = 0; j < neq; j++)
 
100
    for (octave_idx_type i = 0; i < neq; i++)
 
101
      pd [nrowpd * j + i] = tmp_jac (i, j);
 
102
 
 
103
  END_INTERRUPT_WITH_EXCEPTIONS;
 
104
 
 
105
  return 0;
 
106
}
 
107
 
 
108
ColumnVector
 
109
LSODE::do_integrate (double tout)
 
110
{
 
111
  ColumnVector retval;
 
112
 
 
113
  static octave_idx_type nn = 0;
 
114
 
 
115
  if (! initialized || restart || ODEFunc::reset || LSODE_options::reset)
 
116
    {
 
117
      integration_error = false;
 
118
 
 
119
      initialized = true;
 
120
 
 
121
      istate = 1;
 
122
 
 
123
      octave_idx_type n = size ();
 
124
 
 
125
      nn = n;
 
126
 
 
127
      octave_idx_type max_maxord = 0;
 
128
 
 
129
      if (integration_method () == "stiff")
 
130
        {
 
131
          max_maxord = 5;
 
132
 
 
133
          if (jac)
 
134
            method_flag = 21;
 
135
          else
 
136
            method_flag = 22;
 
137
 
 
138
          liw = 20 + n;
 
139
          lrw = 22 + n * (9 + n);
 
140
        }
 
141
      else
 
142
        {
 
143
          max_maxord = 12;
 
144
 
 
145
          method_flag = 10;
 
146
 
 
147
          liw = 20;
 
148
          lrw = 22 + 16 * n;
 
149
        }
 
150
 
 
151
      maxord = maximum_order ();
 
152
 
 
153
      iwork.resize (liw);
 
154
 
 
155
      for (octave_idx_type i = 4; i < 9; i++)
 
156
        iwork(i) = 0;
 
157
 
 
158
      rwork.resize (lrw);
 
159
 
 
160
      for (octave_idx_type i = 4; i < 9; i++)
 
161
        rwork(i) = 0;
 
162
 
 
163
      if (maxord >= 0)
 
164
        {
 
165
          if (maxord > 0 && maxord <= max_maxord)
 
166
            {
 
167
              iwork(4) = maxord;
 
168
              iopt = 1;
 
169
            }     
 
170
          else
 
171
            {
 
172
              (*current_liboctave_error_handler)
 
173
                ("lsode: invalid value for maximum order");
 
174
              integration_error = true;
 
175
              return retval;
 
176
            }
 
177
        }
 
178
 
 
179
      if (stop_time_set)
 
180
        {
 
181
          itask = 4;
 
182
          rwork(0) = stop_time;
 
183
          iopt = 1;
 
184
        }
 
185
      else
 
186
        {
 
187
          itask = 1;
 
188
        }
 
189
 
 
190
      px = x.fortran_vec ();
 
191
 
 
192
      piwork = iwork.fortran_vec ();
 
193
      prwork = rwork.fortran_vec ();
 
194
 
 
195
      restart = false;
 
196
 
 
197
      // ODEFunc
 
198
 
 
199
      // NOTE: this won't work if LSODE passes copies of the state vector.
 
200
      //       In that case we have to create a temporary vector object
 
201
      //       and copy.
 
202
 
 
203
      tmp_x = &x;
 
204
 
 
205
      user_fun = function ();
 
206
      user_jac = jacobian_function ();
 
207
 
 
208
      ColumnVector xdot = (*user_fun) (x, t);
 
209
 
 
210
      if (x.length () != xdot.length ())
 
211
        {
 
212
          (*current_liboctave_error_handler)
 
213
            ("lsode: inconsistent sizes for state and derivative vectors");
 
214
 
 
215
          integration_error = true;
 
216
          return retval;
 
217
        }
 
218
 
 
219
      ODEFunc::reset = false;
 
220
 
 
221
      // LSODE_options
 
222
 
 
223
      rel_tol = relative_tolerance ();
 
224
      abs_tol = absolute_tolerance ();
 
225
 
 
226
      octave_idx_type abs_tol_len = abs_tol.length ();
 
227
 
 
228
      if (abs_tol_len == 1)
 
229
        itol = 1;
 
230
      else if (abs_tol_len == n)
 
231
        itol = 2;
 
232
      else
 
233
        {
 
234
          (*current_liboctave_error_handler)
 
235
            ("lsode: inconsistent sizes for state and absolute tolerance vectors");
 
236
 
 
237
          integration_error = true;
 
238
          return retval;
 
239
        }
 
240
 
 
241
      double iss = initial_step_size ();
 
242
      if (iss >= 0.0)
 
243
        {
 
244
          rwork(4) = iss;
 
245
          iopt = 1;
 
246
        }
 
247
 
 
248
      double maxss = maximum_step_size ();
 
249
      if (maxss >= 0.0)
 
250
        {
 
251
          rwork(5) = maxss;
 
252
          iopt = 1;
 
253
        }
 
254
 
 
255
      double minss = minimum_step_size ();
 
256
      if (minss >= 0.0)
 
257
        {
 
258
          rwork(6) = minss;
 
259
          iopt = 1;
 
260
        }
 
261
 
 
262
      octave_idx_type sl = step_limit ();
 
263
      if (sl > 0)
 
264
        {
 
265
          iwork(5) = sl;
 
266
          iopt = 1;
 
267
        }
 
268
 
 
269
      pabs_tol = abs_tol.fortran_vec ();
 
270
 
 
271
      LSODE_options::reset = false;
 
272
    }
 
273
 
 
274
  F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol,
 
275
                             pabs_tol, itask, istate, iopt, prwork, lrw,
 
276
                             piwork, liw, lsode_j, method_flag));
 
277
 
 
278
  if (f77_exception_encountered)
 
279
    {
 
280
      integration_error = true;
 
281
      (*current_liboctave_error_handler) ("unrecoverable error in lsode");
 
282
    }
 
283
  else
 
284
    {
 
285
      switch (istate)
 
286
        {
 
287
        case 1:  // prior to initial integration step.
 
288
        case 2:  // lsode was successful.
 
289
          retval = x;
 
290
          t = tout;
 
291
          break;
 
292
          
 
293
        case -1:  // excess work done on this call (perhaps wrong mf).
 
294
        case -2:  // excess accuracy requested (tolerances too small).
 
295
        case -3:  // illegal input detected (see printed message).
 
296
        case -4:  // repeated error test failures (check all inputs).
 
297
        case -5:  // repeated convergence failures (perhaps bad jacobian
 
298
                  // supplied or wrong choice of mf or tolerances).
 
299
        case -6:  // error weight became zero during problem. (solution
 
300
                  // component i vanished, and atol or atol(i) = 0.)
 
301
        case -13: // return requested in user-supplied function.
 
302
          integration_error = true;
 
303
          break;
 
304
 
 
305
        default:
 
306
          integration_error = true;
 
307
          (*current_liboctave_error_handler)
 
308
            ("unrecognized value of istate (= %d) returned from lsode",
 
309
             istate);
 
310
          break;
 
311
        }
 
312
    }
 
313
 
 
314
  return retval;
 
315
}
 
316
 
 
317
std::string
 
318
LSODE::error_message (void) const
 
319
{
 
320
  std::string retval;
 
321
 
 
322
  std::ostringstream buf;
 
323
  buf << t;
 
324
  std::string t_curr = buf.str ();
 
325
 
 
326
  switch (istate)
 
327
    {
 
328
    case 1:
 
329
      retval = "prior to initial integration step";
 
330
      break;
 
331
 
 
332
    case 2:
 
333
      retval = "successful exit";
 
334
      break;
 
335
          
 
336
    case 3:
 
337
      retval = "prior to continuation call with modified parameters";
 
338
      break;
 
339
          
 
340
    case -1:
 
341
      retval = std::string ("excess work on this call (t = ")
 
342
        + t_curr + "; perhaps wrong integration method)";
 
343
      break;
 
344
 
 
345
    case -2:
 
346
      retval = "excess accuracy requested (tolerances too small)";
 
347
      break;
 
348
 
 
349
    case -3:
 
350
      retval = "invalid input detected (see printed message)";
 
351
      break;
 
352
 
 
353
    case -4:
 
354
      retval = std::string ("repeated error test failures (t = ")
 
355
        + t_curr + "check all inputs)";
 
356
      break;
 
357
 
 
358
    case -5:
 
359
      retval = std::string ("repeated convergence failures (t = ")
 
360
        + t_curr
 
361
        + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)";
 
362
      break;
 
363
 
 
364
    case -6:
 
365
      retval = std::string ("error weight became zero during problem. (t = ")
 
366
        + t_curr
 
367
        + "; solution component i vanished, and atol or atol(i) == 0)";
 
368
      break;
 
369
 
 
370
    case -13:
 
371
      retval = "return requested in user-supplied function (t = "
 
372
        + t_curr + ")";
 
373
      break;
 
374
 
 
375
    default:
 
376
      retval = "unknown error state";
 
377
      break;
 
378
    }
 
379
 
 
380
  return retval;
 
381
}
 
382
 
 
383
Matrix
 
384
LSODE::do_integrate (const ColumnVector& tout)
 
385
{
 
386
  Matrix retval;
 
387
 
 
388
  octave_idx_type n_out = tout.capacity ();
 
389
  octave_idx_type n = size ();
 
390
 
 
391
  if (n_out > 0 && n > 0)
 
392
    {
 
393
      retval.resize (n_out, n);
 
394
 
 
395
      for (octave_idx_type i = 0; i < n; i++)
 
396
        retval.elem (0, i) = x.elem (i);
 
397
 
 
398
      for (octave_idx_type j = 1; j < n_out; j++)
 
399
        {
 
400
          ColumnVector x_next = do_integrate (tout.elem (j));
 
401
 
 
402
          if (integration_error)
 
403
            return retval;
 
404
 
 
405
          for (octave_idx_type i = 0; i < n; i++)
 
406
            retval.elem (j, i) = x_next.elem (i);
 
407
        }
 
408
    }
 
409
 
 
410
  return retval;
 
411
}
 
412
 
 
413
Matrix
 
414
LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
 
415
{
 
416
  Matrix retval;
 
417
 
 
418
  octave_idx_type n_out = tout.capacity ();
 
419
  octave_idx_type n = size ();
 
420
 
 
421
  if (n_out > 0 && n > 0)
 
422
    {
 
423
      retval.resize (n_out, n);
 
424
 
 
425
      for (octave_idx_type i = 0; i < n; i++)
 
426
        retval.elem (0, i) = x.elem (i);
 
427
 
 
428
      octave_idx_type n_crit = tcrit.capacity ();
 
429
 
 
430
      if (n_crit > 0)
 
431
        {
 
432
          octave_idx_type i_crit = 0;
 
433
          octave_idx_type i_out = 1;
 
434
          double next_crit = tcrit.elem (0);
 
435
          double next_out;
 
436
          while (i_out < n_out)
 
437
            {
 
438
              bool do_restart = false;
 
439
 
 
440
              next_out = tout.elem (i_out);
 
441
              if (i_crit < n_crit)
 
442
                next_crit = tcrit.elem (i_crit);
 
443
 
 
444
              octave_idx_type save_output;
 
445
              double t_out;
 
446
 
 
447
              if (next_crit == next_out)
 
448
                {
 
449
                  set_stop_time (next_crit);
 
450
                  t_out = next_out;
 
451
                  save_output = 1;
 
452
                  i_out++;
 
453
                  i_crit++;
 
454
                  do_restart = true;
 
455
                }
 
456
              else if (next_crit < next_out)
 
457
                {
 
458
                  if (i_crit < n_crit)
 
459
                    {
 
460
                      set_stop_time (next_crit);
 
461
                      t_out = next_crit;
 
462
                      save_output = 0;
 
463
                      i_crit++;
 
464
                      do_restart = true;
 
465
                    }
 
466
                  else
 
467
                    {
 
468
                      clear_stop_time ();
 
469
                      t_out = next_out;
 
470
                      save_output = 1;
 
471
                      i_out++;
 
472
                    }
 
473
                }
 
474
              else
 
475
                {
 
476
                  set_stop_time (next_crit);
 
477
                  t_out = next_out;
 
478
                  save_output = 1;
 
479
                  i_out++;
 
480
                }
 
481
 
 
482
              ColumnVector x_next = do_integrate (t_out);
 
483
 
 
484
              if (integration_error)
 
485
                return retval;
 
486
 
 
487
              if (save_output)
 
488
                {
 
489
                  for (octave_idx_type i = 0; i < n; i++)
 
490
                    retval.elem (i_out-1, i) = x_next.elem (i);
 
491
                }
 
492
 
 
493
              if (do_restart)
 
494
                force_restart ();
 
495
            }
 
496
        }
 
497
      else
 
498
        {
 
499
          retval = do_integrate (tout);
 
500
 
 
501
          if (integration_error)
 
502
            return retval;
 
503
        }
 
504
    }
 
505
 
 
506
  return retval;
 
507
}
 
508
 
 
509
/*
 
510
;;; Local Variables: ***
 
511
;;; mode: C++ ***
 
512
;;; End: ***
 
513
*/