Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.7.0/ALGLIB-3.7.0/diffequations.cs @ 11920

Last change on this file since 11920 was 9268, checked in by mkommend, 12 years ago

#2007: Added AlgLib 3.7.0 to the external libraries.

File size: 34.6 KB
Line 
1/*************************************************************************
2Copyright (c) Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17>>> END OF LICENSE >>>
18*************************************************************************/
19#pragma warning disable 162
20#pragma warning disable 219
21using System;
22
23public partial class alglib
24{
25
26
27    /*************************************************************************
28
29    *************************************************************************/
30    public class odesolverstate
31    {
32        //
33        // Public declarations
34        //
35        public bool needdy { get { return _innerobj.needdy; } set { _innerobj.needdy = value; } }
36        public double[] y { get { return _innerobj.y; } }
37        public double[] dy { get { return _innerobj.dy; } }
38        public double x { get { return _innerobj.x; } set { _innerobj.x = value; } }
39
40        public odesolverstate()
41        {
42            _innerobj = new odesolver.odesolverstate();
43        }
44
45        //
46        // Although some of declarations below are public, you should not use them
47        // They are intended for internal use only
48        //
49        private odesolver.odesolverstate _innerobj;
50        public odesolver.odesolverstate innerobj { get { return _innerobj; } }
51        public odesolverstate(odesolver.odesolverstate obj)
52        {
53            _innerobj = obj;
54        }
55    }
56
57
58    /*************************************************************************
59
60    *************************************************************************/
61    public class odesolverreport
62    {
63        //
64        // Public declarations
65        //
66        public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } }
67        public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } }
68
69        public odesolverreport()
70        {
71            _innerobj = new odesolver.odesolverreport();
72        }
73
74        //
75        // Although some of declarations below are public, you should not use them
76        // They are intended for internal use only
77        //
78        private odesolver.odesolverreport _innerobj;
79        public odesolver.odesolverreport innerobj { get { return _innerobj; } }
80        public odesolverreport(odesolver.odesolverreport obj)
81        {
82            _innerobj = obj;
83        }
84    }
85
86    /*************************************************************************
87    Cash-Karp adaptive ODE solver.
88
89    This subroutine solves ODE  Y'=f(Y,x)  with  initial  conditions  Y(xs)=Ys
90    (here Y may be single variable or vector of N variables).
91
92    INPUT PARAMETERS:
93        Y       -   initial conditions, array[0..N-1].
94                    contains values of Y[] at X[0]
95        N       -   system size
96        X       -   points at which Y should be tabulated, array[0..M-1]
97                    integrations starts at X[0], ends at X[M-1],  intermediate
98                    values at X[i] are returned too.
99                    SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
100        M       -   number of intermediate points + first point + last point:
101                    * M>2 means that you need both Y(X[M-1]) and M-2 values at
102                      intermediate points
103                    * M=2 means that you want just to integrate from  X[0]  to
104                      X[1] and don't interested in intermediate values.
105                    * M=1 means that you don't want to integrate :)
106                      it is degenerate case, but it will be handled correctly.
107                    * M<1 means error
108        Eps     -   tolerance (absolute/relative error on each  step  will  be
109                    less than Eps). When passing:
110                    * Eps>0, it means desired ABSOLUTE error
111                    * Eps<0, it means desired RELATIVE error.  Relative errors
112                      are calculated with respect to maximum values of  Y seen
113                      so far. Be careful to use this criterion  when  starting
114                      from Y[] that are close to zero.
115        H       -   initial  step  lenth,  it  will  be adjusted automatically
116                    after the first  step.  If  H=0,  step  will  be  selected
117                    automatically  (usualy  it  will  be  equal  to  0.001  of
118                    min(x[i]-x[j])).
119
120    OUTPUT PARAMETERS
121        State   -   structure which stores algorithm state between  subsequent
122                    calls of OdeSolverIteration. Used for reverse communication.
123                    This structure should be passed  to the OdeSolverIteration
124                    subroutine.
125
126    SEE ALSO
127        AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
128
129
130      -- ALGLIB --
131         Copyright 01.09.2009 by Bochkanov Sergey
132    *************************************************************************/
133    public static void odesolverrkck(double[] y, int n, double[] x, int m, double eps, double h, out odesolverstate state)
134    {
135        state = new odesolverstate();
136        odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);
137        return;
138    }
139    public static void odesolverrkck(double[] y, double[] x, double eps, double h, out odesolverstate state)
140    {
141        int n;
142        int m;
143
144        state = new odesolverstate();
145        n = ap.len(y);
146        m = ap.len(x);
147        odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);
148
149        return;
150    }
151
152    /*************************************************************************
153    This function provides reverse communication interface
154    Reverse communication interface is not documented or recommended to use.
155    See below for functions which provide better documented API
156    *************************************************************************/
157    public static bool odesolveriteration(odesolverstate state)
158    {
159
160        bool result = odesolver.odesolveriteration(state.innerobj);
161        return result;
162    }
163    /*************************************************************************
164    This function is used to launcn iterations of ODE solver
165
166    It accepts following parameters:
167        diff    -   callback which calculates dy/dx for given y and x
168        obj     -   optional object which is passed to diff; can be NULL
169
170
171      -- ALGLIB --
172         Copyright 01.09.2009 by Bochkanov Sergey
173
174    *************************************************************************/
175    public static void odesolversolve(odesolverstate state, ndimensional_ode_rp diff, object obj)
176    {
177        if( diff==null )
178            throw new alglibexception("ALGLIB: error in 'odesolversolve()' (diff is null)");
179        while( alglib.odesolveriteration(state) )
180        {
181            if( state.needdy )
182            {
183                diff(state.innerobj.y, state.innerobj.x, state.innerobj.dy, obj);
184                continue;
185            }
186            throw new alglibexception("ALGLIB: unexpected error in 'odesolversolve'");
187        }
188    }
189
190
191
192    /*************************************************************************
193    ODE solver results
194
195    Called after OdeSolverIteration returned False.
196
197    INPUT PARAMETERS:
198        State   -   algorithm state (used by OdeSolverIteration).
199
200    OUTPUT PARAMETERS:
201        M       -   number of tabulated values, M>=1
202        XTbl    -   array[0..M-1], values of X
203        YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]
204        Rep     -   solver report:
205                    * Rep.TerminationType completetion code:
206                        * -2    X is not ordered  by  ascending/descending  or
207                                there are non-distinct X[],  i.e.  X[i]=X[i+1]
208                        * -1    incorrect parameters were specified
209                        *  1    task has been solved
210                    * Rep.NFEV contains number of function calculations
211
212      -- ALGLIB --
213         Copyright 01.09.2009 by Bochkanov Sergey
214    *************************************************************************/
215    public static void odesolverresults(odesolverstate state, out int m, out double[] xtbl, out double[,] ytbl, out odesolverreport rep)
216    {
217        m = 0;
218        xtbl = new double[0];
219        ytbl = new double[0,0];
220        rep = new odesolverreport();
221        odesolver.odesolverresults(state.innerobj, ref m, ref xtbl, ref ytbl, rep.innerobj);
222        return;
223    }
224
225}
226public partial class alglib
227{
228    public class odesolver
229    {
230        public class odesolverstate : apobject
231        {
232            public int n;
233            public int m;
234            public double xscale;
235            public double h;
236            public double eps;
237            public bool fraceps;
238            public double[] yc;
239            public double[] escale;
240            public double[] xg;
241            public int solvertype;
242            public bool needdy;
243            public double x;
244            public double[] y;
245            public double[] dy;
246            public double[,] ytbl;
247            public int repterminationtype;
248            public int repnfev;
249            public double[] yn;
250            public double[] yns;
251            public double[] rka;
252            public double[] rkc;
253            public double[] rkcs;
254            public double[,] rkb;
255            public double[,] rkk;
256            public rcommstate rstate;
257            public odesolverstate()
258            {
259                init();
260            }
261            public override void init()
262            {
263                yc = new double[0];
264                escale = new double[0];
265                xg = new double[0];
266                y = new double[0];
267                dy = new double[0];
268                ytbl = new double[0,0];
269                yn = new double[0];
270                yns = new double[0];
271                rka = new double[0];
272                rkc = new double[0];
273                rkcs = new double[0];
274                rkb = new double[0,0];
275                rkk = new double[0,0];
276                rstate = new rcommstate();
277            }
278            public override alglib.apobject make_copy()
279            {
280                odesolverstate _result = new odesolverstate();
281                _result.n = n;
282                _result.m = m;
283                _result.xscale = xscale;
284                _result.h = h;
285                _result.eps = eps;
286                _result.fraceps = fraceps;
287                _result.yc = (double[])yc.Clone();
288                _result.escale = (double[])escale.Clone();
289                _result.xg = (double[])xg.Clone();
290                _result.solvertype = solvertype;
291                _result.needdy = needdy;
292                _result.x = x;
293                _result.y = (double[])y.Clone();
294                _result.dy = (double[])dy.Clone();
295                _result.ytbl = (double[,])ytbl.Clone();
296                _result.repterminationtype = repterminationtype;
297                _result.repnfev = repnfev;
298                _result.yn = (double[])yn.Clone();
299                _result.yns = (double[])yns.Clone();
300                _result.rka = (double[])rka.Clone();
301                _result.rkc = (double[])rkc.Clone();
302                _result.rkcs = (double[])rkcs.Clone();
303                _result.rkb = (double[,])rkb.Clone();
304                _result.rkk = (double[,])rkk.Clone();
305                _result.rstate = (rcommstate)rstate.make_copy();
306                return _result;
307            }
308        };
309
310
311        public class odesolverreport : apobject
312        {
313            public int nfev;
314            public int terminationtype;
315            public odesolverreport()
316            {
317                init();
318            }
319            public override void init()
320            {
321            }
322            public override alglib.apobject make_copy()
323            {
324                odesolverreport _result = new odesolverreport();
325                _result.nfev = nfev;
326                _result.terminationtype = terminationtype;
327                return _result;
328            }
329        };
330
331
332
333
334        public const double odesolvermaxgrow = 3.0;
335        public const double odesolvermaxshrink = 10.0;
336
337
338        /*************************************************************************
339        Cash-Karp adaptive ODE solver.
340
341        This subroutine solves ODE  Y'=f(Y,x)  with  initial  conditions  Y(xs)=Ys
342        (here Y may be single variable or vector of N variables).
343
344        INPUT PARAMETERS:
345            Y       -   initial conditions, array[0..N-1].
346                        contains values of Y[] at X[0]
347            N       -   system size
348            X       -   points at which Y should be tabulated, array[0..M-1]
349                        integrations starts at X[0], ends at X[M-1],  intermediate
350                        values at X[i] are returned too.
351                        SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
352            M       -   number of intermediate points + first point + last point:
353                        * M>2 means that you need both Y(X[M-1]) and M-2 values at
354                          intermediate points
355                        * M=2 means that you want just to integrate from  X[0]  to
356                          X[1] and don't interested in intermediate values.
357                        * M=1 means that you don't want to integrate :)
358                          it is degenerate case, but it will be handled correctly.
359                        * M<1 means error
360            Eps     -   tolerance (absolute/relative error on each  step  will  be
361                        less than Eps). When passing:
362                        * Eps>0, it means desired ABSOLUTE error
363                        * Eps<0, it means desired RELATIVE error.  Relative errors
364                          are calculated with respect to maximum values of  Y seen
365                          so far. Be careful to use this criterion  when  starting
366                          from Y[] that are close to zero.
367            H       -   initial  step  lenth,  it  will  be adjusted automatically
368                        after the first  step.  If  H=0,  step  will  be  selected
369                        automatically  (usualy  it  will  be  equal  to  0.001  of
370                        min(x[i]-x[j])).
371
372        OUTPUT PARAMETERS
373            State   -   structure which stores algorithm state between  subsequent
374                        calls of OdeSolverIteration. Used for reverse communication.
375                        This structure should be passed  to the OdeSolverIteration
376                        subroutine.
377
378        SEE ALSO
379            AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
380
381
382          -- ALGLIB --
383             Copyright 01.09.2009 by Bochkanov Sergey
384        *************************************************************************/
385        public static void odesolverrkck(double[] y,
386            int n,
387            double[] x,
388            int m,
389            double eps,
390            double h,
391            odesolverstate state)
392        {
393            alglib.ap.assert(n>=1, "ODESolverRKCK: N<1!");
394            alglib.ap.assert(m>=1, "ODESolverRKCK: M<1!");
395            alglib.ap.assert(alglib.ap.len(y)>=n, "ODESolverRKCK: Length(Y)<N!");
396            alglib.ap.assert(alglib.ap.len(x)>=m, "ODESolverRKCK: Length(X)<M!");
397            alglib.ap.assert(apserv.isfinitevector(y, n), "ODESolverRKCK: Y contains infinite or NaN values!");
398            alglib.ap.assert(apserv.isfinitevector(x, m), "ODESolverRKCK: Y contains infinite or NaN values!");
399            alglib.ap.assert(math.isfinite(eps), "ODESolverRKCK: Eps is not finite!");
400            alglib.ap.assert((double)(eps)!=(double)(0), "ODESolverRKCK: Eps is zero!");
401            alglib.ap.assert(math.isfinite(h), "ODESolverRKCK: H is not finite!");
402            odesolverinit(0, y, n, x, m, eps, h, state);
403        }
404
405
406        /*************************************************************************
407
408          -- ALGLIB --
409             Copyright 01.09.2009 by Bochkanov Sergey
410        *************************************************************************/
411        public static bool odesolveriteration(odesolverstate state)
412        {
413            bool result = new bool();
414            int n = 0;
415            int m = 0;
416            int i = 0;
417            int j = 0;
418            int k = 0;
419            double xc = 0;
420            double v = 0;
421            double h = 0;
422            double h2 = 0;
423            bool gridpoint = new bool();
424            double err = 0;
425            double maxgrowpow = 0;
426            int klimit = 0;
427            int i_ = 0;
428
429           
430            //
431            // Reverse communication preparations
432            // I know it looks ugly, but it works the same way
433            // anywhere from C++ to Python.
434            //
435            // This code initializes locals by:
436            // * random values determined during code
437            //   generation - on first subroutine call
438            // * values from previous call - on subsequent calls
439            //
440            if( state.rstate.stage>=0 )
441            {
442                n = state.rstate.ia[0];
443                m = state.rstate.ia[1];
444                i = state.rstate.ia[2];
445                j = state.rstate.ia[3];
446                k = state.rstate.ia[4];
447                klimit = state.rstate.ia[5];
448                gridpoint = state.rstate.ba[0];
449                xc = state.rstate.ra[0];
450                v = state.rstate.ra[1];
451                h = state.rstate.ra[2];
452                h2 = state.rstate.ra[3];
453                err = state.rstate.ra[4];
454                maxgrowpow = state.rstate.ra[5];
455            }
456            else
457            {
458                n = -983;
459                m = -989;
460                i = -834;
461                j = 900;
462                k = -287;
463                klimit = 364;
464                gridpoint = false;
465                xc = -338;
466                v = -686;
467                h = 912;
468                h2 = 585;
469                err = 497;
470                maxgrowpow = -271;
471            }
472            if( state.rstate.stage==0 )
473            {
474                goto lbl_0;
475            }
476           
477            //
478            // Routine body
479            //
480           
481            //
482            // prepare
483            //
484            if( state.repterminationtype!=0 )
485            {
486                result = false;
487                return result;
488            }
489            n = state.n;
490            m = state.m;
491            h = state.h;
492            maxgrowpow = Math.Pow(odesolvermaxgrow, 5);
493            state.repnfev = 0;
494           
495            //
496            // some preliminary checks for internal errors
497            // after this we assume that H>0 and M>1
498            //
499            alglib.ap.assert((double)(state.h)>(double)(0), "ODESolver: internal error");
500            alglib.ap.assert(m>1, "ODESolverIteration: internal error");
501           
502            //
503            // choose solver
504            //
505            if( state.solvertype!=0 )
506            {
507                goto lbl_1;
508            }
509           
510            //
511            // Cask-Karp solver
512            // Prepare coefficients table.
513            // Check it for errors
514            //
515            state.rka = new double[6];
516            state.rka[0] = 0;
517            state.rka[1] = (double)1/(double)5;
518            state.rka[2] = (double)3/(double)10;
519            state.rka[3] = (double)3/(double)5;
520            state.rka[4] = 1;
521            state.rka[5] = (double)7/(double)8;
522            state.rkb = new double[6, 5];
523            state.rkb[1,0] = (double)1/(double)5;
524            state.rkb[2,0] = (double)3/(double)40;
525            state.rkb[2,1] = (double)9/(double)40;
526            state.rkb[3,0] = (double)3/(double)10;
527            state.rkb[3,1] = -((double)9/(double)10);
528            state.rkb[3,2] = (double)6/(double)5;
529            state.rkb[4,0] = -((double)11/(double)54);
530            state.rkb[4,1] = (double)5/(double)2;
531            state.rkb[4,2] = -((double)70/(double)27);
532            state.rkb[4,3] = (double)35/(double)27;
533            state.rkb[5,0] = (double)1631/(double)55296;
534            state.rkb[5,1] = (double)175/(double)512;
535            state.rkb[5,2] = (double)575/(double)13824;
536            state.rkb[5,3] = (double)44275/(double)110592;
537            state.rkb[5,4] = (double)253/(double)4096;
538            state.rkc = new double[6];
539            state.rkc[0] = (double)37/(double)378;
540            state.rkc[1] = 0;
541            state.rkc[2] = (double)250/(double)621;
542            state.rkc[3] = (double)125/(double)594;
543            state.rkc[4] = 0;
544            state.rkc[5] = (double)512/(double)1771;
545            state.rkcs = new double[6];
546            state.rkcs[0] = (double)2825/(double)27648;
547            state.rkcs[1] = 0;
548            state.rkcs[2] = (double)18575/(double)48384;
549            state.rkcs[3] = (double)13525/(double)55296;
550            state.rkcs[4] = (double)277/(double)14336;
551            state.rkcs[5] = (double)1/(double)4;
552            state.rkk = new double[6, n];
553           
554            //
555            // Main cycle consists of two iterations:
556            // * outer where we travel from X[i-1] to X[i]
557            // * inner where we travel inside [X[i-1],X[i]]
558            //
559            state.ytbl = new double[m, n];
560            state.escale = new double[n];
561            state.yn = new double[n];
562            state.yns = new double[n];
563            xc = state.xg[0];
564            for(i_=0; i_<=n-1;i_++)
565            {
566                state.ytbl[0,i_] = state.yc[i_];
567            }
568            for(j=0; j<=n-1; j++)
569            {
570                state.escale[j] = 0;
571            }
572            i = 1;
573        lbl_3:
574            if( i>m-1 )
575            {
576                goto lbl_5;
577            }
578           
579            //
580            // begin inner iteration
581            //
582        lbl_6:
583            if( false )
584            {
585                goto lbl_7;
586            }
587           
588            //
589            // truncate step if needed (beyond right boundary).
590            // determine should we store X or not
591            //
592            if( (double)(xc+h)>=(double)(state.xg[i]) )
593            {
594                h = state.xg[i]-xc;
595                gridpoint = true;
596            }
597            else
598            {
599                gridpoint = false;
600            }
601           
602            //
603            // Update error scale maximums
604            //
605            // These maximums are initialized by zeros,
606            // then updated every iterations.
607            //
608            for(j=0; j<=n-1; j++)
609            {
610                state.escale[j] = Math.Max(state.escale[j], Math.Abs(state.yc[j]));
611            }
612           
613            //
614            // make one step:
615            // 1. calculate all info needed to do step
616            // 2. update errors scale maximums using values/derivatives
617            //    obtained during (1)
618            //
619            // Take into account that we use scaling of X to reduce task
620            // to the form where x[0] < x[1] < ... < x[n-1]. So X is
621            // replaced by x=xscale*t, and dy/dx=f(y,x) is replaced
622            // by dy/dt=xscale*f(y,xscale*t).
623            //
624            for(i_=0; i_<=n-1;i_++)
625            {
626                state.yn[i_] = state.yc[i_];
627            }
628            for(i_=0; i_<=n-1;i_++)
629            {
630                state.yns[i_] = state.yc[i_];
631            }
632            k = 0;
633        lbl_8:
634            if( k>5 )
635            {
636                goto lbl_10;
637            }
638           
639            //
640            // prepare data for the next update of YN/YNS
641            //
642            state.x = state.xscale*(xc+state.rka[k]*h);
643            for(i_=0; i_<=n-1;i_++)
644            {
645                state.y[i_] = state.yc[i_];
646            }
647            for(j=0; j<=k-1; j++)
648            {
649                v = state.rkb[k,j];
650                for(i_=0; i_<=n-1;i_++)
651                {
652                    state.y[i_] = state.y[i_] + v*state.rkk[j,i_];
653                }
654            }
655            state.needdy = true;
656            state.rstate.stage = 0;
657            goto lbl_rcomm;
658        lbl_0:
659            state.needdy = false;
660            state.repnfev = state.repnfev+1;
661            v = h*state.xscale;
662            for(i_=0; i_<=n-1;i_++)
663            {
664                state.rkk[k,i_] = v*state.dy[i_];
665            }
666           
667            //
668            // update YN/YNS
669            //
670            v = state.rkc[k];
671            for(i_=0; i_<=n-1;i_++)
672            {
673                state.yn[i_] = state.yn[i_] + v*state.rkk[k,i_];
674            }
675            v = state.rkcs[k];
676            for(i_=0; i_<=n-1;i_++)
677            {
678                state.yns[i_] = state.yns[i_] + v*state.rkk[k,i_];
679            }
680            k = k+1;
681            goto lbl_8;
682        lbl_10:
683           
684            //
685            // estimate error
686            //
687            err = 0;
688            for(j=0; j<=n-1; j++)
689            {
690                if( !state.fraceps )
691                {
692                   
693                    //
694                    // absolute error is estimated
695                    //
696                    err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j]));
697                }
698                else
699                {
700                   
701                    //
702                    // Relative error is estimated
703                    //
704                    v = state.escale[j];
705                    if( (double)(v)==(double)(0) )
706                    {
707                        v = 1;
708                    }
709                    err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j])/v);
710                }
711            }
712           
713            //
714            // calculate new step, restart if necessary
715            //
716            if( (double)(maxgrowpow*err)<=(double)(state.eps) )
717            {
718                h2 = odesolvermaxgrow*h;
719            }
720            else
721            {
722                h2 = h*Math.Pow(state.eps/err, 0.2);
723            }
724            if( (double)(h2)<(double)(h/odesolvermaxshrink) )
725            {
726                h2 = h/odesolvermaxshrink;
727            }
728            if( (double)(err)>(double)(state.eps) )
729            {
730                h = h2;
731                goto lbl_6;
732            }
733           
734            //
735            // advance position
736            //
737            xc = xc+h;
738            for(i_=0; i_<=n-1;i_++)
739            {
740                state.yc[i_] = state.yn[i_];
741            }
742           
743            //
744            // update H
745            //
746            h = h2;
747           
748            //
749            // break on grid point
750            //
751            if( gridpoint )
752            {
753                goto lbl_7;
754            }
755            goto lbl_6;
756        lbl_7:
757           
758            //
759            // save result
760            //
761            for(i_=0; i_<=n-1;i_++)
762            {
763                state.ytbl[i,i_] = state.yc[i_];
764            }
765            i = i+1;
766            goto lbl_3;
767        lbl_5:
768            state.repterminationtype = 1;
769            result = false;
770            return result;
771        lbl_1:
772            result = false;
773            return result;
774           
775            //
776            // Saving state
777            //
778        lbl_rcomm:
779            result = true;
780            state.rstate.ia[0] = n;
781            state.rstate.ia[1] = m;
782            state.rstate.ia[2] = i;
783            state.rstate.ia[3] = j;
784            state.rstate.ia[4] = k;
785            state.rstate.ia[5] = klimit;
786            state.rstate.ba[0] = gridpoint;
787            state.rstate.ra[0] = xc;
788            state.rstate.ra[1] = v;
789            state.rstate.ra[2] = h;
790            state.rstate.ra[3] = h2;
791            state.rstate.ra[4] = err;
792            state.rstate.ra[5] = maxgrowpow;
793            return result;
794        }
795
796
797        /*************************************************************************
798        ODE solver results
799
800        Called after OdeSolverIteration returned False.
801
802        INPUT PARAMETERS:
803            State   -   algorithm state (used by OdeSolverIteration).
804
805        OUTPUT PARAMETERS:
806            M       -   number of tabulated values, M>=1
807            XTbl    -   array[0..M-1], values of X
808            YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]
809            Rep     -   solver report:
810                        * Rep.TerminationType completetion code:
811                            * -2    X is not ordered  by  ascending/descending  or
812                                    there are non-distinct X[],  i.e.  X[i]=X[i+1]
813                            * -1    incorrect parameters were specified
814                            *  1    task has been solved
815                        * Rep.NFEV contains number of function calculations
816
817          -- ALGLIB --
818             Copyright 01.09.2009 by Bochkanov Sergey
819        *************************************************************************/
820        public static void odesolverresults(odesolverstate state,
821            ref int m,
822            ref double[] xtbl,
823            ref double[,] ytbl,
824            odesolverreport rep)
825        {
826            double v = 0;
827            int i = 0;
828            int i_ = 0;
829
830            m = 0;
831            xtbl = new double[0];
832            ytbl = new double[0,0];
833
834            rep.terminationtype = state.repterminationtype;
835            if( rep.terminationtype>0 )
836            {
837                m = state.m;
838                rep.nfev = state.repnfev;
839                xtbl = new double[state.m];
840                v = state.xscale;
841                for(i_=0; i_<=state.m-1;i_++)
842                {
843                    xtbl[i_] = v*state.xg[i_];
844                }
845                ytbl = new double[state.m, state.n];
846                for(i=0; i<=state.m-1; i++)
847                {
848                    for(i_=0; i_<=state.n-1;i_++)
849                    {
850                        ytbl[i,i_] = state.ytbl[i,i_];
851                    }
852                }
853            }
854            else
855            {
856                rep.nfev = 0;
857            }
858        }
859
860
861        /*************************************************************************
862        Internal initialization subroutine
863        *************************************************************************/
864        private static void odesolverinit(int solvertype,
865            double[] y,
866            int n,
867            double[] x,
868            int m,
869            double eps,
870            double h,
871            odesolverstate state)
872        {
873            int i = 0;
874            double v = 0;
875            int i_ = 0;
876
877           
878            //
879            // Prepare RComm
880            //
881            state.rstate.ia = new int[5+1];
882            state.rstate.ba = new bool[0+1];
883            state.rstate.ra = new double[5+1];
884            state.rstate.stage = -1;
885            state.needdy = false;
886           
887            //
888            // check parameters.
889            //
890            if( (n<=0 || m<1) || (double)(eps)==(double)(0) )
891            {
892                state.repterminationtype = -1;
893                return;
894            }
895            if( (double)(h)<(double)(0) )
896            {
897                h = -h;
898            }
899           
900            //
901            // quick exit if necessary.
902            // after this block we assume that M>1
903            //
904            if( m==1 )
905            {
906                state.repnfev = 0;
907                state.repterminationtype = 1;
908                state.ytbl = new double[1, n];
909                for(i_=0; i_<=n-1;i_++)
910                {
911                    state.ytbl[0,i_] = y[i_];
912                }
913                state.xg = new double[m];
914                for(i_=0; i_<=m-1;i_++)
915                {
916                    state.xg[i_] = x[i_];
917                }
918                return;
919            }
920           
921            //
922            // check again: correct order of X[]
923            //
924            if( (double)(x[1])==(double)(x[0]) )
925            {
926                state.repterminationtype = -2;
927                return;
928            }
929            for(i=1; i<=m-1; i++)
930            {
931                if( ((double)(x[1])>(double)(x[0]) && (double)(x[i])<=(double)(x[i-1])) || ((double)(x[1])<(double)(x[0]) && (double)(x[i])>=(double)(x[i-1])) )
932                {
933                    state.repterminationtype = -2;
934                    return;
935                }
936            }
937           
938            //
939            // auto-select H if necessary
940            //
941            if( (double)(h)==(double)(0) )
942            {
943                v = Math.Abs(x[1]-x[0]);
944                for(i=2; i<=m-1; i++)
945                {
946                    v = Math.Min(v, Math.Abs(x[i]-x[i-1]));
947                }
948                h = 0.001*v;
949            }
950           
951            //
952            // store parameters
953            //
954            state.n = n;
955            state.m = m;
956            state.h = h;
957            state.eps = Math.Abs(eps);
958            state.fraceps = (double)(eps)<(double)(0);
959            state.xg = new double[m];
960            for(i_=0; i_<=m-1;i_++)
961            {
962                state.xg[i_] = x[i_];
963            }
964            if( (double)(x[1])>(double)(x[0]) )
965            {
966                state.xscale = 1;
967            }
968            else
969            {
970                state.xscale = -1;
971                for(i_=0; i_<=m-1;i_++)
972                {
973                    state.xg[i_] = -1*state.xg[i_];
974                }
975            }
976            state.yc = new double[n];
977            for(i_=0; i_<=n-1;i_++)
978            {
979                state.yc[i_] = y[i_];
980            }
981            state.solvertype = solvertype;
982            state.repterminationtype = 0;
983           
984            //
985            // Allocate arrays
986            //
987            state.y = new double[n];
988            state.dy = new double[n];
989        }
990
991
992    }
993}
994
Note: See TracBrowser for help on using the repository browser.