integrator_bdf.cpp
Go to the documentation of this file.
00001 /*
00002  *    This file is part of ACADO Toolkit.
00003  *
00004  *    ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
00005  *    Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
00006  *    Milan Vukov, Rien Quirynen, KU Leuven.
00007  *    Developed within the Optimization in Engineering Center (OPTEC)
00008  *    under supervision of Moritz Diehl. All rights reserved.
00009  *
00010  *    ACADO Toolkit is free software; you can redistribute it and/or
00011  *    modify it under the terms of the GNU Lesser General Public
00012  *    License as published by the Free Software Foundation; either
00013  *    version 3 of the License, or (at your option) any later version.
00014  *
00015  *    ACADO Toolkit is distributed in the hope that it will be useful,
00016  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  *    Lesser General Public License for more details.
00019  *
00020  *    You should have received a copy of the GNU Lesser General Public
00021  *    License along with ACADO Toolkit; if not, write to the Free Software
00022  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00023  *
00024  */
00025 
00026 
00027 
00034 #include <acado/utils/acado_utils.hpp>
00035 #include <acado/matrix_vector/matrix_vector.hpp>
00036 #include <acado/symbolic_expression/symbolic_expression.hpp>
00037 #include <acado/function/function_.hpp>
00038 #include <acado/integrator/integrator.hpp>
00039 
00040 using namespace std;
00041 
00042 BEGIN_NAMESPACE_ACADO
00043 
00044 
00045 //
00046 // PUBLIC MEMBER FUNCTIONS:
00047 //
00048 
00049 IntegratorBDF::IntegratorBDF( )
00050               :Integrator( ){
00051 
00052    initializeVariables();
00053 }
00054 
00055 
00056 IntegratorBDF::IntegratorBDF( const DifferentialEquation& rhs_ )
00057               :Integrator( ){
00058 
00059     init( rhs_ );
00060 }
00061 
00062 IntegratorBDF::IntegratorBDF( const IntegratorBDF& arg )
00063               :Integrator( arg ){
00064 
00065     constructAll( arg );
00066 }
00067 
00068 
00069 
00070 IntegratorBDF::~IntegratorBDF( ){
00071 
00072     deleteAll();
00073 }
00074 
00075 
00076 IntegratorBDF& IntegratorBDF::operator=( const IntegratorBDF& arg ){
00077 
00078     if ( this != &arg ){
00079         deleteAll();
00080         Integrator::operator=( arg );
00081         constructAll(arg);
00082     }
00083     return *this;
00084 }
00085 
00086 
00087 
00088 returnValue IntegratorBDF::init( const DifferentialEquation &rhs_ ){
00089 
00090     // RHS:
00091     // ---------
00092     rhs = new DifferentialEquation( rhs_ );
00093     m   = rhs->getDim ();
00094     ma  = rhs->getNXA ();
00095     mdx = rhs->getNDX ();
00096     mn  = rhs->getN   ();
00097     mu  = rhs->getNU  ();
00098     mui = rhs->getNUI ();
00099     mp  = rhs->getNP  ();
00100     mpi = rhs->getNPI ();
00101     mw  = rhs->getNW  ();
00102     md  = rhs->getNumDynamicEquations();
00103 
00104     rhs->makeImplicit();
00105     allocateMemory();
00106 
00107     return SUCCESSFUL_RETURN;
00108 }
00109 
00110 
00111 void IntegratorBDF::allocateMemory( ){
00112 
00113 
00114     int run1, run2, run3;
00115 
00116     if( m < 1 ){
00117         ACADOERROR(RET_TRIVIAL_RHS);
00118         return;
00119     }
00120 
00121     if( mdx > m-ma ){
00122         ACADOERROR(RET_TO_MANY_DIFFERENTIAL_STATE_DERIVATIVES);
00123         return;
00124     }
00125 
00126     if( m < md+ma ) md = m - ma;
00127 
00128 
00129     initializeVariables();
00130 
00131     initialAlgebraicResiduum = new double[m];
00132     relaxationConstant       = 5.0;
00133 
00134 
00135     // BUTCHER TABLEAU:
00136     // ----------------
00137     dim = 7;
00138     A  = new double*[dim];
00139     b4 = new double [dim];
00140     b5 = new double [dim];
00141     c  = new double [dim];
00142 
00143     for( run1 = 0; run1 < dim; run1++ ){
00144         A[run1] = new double[dim];
00145     }
00146 
00147     initializeButcherTableau();
00148 
00149     // RK-STARTER:
00150     // -----------
00151     ndir = rhs->getNumberOfVariables() + 1 + 2*md;
00152 
00153     eta4  = new double [m];
00154     eta5  = new double [m];
00155     eta4_ = new double [m];
00156     eta5_ = new double [m];
00157 
00158     for( run1 = 0; run1 < m; run1++ ){
00159 
00160         eta4 [run1] = 0.0;
00161         eta5 [run1] = 0.0;
00162         eta4_[run1] = 0.0;
00163         eta5_[run1] = 0.0;
00164     }
00165 
00166 
00167     k     = new double**[4];
00168     k2    = new double**[4];
00169     l     = new double**[4];
00170     l2    = new double**[4];
00171 
00172     for( run3 = 0; run3 < 4; run3++ ){
00173 
00174         k   [run3]  = new double*[dim];
00175         k2  [run3]  = new double*[dim];
00176         l   [run3]  = new double*[dim];
00177         l2  [run3]  = new double*[dim];
00178 
00179         for( run1 = 0; run1 < dim; run1++ ){
00180             k    [run3][run1] = new double[m];
00181             k2   [run3][run1] = new double[m];
00182 
00183             for( run2 = 0; run2 < m; run2++ ){
00184                 k [run3][run1][run2] = 0.0;
00185                 k2[run3][run1][run2] = 0.0;
00186             }
00187 
00188             l    [run3][run1] = new double[ndir];
00189             l2   [run3][run1] = new double[ndir];
00190 
00191             for( run2 = 0; run2 < ndir; run2++ ){
00192                 l [run3][run1][run2] = 0.0;
00193                 l2[run3][run1][run2] = 0.0;
00194             }
00195         }
00196     }
00197 
00198 
00199     iseed = new double[ndir];
00200 
00201     x     = new double [ndir];
00202 
00203     for( run1 = 0; run1 < ndir; run1++ ){
00204         x[run1]     = 0.0;
00205         iseed[run1] = 0.0;
00206     }
00207 
00208     t = 0.0;
00209 
00210 
00211     // BDF-METHOD:
00212     // -----------
00213     nstep     = 5;
00214     psi       = (double**)calloc(1,sizeof(double*));
00215     psi_      = new double[nstep];
00216     gamma     = (double**)calloc(1,sizeof(double*));
00217 
00218     c2.init(m);
00219     c2.setZero();
00220 
00221     nOfNewtonSteps = (int*)calloc(1,sizeof(int));
00222     eta            = new double*[4];
00223     eta2           = new double*[4];
00224 
00225     nOfNewtonSteps[0] = 0;
00226 
00227     psi   [0] = (double*)calloc(nstep,sizeof(double));
00228     gamma [0] = (double*)calloc(nstep,sizeof(double));
00229 
00230     nablaY.init( nstep, m );
00231     nablaY.setZero();
00232 
00233     nablaY_.init( nstep, m );
00234     nablaY_.setZero();
00235 
00236     phi.init( nstep, m );
00237     phi.setZero();
00238 
00239     delta.init( nstep, m );
00240     delta.setZero();
00241 
00242 
00243     for( run1 = 0; run1 < nstep; run1++ ){
00244          psi_    [run1] = 0.0;
00245          psi  [0][run1] = 0.0;
00246          gamma[0][run1] = 0.0;
00247     }
00248 
00249     maxNM = 1;
00250     M     = (DMatrix**)calloc(maxNM,sizeof(DMatrix*));
00251 //    qr.resize( maxNM );
00252     M_index   = (int*)calloc(maxNM,sizeof(int));
00253 
00254     M_index[0] = 0;
00255     M      [0] = 0;
00256     nOfM       = 0;
00257 
00258 
00259     for( run1 = 0; run1 < 4; run1++ ){
00260         eta [run1] = new double[m];
00261         eta2[run1] = new double[m];
00262         for( run2 = 0; run2 < m; run2++ ){
00263             eta [run1][run2] = 0.0;
00264             eta2[run1][run2] = 0.0;
00265         }
00266     }
00267 
00268 
00269     F  = new double[m];
00270     F2 = new double[m];
00271 
00272 
00273     // INTERNAL INDEX LISTS:
00274     // ---------------------
00275     diff_index = new int[m];
00276 
00277     for( run1 = 0; run1 < md; run1++ ){
00278         diff_index[run1] = rhs->getStateEnumerationIndex( run1 );
00279         if( diff_index[run1] == rhs->getNumberOfVariables() ){
00280             diff_index[run1] = diff_index[run1] + 1 + run1;
00281         }
00282     }
00283 
00284     ddiff_index = new int[md];
00285 
00286     for( run1 = 0; run1 < md; run1++ ){
00287         ddiff_index[run1] = rhs->index( VT_DDIFFERENTIAL_STATE, run1 );
00288         if( ddiff_index[run1] == rhs->getNumberOfVariables() ){
00289             ddiff_index[run1] = ddiff_index[run1] + 1 + md + run1;
00290         }
00291     }
00292 
00293 
00294     alg_index = new int[ma];
00295 
00296     for( run1 = 0; run1 < ma; run1++ ){
00297         alg_index [run1]    = rhs->index( VT_ALGEBRAIC_STATE, run1 );
00298         diff_index[md+run1] = alg_index [run1];
00299     }
00300 
00301     control_index       = new int[mu ];
00302 
00303     for( run1 = 0; run1 < mu; run1++ ){
00304         control_index[run1] = rhs->index( VT_CONTROL, run1 );
00305     }
00306 
00307     parameter_index     = new int[mp ];
00308 
00309     for( run1 = 0; run1 < mp; run1++ ){
00310         parameter_index[run1] = rhs->index( VT_PARAMETER, run1 );
00311     }
00312 
00313     int_control_index   = new int[mui];
00314 
00315     for( run1 = 0; run1 < mui; run1++ ){
00316         int_control_index[run1] = rhs->index( VT_INTEGER_CONTROL, run1 );
00317     }
00318 
00319     int_parameter_index = new int[mpi];
00320 
00321     for( run1 = 0; run1 < mpi; run1++ ){
00322         int_parameter_index[run1] = rhs->index( VT_INTEGER_PARAMETER, run1 );
00323     }
00324 
00325     disturbance_index   = new int[mw ];
00326 
00327     for( run1 = 0; run1 < mw; run1++ ){
00328         disturbance_index[run1] = rhs->index( VT_DISTURBANCE, run1 );
00329     }
00330 
00331     time_index = rhs->index( VT_TIME, 0 );
00332 
00333 
00334     // OTHERS:
00335     // -------
00336     diff_scale.init(m);
00337 
00338     for( run1 = 0; run1 < md; run1++ ){
00339         diff_scale(run1) = rhs->scale( VT_DIFFERENTIAL_STATE, run1 );
00340     }
00341     for( run1 = 0; run1 < ma; run1++ ){
00342         diff_scale(md+run1) = rhs->scale( VT_ALGEBRAIC_STATE, run1 );
00343     }
00344 
00345     initial_guess    = new double[m];
00346 
00347     for( run1 = 0; run1 < m; run1++ ){
00348 
00349         initial_guess[run1] = 0.0;
00350     }
00351 
00352 
00353     // SENSITIVITIES:
00354     // --------------
00355     c2G .init(md);
00356     c2G2.init(md);
00357     c2G3.init(md);
00358 
00359     phiG.init(nstep,m);
00360     phiG2.init(nstep,m);
00361     phiG3.init(nstep,m);
00362 
00363     deltaG.init(nstep,m);
00364     deltaG2.init(nstep,m);
00365     deltaG3.init(nstep,m);
00366 
00367     kH         = new double**[dim];
00368     for( run1 = 0; run1 < dim; run1++ ){
00369          kH[run1] = new double*[nstep];
00370          for( run2 = 0; run2 < nstep; run2++ ){
00371              kH[run1][run2] = new double[ndir];
00372          }
00373     }
00374     zH         = new double*[6];
00375     for( run1 = 0; run1 < 6; run1++ ){
00376          zH[run1] = new double[ndir];
00377     }
00378 
00379     kH2        = new double**[dim];
00380     for( run1 = 0; run1 < dim; run1++ ){
00381          kH2[run1] = new double*[nstep];
00382          for( run2 = 0; run2 < nstep; run2++ ){
00383              kH2[run1][run2] = new double[ndir];
00384          }
00385     }
00386     zH2        = new double*[6];
00387     for( run1 = 0; run1 < 6; run1++ ){
00388          zH2[run1] = new double[ndir];
00389     }
00390     kH3        = new double**[dim];
00391     for( run1 = 0; run1 < dim; run1++ ){
00392          kH3[run1] = new double*[nstep];
00393          for( run2 = 0; run2 < nstep; run2++ ){
00394              kH3[run1][run2] = new double[ndir];
00395          }
00396     }
00397     zH3        = new double*[6];
00398     for( run1 = 0; run1 < 6; run1++ ){
00399          zH3[run1] = new double[ndir];
00400     }
00401 
00402 
00403     // STORAGE:
00404     // --------
00405     maxAlloc = 1;
00406 }
00407 
00408 
00409 void IntegratorBDF::initializeVariables(){
00410 
00411     initialAlgebraicResiduum = 0  ;
00412     relaxationConstant       = 0.0;
00413 
00414     dim = 0; A = 0; b4 = 0; b5 = 0; c = 0;
00415 
00416     eta4  = 0; eta5  = 0;
00417     eta4_ = 0; eta5_ = 0;
00418 
00419     k = 0; k2 = 0; l = 0; l2 = 0;
00420 
00421     iseed = 0; nstep = 0; psi   = 0;
00422     psi_  = 0; gamma = 0; eta   = 0;
00423     eta2  = 0; x     = 0; t     = 0;
00424 
00425     nOfNewtonSteps = 0;
00426     maxNM = 0; M = 0; M_index = 0; nOfM = 0;
00427 
00428     F  = 0; F2 = 0;
00429 
00430     initial_guess = 0;
00431 
00432     ndir = 0; G = 0; etaG = 0;
00433 
00434     G2  = 0; G3   = 0; etaG2 = 0; etaG3 = 0;
00435     H   = 0; etaH = 0; kH    = 0; zH    = 0;
00436     H2  = 0; H3   = 0; etaH2 = 0; etaH3 = 0;
00437     kH2 = 0; zH2  = 0; kH3   = 0; zH3   = 0;
00438 
00439     maxAlloc = 0;
00440 
00441     nFcnEvaluations = 0;
00442     nJacEvaluations = 0;
00443 }
00444 
00445 
00446 void IntegratorBDF::constructAll( const IntegratorBDF& arg ){
00447 
00448     int run1, run2, run3;
00449 
00450     rhs = new DifferentialEquation( *arg.rhs );
00451 
00452     m   = arg.m;
00453     ma  = arg.ma;
00454     mdx = arg.mdx;
00455     mn  = arg.mn;
00456     mu  = arg.mu;
00457     mui = arg.mui;
00458     mp  = arg.mp;
00459     mpi = arg.mpi;
00460     mw  = arg.mw;
00461     md  = m-ma;
00462 
00463     initialAlgebraicResiduum = new double[m]         ;
00464     relaxationConstant       = arg.relaxationConstant;
00465 
00466 
00467     // BUTCHER TABLEAU:
00468     // ----------------
00469     dim = arg.dim;
00470     A  = new double*[dim];
00471     b4 = new double [dim];
00472     b5 = new double [dim];
00473     c  = new double [dim];
00474 
00475     for( run1 = 0; run1 < dim; run1++ ){
00476         A[run1] = new double[dim];
00477     }
00478 
00479     initializeButcherTableau();
00480 
00481 
00482     // RK-STARTER:
00483     // -----------
00484     ndir = rhs->getNumberOfVariables() + 1 + 2*md;
00485 
00486     eta4  = new double [m];
00487     eta5  = new double [m];
00488     eta4_ = new double [m];
00489     eta5_ = new double [m];
00490 
00491     for( run1 = 0; run1 < m; run1++ ){
00492 
00493         eta4 [run1] = 0.0;
00494         eta5 [run1] = 0.0;
00495         eta4_[run1] = 0.0;
00496         eta5_[run1] = 0.0;
00497     }
00498 
00499 
00500     k     = new double**[4];
00501     k2    = new double**[4];
00502     l     = new double**[4];
00503     l2    = new double**[4];
00504 
00505     for( run3 = 0; run3 < 4; run3++ ){
00506 
00507         k   [run3]  = new double*[dim];
00508         k2  [run3]  = new double*[dim];
00509         l   [run3]  = new double*[dim];
00510         l2  [run3]  = new double*[dim];
00511 
00512         for( run1 = 0; run1 < dim; run1++ ){
00513             k    [run3][run1] = new double[m];
00514             k2   [run3][run1] = new double[m];
00515 
00516             for( run2 = 0; run2 < m; run2++ ){
00517                 k [run3][run1][run2] = 0.0;
00518                 k2[run3][run1][run2] = 0.0;
00519             }
00520 
00521             l    [run3][run1] = new double[ndir];
00522             l2   [run3][run1] = new double[ndir];
00523 
00524             for( run2 = 0; run2 < ndir; run2++ ){
00525                 l [run3][run1][run2] = 0.0;
00526                 l2[run3][run1][run2] = 0.0;
00527             }
00528         }
00529     }
00530 
00531     iseed = new double[ndir];
00532 
00533     x     = new double [ndir];
00534 
00535     for( run1 = 0; run1 < ndir; run1++ ){
00536         x[run1]     = 0.0;
00537         iseed[run1] = 0.0;
00538     }
00539 
00540     t     = 0.0;
00541 
00542 
00543     // BDF-METHOD:
00544     // -----------
00545     nstep     = 5;
00546     psi       = (double**)calloc(1,sizeof(double*));
00547     psi_      = new double[nstep];
00548     gamma     = (double**)calloc(1,sizeof(double*));
00549 
00550     c2.init(m);
00551     c2.setZero();
00552 
00553     nOfNewtonSteps = (int*)calloc(1,sizeof(int));
00554     eta            = new double*[4];
00555     eta2           = new double*[4];
00556 
00557     nOfNewtonSteps[0] = 0;
00558 
00559     psi   [0] = (double*)calloc(nstep,sizeof(double));
00560     gamma [0] = (double*)calloc(nstep,sizeof(double));
00561 
00562     nablaY.init( nstep, m );
00563     nablaY.setZero();
00564 
00565     nablaY_.init( nstep, m );
00566     nablaY_.setZero();
00567 
00568     phi.init( nstep, m );
00569     phi.setZero();
00570 
00571     delta.init( nstep, m );
00572     delta.setZero();
00573 
00574     for( run1 = 0; run1 < nstep; run1++ ){
00575 
00576          psi_    [run1] = 0.0;
00577          psi  [0][run1] = 0.0;
00578          gamma[0][run1] = 0.0;
00579     }
00580 
00581 
00582     maxNM = 1;
00583     M     = (DMatrix**)calloc(maxNM,sizeof(DMatrix*));
00584     M_index   = (int*)calloc(maxNM,sizeof(int));
00585 
00586     M_index[0] = 0;
00587     M      [0] = 0;
00588     nOfM       = 0;
00589 
00590     las = arg.las;
00591 
00592     for( run1 = 0; run1 < 4; run1++ ){
00593         eta [run1] = new double[m];
00594         eta2[run1] = new double[m];
00595         for( run2 = 0; run2 < m; run2++ ){
00596             eta [run1][run2] = 0.0;
00597             eta2[run1][run2] = 0.0;
00598         }
00599     }
00600 
00601     F  = new double[m];
00602     F2 = new double[m];
00603 
00604 
00605     // SETTINGS:
00606     // ---------
00607     h    = (double*)calloc(1,sizeof(double));
00608     h[0]  = 0.001    ;
00609     hini  = 0.001    ;
00610     hmin = 0.000001  ;
00611     hmax = 1.0e10    ;
00612 
00613     tune  = 0.5      ;
00614     TOL   = 0.000001 ;
00615 
00616 
00617     // INTERNAL INDEX LISTS:
00618     // ---------------------
00619     diff_index = new int[m];
00620 
00621     for( run1 = 0; run1 < md; run1++ ){
00622         diff_index[run1] = rhs->index( VT_DIFFERENTIAL_STATE, run1 );
00623         if( diff_index[run1] == rhs->getNumberOfVariables() ){
00624             diff_index[run1] = diff_index[run1] + 1 + run1;
00625         }
00626     }
00627 
00628     ddiff_index = new int[md];
00629 
00630     for( run1 = 0; run1 < md; run1++ ){
00631         ddiff_index[run1] = rhs->index( VT_DDIFFERENTIAL_STATE, run1 );
00632         if( ddiff_index[run1] == rhs->getNumberOfVariables() ){
00633             ddiff_index[run1] = ddiff_index[run1] + 1 + md + run1;
00634         }
00635     }
00636 
00637     alg_index = new int[ma];
00638 
00639     for( run1 = 0; run1 < ma; run1++ ){
00640         alg_index [run1]    = rhs->index( VT_ALGEBRAIC_STATE, run1 );
00641         diff_index[md+run1] = alg_index [run1];
00642     }
00643 
00644     control_index       = new int[mu ];
00645 
00646     for( run1 = 0; run1 < mu; run1++ ){
00647         control_index[run1] = rhs->index( VT_CONTROL, run1 );
00648     }
00649 
00650     parameter_index     = new int[mp ];
00651 
00652     for( run1 = 0; run1 < mp; run1++ ){
00653         parameter_index[run1] = rhs->index( VT_PARAMETER, run1 );
00654     }
00655 
00656     int_control_index   = new int[mui];
00657 
00658     for( run1 = 0; run1 < mui; run1++ ){
00659         int_control_index[run1] = rhs->index( VT_INTEGER_CONTROL, run1 );
00660     }
00661 
00662     int_parameter_index = new int[mpi];
00663 
00664     for( run1 = 0; run1 < mpi; run1++ ){
00665         int_parameter_index[run1] = rhs->index( VT_INTEGER_PARAMETER, run1 );
00666     }
00667 
00668     disturbance_index   = new int[mw ];
00669 
00670     for( run1 = 0; run1 < mw; run1++ ){
00671         disturbance_index[run1] = rhs->index( VT_DISTURBANCE, run1 );
00672     }
00673 
00674     time_index = rhs->index( VT_TIME, 0 );
00675 
00676 
00677     // OTHERS:
00678     // -------
00679     maxNumberOfSteps = 1000;
00680     count            = 0   ;
00681     count2           = 0   ;
00682     count3           = 0   ;
00683 
00684     diff_scale.init(m);
00685 
00686     for( run1 = 0; run1 < md; run1++ ){
00687         diff_scale(run1) = rhs->scale( VT_DIFFERENTIAL_STATE, run1 );
00688     }
00689     for( run1 = 0; run1 < ma; run1++ ){
00690         diff_scale(md+run1) = rhs->scale( VT_ALGEBRAIC_STATE, run1 );
00691     }
00692 
00693     initial_guess    = new double[m];
00694 
00695     for( run1 = 0; run1 < m; run1++ ){
00696 
00697         initial_guess[run1] = 0.0;
00698     }
00699 
00700 
00701     // PRINT-LEVEL:
00702     // ------------
00703     PrintLevel = LOW;
00704 
00705 
00706     // SENSITIVITIES:
00707     // --------------
00708     nFDirs     = 0   ;
00709     nBDirs     = 0   ;
00710 
00711     nFDirs2    = 0   ;
00712     nBDirs2    = 0   ;
00713 
00714     G          = NULL;
00715     etaG       = NULL;
00716 
00717     G2         = NULL;
00718     G3         = NULL;
00719     etaG2      = NULL;
00720     etaG3      = NULL;
00721 
00722     c2G.init(md);
00723     c2G2.init(md);
00724     c2G3.init(md);
00725 
00726     phiG.init(nstep,m);
00727     phiG2.init(nstep,m);
00728     phiG3.init(nstep,m);
00729 
00730     deltaG.init(nstep,m);
00731     deltaG2.init(nstep,m);
00732     deltaG3.init(nstep,m);
00733 
00734 
00735     H          = NULL;
00736     etaH       = NULL;
00737 
00738     kH         = new double**[dim];
00739     for( run1 = 0; run1 < dim; run1++ ){
00740          kH[run1] = new double*[nstep];
00741          for( run2 = 0; run2 < nstep; run2++ ){
00742              kH[run1][run2] = new double[ndir];
00743          }
00744     }
00745     zH         = new double*[6];
00746     for( run1 = 0; run1 < 6; run1++ ){
00747          zH[run1] = new double[ndir];
00748     }
00749 
00750     H2         = NULL;
00751     H3         = NULL;
00752     etaH2      = NULL;
00753     etaH3      = NULL;
00754 
00755     kH2        = new double**[dim];
00756     for( run1 = 0; run1 < dim; run1++ ){
00757          kH2[run1] = new double*[nstep];
00758          for( run2 = 0; run2 < nstep; run2++ ){
00759              kH2[run1][run2] = new double[ndir];
00760          }
00761     }
00762     zH2        = new double*[6];
00763     for( run1 = 0; run1 < 6; run1++ ){
00764          zH2[run1] = new double[ndir];
00765     }
00766     kH3        = new double**[dim];
00767     for( run1 = 0; run1 < dim; run1++ ){
00768          kH3[run1] = new double*[nstep];
00769          for( run2 = 0; run2 < nstep; run2++ ){
00770              kH3[run1][run2] = new double[ndir];
00771          }
00772     }
00773     zH3        = new double*[6];
00774     for( run1 = 0; run1 < 6; run1++ ){
00775          zH3[run1] = new double[ndir];
00776     }
00777 
00778     // THE STATE OF AGGREGATION:
00779     // -------------------------
00780     soa        = SOA_UNFROZEN;
00781 
00782 
00783     // STORAGE:
00784     // --------
00785     maxAlloc = 1;
00786 
00787     nFcnEvaluations = 0;
00788     nJacEvaluations = 0;
00789 }
00790 
00791 
00792 Integrator* IntegratorBDF::clone() const{
00793 
00794     return new IntegratorBDF(*this);
00795 }
00796 
00797 
00798 void IntegratorBDF::deleteAll(){
00799 
00800     int run1, run2;
00801 
00802     if( initialAlgebraicResiduum != 0 )
00803         delete[] initialAlgebraicResiduum;
00804 
00805 
00806     // BUTCHER-
00807     // TABLEAU:
00808     // ----------
00809     for( run1 = 0; run1 < dim; run1++ ){
00810         delete[] A[run1];
00811     }
00812 
00813     delete[] A;
00814     delete[] b4;
00815     delete[] b5;
00816     delete[] c ;
00817 
00818 
00819     // RK-ALGORITHM:
00820     // -------------
00821     if( eta4 != NULL ){
00822         delete[] eta4;
00823     }
00824     if( eta4_ != NULL ){
00825         delete[] eta4_;
00826     }
00827     if( eta5 != NULL ){
00828         delete[] eta5;
00829     }
00830     if( eta5_ != NULL ){
00831         delete[] eta5_;
00832     }
00833 
00834 
00835 
00836     for( run2 = 0; run2 < 4; run2++ ){
00837 
00838          for( run1 = 0; run1 < dim; run1++ ){
00839             if( k[run2]  != NULL )
00840                 delete[] k[run2][run1] ;
00841             if( k2[run2] != NULL )
00842                 delete[] k2[run2][run1];
00843             if( l[run2]  != NULL )
00844                 delete[] l[run2][run1] ;
00845             if( l2[run2] != NULL )
00846                 delete[] l2[run2][run1];
00847          }
00848 
00849         if( k != NULL )
00850             delete[] k[run2] ;
00851 
00852         if( k2 != NULL )
00853             delete[] k2[run2];
00854 
00855         if( l != NULL )
00856             delete[] l[run2] ;
00857 
00858         if( l2 != NULL )
00859             delete[] l2[run2];
00860    }
00861 
00862 
00863     if( k != NULL )
00864         delete[] k;
00865 
00866     if( k2!= NULL )
00867         delete[] k2;
00868 
00869     if( l != NULL )
00870         delete[] l ;
00871 
00872     if( l2!= NULL )
00873         delete[] l2;
00874 
00875     if( iseed != NULL ){
00876         delete[] iseed;
00877     }
00878 
00879     if( x != NULL )
00880         delete[] x;
00881 
00882 
00883 
00884     for( run1 = 0; run1 < maxAlloc; run1++ ){
00885 
00886         if( psi[run1] != NULL ){
00887             free(psi[run1]);
00888         }
00889         if( gamma[run1] != NULL ){
00890             free(gamma[run1]);
00891         }
00892     }
00893 
00894     if( psi != NULL ){
00895         free(psi);
00896     }
00897 
00898     if( psi_ != NULL ){
00899         delete[] psi_;
00900     }
00901 
00902     if( gamma != NULL )
00903         free(gamma);
00904 
00905     if( nOfNewtonSteps != NULL )
00906         free(nOfNewtonSteps);
00907 
00908     for( run1 = 0; run1 < 4; run1++ ){
00909 
00910         if( eta != NULL )
00911             delete[] eta[run1];
00912         if( eta2 != NULL )
00913             delete[] eta2[run1];
00914     }
00915 
00916     delete[] eta    ;
00917     delete[] eta2   ;
00918 
00919     for( run1 = 0; run1 < maxNM; run1++ ){
00920          if( M[run1] != 0  )
00921              delete M[run1];
00922     }
00923 
00924     if( M != NULL ){
00925         free(M);
00926         free(M_index);
00927     }
00928 
00929     if( F != NULL )
00930         delete[] F;
00931     if( F2 != NULL )
00932         delete[] F2;
00933 
00934 
00935     // OTHERS:
00936     // -------
00937     if( initial_guess != NULL ){
00938         delete[] initial_guess;
00939     }
00940 
00941     // SENSITIVITIES:
00942     // --------------
00943 
00944     if( G  != NULL )
00945         delete[] G;
00946 
00947     if( etaG  != NULL )
00948         delete[] etaG;
00949 
00950     if( G2  != NULL )
00951         delete[] G2;
00952 
00953     if( G3  != NULL )
00954         delete[] G3;
00955 
00956     if( etaG2  != NULL )
00957         delete[] etaG2;
00958 
00959     if( etaG3  != NULL )
00960         delete[] etaG3;
00961 
00962 
00963     // ----------------------------------------
00964 
00965     if( H  != NULL )
00966         delete[] H;
00967 
00968     // ----------------------------------------
00969 
00970     if( H2  != NULL )
00971         delete[] H2;
00972     if( H3  != NULL )
00973         delete[] H3;
00974 
00975     for( run1 = 0; run1 < nstep; run1++ ){
00976         if( etaH != NULL )
00977             delete[] etaH[run1];
00978         if( etaH2 != NULL )
00979             delete[] etaH2[run1];
00980         if( etaH3 != NULL )
00981             delete[] etaH3[run1];
00982     }
00983     if( etaH != NULL )
00984         delete[] etaH;
00985     if( etaH2 != NULL )
00986         delete[] etaH2;
00987     if( etaH3 != NULL )
00988         delete[] etaH3;
00989 
00990 
00991     for( run1 = 0; run1 < dim; run1++ ){
00992         for( run2 = 0; run2 < nstep; run2++ ){
00993             if( kH != NULL )
00994                 delete[] kH[run1][run2];
00995             if( kH2 != NULL )
00996                 delete[] kH2[run1][run2];
00997             if( kH3 != NULL )
00998                 delete[] kH3[run1][run2];
00999         }
01000         if( kH != NULL )
01001             delete[] kH[run1];
01002         if( kH2 != NULL )
01003             delete[] kH2[run1];
01004         if( kH3 != NULL )
01005             delete[] kH3[run1];
01006     }
01007     if( kH != NULL )
01008         delete[] kH;
01009     if( kH2 != NULL )
01010         delete[] kH2;
01011     if( kH3 != NULL )
01012         delete[] kH3;
01013 
01014 
01015     for( run1 = 0; run1 < 6; run1++ ){
01016         if( zH != NULL )
01017             delete[] zH[run1];
01018         if( zH2 != NULL )
01019             delete[] zH2[run1];
01020         if( zH3 != NULL )
01021             delete[] zH3[run1];
01022     }
01023     if( zH != NULL )
01024         delete[] zH;
01025     if( zH2 != NULL )
01026         delete[] zH2;
01027     if( zH3 != NULL )
01028         delete[] zH3;
01029 }
01030 
01031 
01032 returnValue IntegratorBDF::freezeMesh(){
01033 
01034 
01035     if( soa != SOA_UNFROZEN ){
01036        if( PrintLevel != NONE ){
01037            return ACADOWARNING(RET_ALREADY_FROZEN);
01038        }
01039        return RET_ALREADY_FROZEN;
01040     }
01041 
01042     soa = SOA_FREEZING_MESH;
01043     return SUCCESSFUL_RETURN;
01044 }
01045 
01046 
01047 returnValue IntegratorBDF::freezeAll(){
01048 
01049     if( soa != SOA_UNFROZEN ){
01050        if( PrintLevel != NONE ){
01051            return ACADOWARNING(RET_ALREADY_FROZEN);
01052        }
01053        return RET_ALREADY_FROZEN;
01054     }
01055 
01056     soa = SOA_FREEZING_ALL;
01057     return SUCCESSFUL_RETURN;
01058 }
01059 
01060 
01061 returnValue IntegratorBDF::unfreeze(){
01062 
01063     int run1, run2;
01064 
01065     for( run1 = 1; run1 < maxAlloc; run1++ ){
01066 
01067         if( psi[run1] != NULL ){
01068             free(psi[run1]);
01069         }
01070         if( gamma[run1] != NULL ){
01071             free(gamma[run1]);
01072         }
01073     }
01074 
01075     maxAlloc = 1;
01076 
01077     psi            = (double**)realloc(psi,maxAlloc*sizeof(double*));
01078     gamma          = (double**)realloc(gamma,maxAlloc*sizeof(double*));
01079     nOfNewtonSteps = (int*)realloc(nOfNewtonSteps,(maxAlloc)*sizeof(int));
01080 
01081     for( run1 = 0; run1 < maxAlloc; run1++ ){
01082         nOfNewtonSteps[run1] = 0;
01083     }
01084 
01085     for( run1 = 0; run1 < maxAlloc; run1++ ){
01086 
01087         for( run2 = 0; run2 < 5; run2++ ){
01088              psi  [run1][run2] = 0.0;
01089              gamma[run1][run2] = 0.0;
01090         }
01091     }
01092 
01093     for( run1 = 0; run1 < maxNM; run1++ ){
01094          if( M[run1] != 0  )
01095              delete M[run1];
01096          M[run1] = 0;
01097     }
01098 
01099     maxNM = 1;
01100     nOfM  = 0;
01101     M       = (DMatrix**)realloc(M,maxNM*sizeof(DMatrix*));
01102     M_index = (int*)realloc(M_index,maxAlloc*sizeof(int));
01103 
01104     h = (double*)realloc(h,maxAlloc*sizeof(double));
01105 
01106     soa = SOA_UNFROZEN;
01107 
01108     return SUCCESSFUL_RETURN;
01109 }
01110 
01111 
01112 returnValue IntegratorBDF::evaluate( const DVector &x0  ,
01113                                      const DVector &xa  ,
01114                                      const DVector &p   ,
01115                                      const DVector &u   ,
01116                                      const DVector &w   ,
01117                                      const Grid   &t_    ){
01118 
01119     int         run1;
01120     returnValue returnvalue;
01121 
01122     if( rhs == NULL ){
01123         return ACADOERROR(RET_TRIVIAL_RHS);
01124     }
01125 
01126     Integrator::initializeOptions();
01127 
01128     timeInterval = t_;
01129 
01130     xStore.init( md+ma, timeInterval );
01131     iStore.init( mn   , timeInterval );
01132 
01133     t             = timeInterval.getFirstTime();
01134     x[time_index] = timeInterval.getFirstTime();
01135 
01136     if( x0.isEmpty() == BT_TRUE ) return ACADOERROR(RET_MISSING_INPUTS);
01137 
01138     if( (int) x0.getDim() < md )
01139         return ACADOERROR(RET_INPUT_HAS_WRONG_DIMENSION);
01140 
01141     for( run1 = 0; run1 < md; run1++ ){
01142         eta4[run1]          = x0(run1);
01143         eta5[run1]          = x0(run1);
01144         xStore(0,run1)      = x0(run1);
01145         x[diff_index[run1]] = x0(run1);
01146     }
01147 
01148     if( soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN  ){
01149        h[0] = hini;
01150 
01151        if( timeInterval.getIntervalLength() - h[0] < EPS ){
01152            h[0] = timeInterval.getIntervalLength();
01153        }
01154 
01155        if( h[0] < 10.0*EPS )
01156            return ACADOERROR(RET_TO_SMALL_OR_NEGATIVE_TIME_INTERVAL);
01157     }
01158 
01159     if( ma > 0 ){
01160 
01161         if( (int) xa.getDim() < ma )
01162             return ACADOERROR(RET_INPUT_HAS_WRONG_DIMENSION);
01163 
01164         for( run1 = 0; run1 < ma; run1++ ){
01165             eta4[md+run1]          = xa(run1);
01166             eta5[md+run1]          = xa(run1);
01167             xStore(0,md+run1)      = xa(run1);
01168             x[diff_index[md+run1]] = xa(run1);
01169             initial_guess[md+run1] = xa(run1);
01170         }
01171     }
01172 
01173     if( nFDirs != 0 )
01174         for( run1 = 0; run1 < md; run1++ )
01175             etaG[run1] = fseed(diff_index[run1]);
01176 
01177     if( mp > 0 ){
01178 
01179         if( (int) p.getDim() < mp )
01180             return ACADOERROR(RET_INPUT_HAS_WRONG_DIMENSION);
01181 
01182         for( run1 = 0; run1 < mp; run1++ ){
01183             x[parameter_index[run1]] = p(run1);
01184         }
01185     }
01186 
01187     if( mu > 0 ){
01188 
01189         if( (int) u.getDim() < mu )
01190             return ACADOERROR(RET_INPUT_HAS_WRONG_DIMENSION);
01191 
01192         for( run1 = 0; run1 < mu; run1++ ){
01193             x[control_index[run1]] = u(run1);
01194         }
01195     }
01196 
01197 
01198     if( mw > 0 ){
01199 
01200         if( (int) w.getDim() < mw )
01201             return ACADOERROR(RET_INPUT_HAS_WRONG_DIMENSION);
01202 
01203         for( run1 = 0; run1 < mw; run1++ ){
01204             x[disturbance_index[run1]] = w(run1);
01205         }
01206     }
01207 
01208         // log initial step
01209         logCurrentIntegratorStep( x0,xa );
01210 
01211      // start the time measurement:
01212      // ---------------------------
01213 
01214         totalTime.start();
01215         nFcnEvaluations = 0;
01216         nJacEvaluations = 0;
01217 
01218 
01219      // initialize the scaling based on the initial states:
01220      // ---------------------------------------------------
01221 
01222         double atol;
01223         get( ABSOLUTE_TOLERANCE, atol );
01224 
01225         for( run1 = 0; run1 < m; run1++ )
01226             diff_scale(run1) = fabs(eta4[run1]) + atol/TOL;
01227 
01228 
01229     returnvalue = rhs[0].evaluate( 0, x, initialAlgebraicResiduum );
01230 
01231     if( returnvalue != SUCCESSFUL_RETURN ){
01232 
01233         totalTime.stop();
01234         return ACADOWARNING(returnvalue);
01235     }
01236 
01237      // PRINTING:
01238      // ---------
01239         if( PrintLevel == MEDIUM ){
01240             acadoPrintCopyrightNotice( "IntegratorBDF -- A BDF integrator (order 4)." );
01241         }
01242         if( PrintLevel == HIGH ){
01243             cout << "BDF: t = " << t << "\t";
01244             for( run1 = 0; run1 < md; run1++ ){
01245                 cout << "x[" << run1 << "] = " << scientific << eta4[run1] << "  ";
01246             }
01247             for( run1 = 0; run1 < ma; run1++ ){
01248                 cout << "xq[" << run1 << "] = " << eta4[md+run1] << "  ";
01249             }
01250             cout << endl;
01251         }
01252 
01253 
01254     count3 = 0;
01255 
01256     returnvalue = rk_start();
01257 
01258     if( returnvalue != SUCCESSFUL_RETURN ){
01259 
01260         totalTime.stop();
01261 
01262         if( PrintLevel != NONE )
01263             return ACADOERROR(returnvalue);
01264 
01265         return returnvalue;
01266     }
01267 
01268     if( soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN  ){
01269 
01270        if( timeInterval.getLastTime() - t - h[0] < EPS ){
01271            h[0] = timeInterval.getLastTime() - t;
01272        }
01273     }
01274 
01275     returnvalue = RET_FINAL_STEP_NOT_PERFORMED_YET;
01276 
01277     count = 7;
01278 
01279     while( returnvalue == RET_FINAL_STEP_NOT_PERFORMED_YET && count <= maxNumberOfSteps ){
01280 
01281         returnvalue = step(count);
01282         count++;
01283     }
01284 
01285         // log final step
01286         logCurrentIntegratorStep( );
01287 
01288     count2 = count-1;
01289 
01290 
01291     for( run1 = 0; run1 < mn; run1++ )
01292         iStore( 0, run1 ) = iStore( 1, run1 );
01293 
01294 
01295     // stop the measurement of total time:
01296     // -----------------------------------
01297        totalTime.stop();
01298 
01299            
01300 //         cout <<"numIntSteps = %d\n",count-1);
01301            
01302 
01303     // SET THE LOGGING INFORMATION:
01304     // ----------------------------------------------------------------------------------------
01305 
01306        setLast( LOG_TIME_INTEGRATOR                              , totalTime.getTime()           );
01307        setLast( LOG_NUMBER_OF_INTEGRATOR_STEPS                   , count-1                       );
01308        setLast( LOG_NUMBER_OF_INTEGRATOR_REJECTED_STEPS          , getNumberOfRejectedSteps()    );
01309        setLast( LOG_NUMBER_OF_INTEGRATOR_FUNCTION_EVALUATIONS    , nFcnEvaluations               );
01310        setLast( LOG_NUMBER_OF_BDF_INTEGRATOR_JACOBIAN_EVALUATIONS, nJacEvaluations               );
01311        setLast( LOG_TIME_INTEGRATOR_FUNCTION_EVALUATIONS         , functionEvaluation.getTime()  );
01312        setLast( LOG_TIME_BDF_INTEGRATOR_JACOBIAN_EVALUATION      , jacComputation.getTime()      );
01313        setLast( LOG_TIME_BDF_INTEGRATOR_JACOBIAN_DECOMPOSITION   , jacDecomposition.getTime()    );
01314 
01315     // ----------------------------------------------------------------------------------------
01316 
01317 
01318 
01319     if( count > maxNumberOfSteps ){
01320 
01321         if( PrintLevel != NONE )
01322             return ACADOERROR(RET_MAX_NUMBER_OF_STEPS_EXCEEDED);
01323         return RET_MAX_NUMBER_OF_STEPS_EXCEEDED;
01324     }
01325 
01326 
01327      // PRINTING:
01328      // ---------
01329         if( PrintLevel == MEDIUM ){
01330 
01331             cout << "\n Results at  t =  " << t << "   : \n\n";
01332             for( run1 = 0; run1 < md; run1++ ){
01333                 cout << "x[" << run1 << "] = " << scientific << nablaY(0,run1) << "  ";
01334             }
01335             for( run1 = 0; run1 < ma; run1++ ){
01336                 cout << "x[" << run1 << "] = " << scientific << nablaY(0,md+run1) << "  ";
01337             }
01338             cout << endl;
01339             printBDFfinalResults();
01340         }
01341                 
01342         int printIntegratorProfile = 0;
01343         get( PRINT_INTEGRATOR_PROFILE,printIntegratorProfile );
01344         
01345         if ( (BooleanType)printIntegratorProfile == BT_TRUE )
01346         {
01347                 printRunTimeProfile( );
01348         }
01349         else
01350         {
01351                 if( PrintLevel == MEDIUM  || PrintLevel == HIGH )
01352                         cout << "BDF: number of steps:  " << count - 1 << endl;
01353         }
01354 
01355     return returnvalue;
01356 }
01357 
01358 
01359 returnValue IntegratorBDF::setProtectedForwardSeed( const DVector &xSeed ,
01360                                                     const DVector &pSeed ,
01361                                                     const DVector &uSeed ,
01362                                                     const DVector &wSeed ,
01363                                                     const int    &order   ){
01364 
01365     if( order == 2 ){
01366         return setForwardSeed2( xSeed, pSeed, uSeed, wSeed );
01367     }
01368     if( order < 1 || order > 2 ){
01369         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
01370     }
01371 
01372     if( nBDirs > 0 ){
01373         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
01374     }
01375 
01376     int run2;
01377 
01378     if( G  != NULL ){
01379         delete[] G;
01380         G = NULL;
01381     }
01382 
01383     if( etaG  != NULL ){
01384         delete[] etaG;
01385         etaG = NULL;
01386     }
01387 
01388     nFDirs = 1;
01389 
01390     fseed.init(ndir);
01391     fseed.setZero();
01392 
01393     G = new double[ndir];
01394 
01395     for( run2 = 0; run2 < (ndir); run2++ )
01396          G[run2] = 0.0;
01397 
01398     etaG = new double[m];
01399 
01400     nablaG.init( nstep, m );
01401 
01402 
01403     if( xSeed.getDim() != 0 )
01404         for( run2 = 0; run2 < md; run2++ )
01405             fseed(diff_index[run2]) = xSeed(run2);
01406 
01407     if( pSeed.getDim() != 0 )
01408         for( run2 = 0; run2 < mp; run2++ ){
01409             fseed(parameter_index[run2]) = pSeed(run2);
01410             G    [parameter_index[run2]] = pSeed(run2);
01411         }
01412 
01413     if( uSeed.getDim() != 0 )
01414         for( run2 = 0; run2 < mu; run2++ ){
01415             fseed(control_index[run2]) = uSeed(run2);
01416             G    [control_index[run2]] = uSeed(run2);
01417         }
01418 
01419     if( wSeed.getDim() != 0 )
01420         for( run2 = 0; run2 < mw; run2++ ){
01421             fseed(disturbance_index[run2]) = wSeed(run2);
01422             G    [disturbance_index[run2]] = wSeed(run2);
01423         }
01424 
01425     return SUCCESSFUL_RETURN;
01426 }
01427 
01428 
01429 returnValue IntegratorBDF::setForwardSeed2( const DVector &xSeed ,
01430                                             const DVector &pSeed ,
01431                                             const DVector &uSeed ,
01432                                             const DVector &wSeed   ){
01433 
01434     int run2;
01435 
01436 
01437     if( G2  != NULL ){
01438         delete[] G2;
01439         G2 = NULL;
01440     }
01441 
01442     if( G3  != NULL ){
01443         delete[] G3;
01444         G3 = NULL;
01445     }
01446 
01447     if( etaG2  != NULL ){
01448         delete[] etaG2;
01449         etaG2 = NULL;
01450     }
01451 
01452     if( etaG3  != NULL ){
01453         delete[] etaG3;
01454         etaG3 = NULL;
01455     }
01456 
01457 
01458     nFDirs2 = 1;
01459 
01460     fseed2.init(ndir);
01461     fseed2.setZero();
01462 
01463     G2 = new double[ndir];
01464     G3 = new double[ndir];
01465 
01466     for( run2 = 0; run2 < (ndir); run2++ ){
01467          G2[run2] = 0.0;
01468          G3[run2] = 0.0;
01469     }
01470 
01471     etaG2 = new double[m];
01472     etaG3 = new double[m];
01473 
01474     nablaG2.init( nstep, m );
01475     nablaG3.init( nstep, m );
01476 
01477     if( xSeed.getDim() != 0 )
01478         for( run2 = 0; run2 < md; run2++ )
01479             fseed2(diff_index[run2]) = xSeed(run2);
01480 
01481     if( pSeed.getDim() != 0 )
01482         for( run2 = 0; run2 < mp; run2++ ){
01483             fseed2(parameter_index[run2]) = pSeed(run2);
01484             G2    [parameter_index[run2]] = pSeed(run2);
01485         }
01486 
01487     if( uSeed.getDim() != 0 )
01488         for( run2 = 0; run2 < mu; run2++ ){
01489             fseed2(control_index[run2]) = uSeed(run2);
01490             G2    [control_index[run2]] = uSeed(run2);
01491         }
01492 
01493     if( wSeed.getDim() != 0 )
01494         for( run2 = 0; run2 < mw; run2++ ){
01495             fseed2(disturbance_index[run2]) = wSeed(run2);
01496             G2    [disturbance_index[run2]] = wSeed(run2);
01497         }
01498 
01499     return SUCCESSFUL_RETURN;
01500 }
01501 
01502 
01503 returnValue IntegratorBDF::setProtectedBackwardSeed( const DVector &seed, const int &order ){
01504 
01505     if( order == 2 ){
01506         return setBackwardSeed2(seed);
01507     }
01508     if( order < 1 || order > 2 ){
01509         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
01510     }
01511 
01512     if( nFDirs > 0 ){
01513         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
01514     }
01515 
01516     int run1, run2, run3;
01517 
01518     if( H  != NULL ){
01519         delete[] H;
01520         H = NULL;
01521     }
01522 
01523     for( run1 = 0; run1 < nstep; run1++ ){
01524         if( etaH != NULL )
01525             delete[] etaH[run1];
01526     }
01527     if( etaH != NULL ){
01528         delete[] etaH;
01529         etaH = 0;
01530     }
01531 
01532     nBDirs = 1;
01533 
01534     bseed.init(m);
01535     bseed.setZero();
01536     H = new double[m];
01537 
01538     if( seed.getDim() != 0 )
01539         for( run2 = 0; run2 < md; run2++ )
01540             bseed(run2) = seed(run2);
01541 
01542     etaH = new double*[nstep];
01543     for( run3 = 0; run3 < nstep; run3++ ){
01544         etaH[run3] = new double[ndir];
01545         for( run2 = 0; run2 < ndir; run2++ ){
01546             etaH[run3][run2] = 0.0;
01547         }
01548     }
01549 
01550     nablaH.init( nstep, ndir );
01551     nablaH.setZero();
01552 
01553     nablaH_.init( nstep, ndir );
01554     nablaH_.setZero();
01555 
01556     deltaH.init( nstep, ndir );
01557     deltaH.setZero();
01558 
01559     c2H.init(ndir);
01560     c2H.setZero();
01561 
01562     return SUCCESSFUL_RETURN;
01563 }
01564 
01565 
01566 returnValue IntegratorBDF::setBackwardSeed2( const DVector &seed ){
01567 
01568     int run1, run2, run3;
01569 
01570     if( H2  != NULL ){
01571         delete[] H2;
01572         H2 = NULL;
01573     }
01574     if( H3  != NULL ){
01575         delete[] H3;
01576         H3 = NULL;
01577     }
01578 
01579     for( run1 = 0; run1 < nstep; run1++ ){
01580         if( etaH2 != NULL )
01581             delete[] etaH2[run1];
01582         if( etaH3 != NULL )
01583             delete[] etaH3[run1];
01584     }
01585     if( etaH2 != NULL ){
01586         delete[] etaH2;
01587         etaH2 = 0;
01588     }
01589     if( etaH3 != NULL ){
01590         delete[] etaH3;
01591         etaH3 = 0;
01592     }
01593 
01594     nBDirs2 = 1;
01595 
01596     bseed2.init(m);
01597 
01598     H2 = new double[m];
01599     H3 = new double[m];
01600 
01601     if( seed.getDim() != 0 ){
01602         for( run2 = 0; run2 < md; run2++ ){
01603             bseed2(run2) = seed(run2);
01604         }
01605     }
01606 
01607     etaH2 = new double*[nstep];
01608     for( run3 = 0; run3 < nstep; run3++ ){
01609         etaH2[run3] = new double[ndir];
01610         for( run2 = 0; run2 < ndir; run2++ ){
01611              etaH2[run3][run2] = 0.0;
01612         }
01613     }
01614     etaH3 = new double*[nstep];
01615     for( run3 = 0; run3 < nstep; run3++ ){
01616         etaH3[run3] = new double[ndir];
01617         for( run2 = 0; run2 < ndir; run2++ ){
01618              etaH3[run3][run2] = 0.0;
01619         }
01620     }
01621 
01622     nablaH2.init( nstep, ndir );
01623     nablaH2.setZero();
01624 
01625     nablaH3.init( nstep, ndir );
01626     nablaH3.setZero();
01627 
01628     nablaH2_.init( nstep, ndir );
01629     nablaH2_.setZero();
01630 
01631     nablaH3_.init( nstep, ndir );
01632     nablaH3_.setZero();
01633 
01634     deltaH2.init( nstep, ndir );
01635     deltaH2.setZero();
01636 
01637     deltaH3.init( nstep, ndir );
01638     deltaH3.setZero();
01639 
01640     c2H2.init( ndir );
01641     c2H2.setZero();
01642 
01643     c2H3.init( ndir );
01644     c2H3.setZero();
01645 
01646     return SUCCESSFUL_RETURN;
01647 }
01648 
01649 
01650 returnValue IntegratorBDF::evaluateSensitivities(){
01651 
01652     int         run1;
01653     returnValue returnvalue;
01654 
01655     if( rhs == NULL ){
01656         return ACADOERROR(RET_TRIVIAL_RHS);
01657     }
01658 
01659     if( soa != SOA_EVERYTHING_FROZEN ){
01660         return ACADOERROR(RET_NOT_FROZEN);
01661     }
01662 
01663     if( nBDirs2 == 0 && nFDirs != 0 ){
01664         t = timeInterval.getFirstTime();
01665         dxStore.init ( md+ma, timeInterval );
01666         for( run1 = 0; run1 < md; run1++ ){
01667              etaG[run1] = fseed(diff_index[run1]);
01668         }
01669     }
01670 
01671     nablaH.setZero();
01672 
01673     if( nBDirs != 0 ){
01674         for( run1 = 0; run1 < md; run1++ ){
01675             nablaH(0,diff_index[run1]) = bseed(run1);
01676         }
01677     }
01678 
01679     if( nFDirs2 != 0 ){
01680         t = timeInterval.getFirstTime();
01681         ddxStore.init( md+ma, timeInterval );
01682         for( run1 = 0; run1 < md; run1++ ){
01683             etaG2[run1] = fseed2(diff_index[run1]);
01684             etaG3[run1] = 0.0;
01685         }
01686     }
01687 
01688     nablaH2.setZero();
01689     nablaH3.setZero();
01690 
01691     if( nBDirs2 != 0 ){
01692         for( run1 = 0; run1 < md; run1++ ){
01693             nablaH2(0,diff_index[run1]) = bseed2(run1);
01694         }
01695     }
01696 
01697 
01698     returnvalue = RET_FINAL_STEP_NOT_PERFORMED_YET;
01699 
01700     if( nBDirs > 0 || nBDirs2 > 0 ){
01701 
01702         t = timeInterval.getLastTime();
01703         h[0] = 1.0;
01704 
01705         int oldCount = count;
01706 
01707         count--;
01708         while( returnvalue == RET_FINAL_STEP_NOT_PERFORMED_YET && count >= 7 ){
01709 
01710             returnvalue = step( count );
01711             count--;
01712         }
01713 
01714         if( returnvalue != SUCCESSFUL_RETURN &&
01715             returnvalue != RET_FINAL_STEP_NOT_PERFORMED_YET ){
01716             count = oldCount;
01717             if( PrintLevel != NONE )
01718                 return ACADOERROR(returnvalue);
01719             return returnvalue;
01720         }
01721 
01722         returnvalue = rk_start();
01723 
01724         if( returnvalue == SUCCESSFUL_RETURN ){
01725 
01726            // PRINTING:
01727            // ---------
01728               if( PrintLevel == MEDIUM )
01729                   printBDFfinalResults();
01730 
01731             count = oldCount;
01732 
01733             return SUCCESSFUL_RETURN;
01734         }
01735         count = oldCount;
01736     }
01737     else{
01738 
01739         t = timeInterval.getFirstTime();
01740 
01741         returnvalue = rk_start();
01742 
01743         if( returnvalue != SUCCESSFUL_RETURN ){
01744             if( PrintLevel != NONE )
01745                 return ACADOERROR(returnvalue);
01746             return returnvalue;
01747         }
01748 
01749         returnvalue = RET_FINAL_STEP_NOT_PERFORMED_YET;
01750         count = 7;
01751         while( returnvalue == RET_FINAL_STEP_NOT_PERFORMED_YET &&
01752                count <= maxNumberOfSteps ){
01753 
01754             returnvalue = step(count);
01755             count++;
01756         }
01757 
01758         if( nBDirs2 == 0 && nFDirs != 0 )
01759             for( run1 = 0; run1 < m; run1++ )
01760                 dxStore( 0, run1 ) = dxStore( 1, run1 );
01761 
01762         if( nFDirs2 != 0 )
01763             for( run1 = 0; run1 < m; run1++ )
01764                 ddxStore( 0, run1 ) = ddxStore( 1, run1 );
01765 
01766         if( count > maxNumberOfSteps ){
01767             if( PrintLevel != NONE )
01768                 return ACADOERROR(RET_MAX_NUMBER_OF_STEPS_EXCEEDED);
01769             return RET_MAX_NUMBER_OF_STEPS_EXCEEDED;
01770         }
01771 
01772         // PRINTING:
01773         // ---------
01774            if( PrintLevel == MEDIUM )
01775                printBDFfinalResults();
01776 
01777         return SUCCESSFUL_RETURN;
01778     }
01779 
01780     if( PrintLevel != NONE )
01781         return ACADOERROR(returnvalue);
01782 
01783     return returnvalue;
01784 }
01785 
01786 
01787 returnValue IntegratorBDF::step(int number_){
01788 
01789     int run1;
01790     double E = EPS;
01791 
01792      if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
01793 
01794          int oldAlloc = maxAlloc;
01795 
01796          if( number_ >= maxAlloc ){
01797              maxAlloc = 2*maxAlloc;
01798 
01799              psi            = (double**)realloc(psi,maxAlloc*sizeof(double*));
01800              gamma          = (double**)realloc(gamma,maxAlloc*sizeof(double*));
01801              nOfNewtonSteps = (int*)realloc(nOfNewtonSteps,(maxAlloc)*sizeof(int));
01802 
01803              for( run1 = oldAlloc; run1 < maxAlloc; run1++ ){
01804                  nOfNewtonSteps[run1] = 0;
01805              }
01806              for( run1 = oldAlloc; run1 < maxAlloc; run1++ ){
01807                  psi  [run1] = (double*)calloc(nstep,sizeof(double));
01808                  gamma[run1] = (double*)calloc(nstep,sizeof(double));
01809              }
01810 
01811              M_index      = (int*)realloc(M_index,maxAlloc*sizeof(int));
01812              h = (double*)realloc(h,maxAlloc*sizeof(double));
01813          }
01814      }
01815 
01816 
01817     if( soa == SOA_EVERYTHING_FROZEN || soa == SOA_MESH_FROZEN ){
01818         rel_time_scale = h[0];
01819         h[0] = h[number_];
01820     }
01821 
01822     if( soa == SOA_FREEZING_MESH ||
01823         soa == SOA_MESH_FROZEN   ||
01824         soa == SOA_UNFROZEN      ){
01825 
01826         if( number_ == 7 ){
01827            E = determinePredictor(7, BT_TRUE );
01828         }
01829         else{
01830            E = determinePredictor(7, BT_FALSE );
01831         }
01832     }
01833     if( soa == SOA_FREEZING_ALL ){
01834 
01835         if( number_ == 7 ){
01836            E = determinePredictor(number_, BT_TRUE );
01837         }
01838         else{
01839            E = determinePredictor(number_, BT_FALSE );
01840         }
01841     }
01842 
01843 
01844     if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
01845 
01846         int number_of_rejected_steps = 0;
01847 
01848         if( E < 0.0 ){
01849             return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
01850         }
01851 
01852         // REJECT THE STEP IF GIVEN TOLERANCE IS NOT ACHIEVED:
01853         // ---------------------------------------------------
01854         while( E >= TOL ){
01855 
01856             if( PrintLevel == HIGH ){
01857 
01858                 cout << "STEP REJECTED: error estimate           = " << scientific << E
01859                          << "               required local tolerance = " <<  TOL << endl;
01860             }
01861 
01862             number_of_rejected_steps++;
01863 
01864             if( soa == SOA_FREEZING_MESH ||
01865                 soa == SOA_UNFROZEN      ){
01866 
01867                 psi[7][0] = psi_[0];
01868                 psi[7][1] = psi_[1];
01869                 psi[7][2] = psi_[2];
01870                 psi[7][3] = psi_[3];
01871             }
01872             if( soa == SOA_FREEZING_ALL ){
01873 
01874                 psi[number_][0] = psi_[0];
01875                 psi[number_][1] = psi_[1];
01876                 psi[number_][2] = psi_[2];
01877                 psi[number_][3] = psi_[3];
01878             }
01879 
01880 
01881             if( h[0] <= hmin + EPS ){
01882                 return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
01883             }
01884             h[0] = 0.5*h[0];
01885             if( E > 0.9*INFTY ) h[0] = 0.2*h[0];
01886 
01887             if( h[0] < hmin ){
01888                 h[0] = hmin;
01889             }
01890 
01891             if( soa == SOA_FREEZING_MESH ||
01892                 soa == SOA_UNFROZEN      ){
01893 
01894                 E = determinePredictor( 7, BT_FALSE );
01895             }
01896             if( soa == SOA_FREEZING_ALL ){
01897 
01898                 E = determinePredictor( number_, BT_FALSE );
01899             }
01900 
01901             if( E < 0.0 ){
01902                 return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
01903             }
01904         }
01905 
01906         count3 += number_of_rejected_steps;
01907     }
01908 
01909     // PROCEED IF THE STEP IS ACCEPTED:
01910     // --------------------------------
01911 
01912     if( soa != SOA_EVERYTHING_FROZEN ) nablaY = nablaY_;
01913 
01914         // log current step
01915         logCurrentIntegratorStep( );
01916 
01917 
01918      // compute forward derivatives if requested:
01919      // ------------------------------------------
01920 
01921      if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 ){
01922 
01923          if( nBDirs != 0 ){
01924              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
01925          }
01926 
01927          if( soa == SOA_FREEZING_ALL || soa == SOA_EVERYTHING_FROZEN ){
01928              determineBDFEtaGForward(number_);
01929          }
01930          else{
01931              determineBDFEtaGForward(7);
01932          }
01933      }
01934      if( nBDirs > 0 ){
01935 
01936          if( soa != SOA_EVERYTHING_FROZEN ){
01937              return ACADOERROR(RET_NOT_FROZEN);
01938          }
01939          if( nFDirs != 0 || nBDirs2 != 0 || nFDirs2 != 0 ){
01940              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
01941          }
01942          determineBDFEtaHBackward(number_);
01943      }
01944      if( nFDirs2 > 0 ){
01945 
01946          if( soa != SOA_EVERYTHING_FROZEN ){
01947              return ACADOERROR(RET_NOT_FROZEN);
01948          }
01949          if( nBDirs != 0 || nBDirs2 != 0 || nFDirs != 1 ){
01950              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
01951          }
01952          determineBDFEtaGForward2(number_);
01953      }
01954      if( nBDirs2 > 0 ){
01955 
01956          if( soa != SOA_EVERYTHING_FROZEN ){
01957              return ACADOERROR(RET_NOT_FROZEN);
01958          }
01959          if( nBDirs != 0 || nFDirs2 != 0 || nFDirs != 1 ){
01960              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
01961          }
01962 
01963          determineBDFEtaHBackward2(number_);
01964      }
01965 
01966 
01967      // increase the time:
01968      // ----------------------------------------------
01969 
01970      if( nBDirs > 0 || nBDirs2 > 0 ){
01971 
01972          t = t - h[0];
01973      }
01974      else{
01975 
01976          t = t + h[0];
01977      }
01978 //      printf( "t = %e,  stepsize = %e\n", t,h[0] );
01979 
01980      // PRINTING:
01981      // ---------
01982         printBDFIntermediateResults();
01983 
01984 
01985      // STORAGE:
01986      // --------
01987 
01988      if( soa == SOA_FREEZING_MESH || soa == SOA_FREEZING_ALL ){
01989 
01990          if( number_ >= maxAlloc){
01991 
01992              maxAlloc = 2*maxAlloc;
01993              h = (double*)realloc(h,maxAlloc*sizeof(double));
01994          }
01995 
01996          h[number_] = h[0];
01997      }
01998 
01999      if( nFDirs == 0 && nBDirs  == 0 && nFDirs2 == 0 && nBDirs == 0 ){
02000 
02001          interpolate( number_, nablaY, xStore );
02002 
02003          int i1 = timeInterval.getFloorIndex( t             );
02004          int i2 = timeInterval.getFloorIndex( t + c[6]*h[0] );
02005          int jj;
02006 
02007          for( jj = i1+1; jj <= i2; jj++ )
02008              for( run1 = 0; run1 < mn; run1++ )
02009                  iStore( jj, run1 ) = x[rhs->index( VT_INTERMEDIATE_STATE, run1 )];
02010      }
02011      if( nFDirs  > 0 && nBDirs2 == 0 && nFDirs2 == 0 ) interpolate( number_, nablaG , dxStore  );
02012      if( nFDirs2 > 0                                 ) interpolate( number_, nablaG3, ddxStore );
02013 
02014 
02015      if( nBDirs == 0 || nBDirs2 == 0 ){
02016 
02017      // Stop the algorithm if  t >= te:
02018      // ----------------------------------------------
02019         if( t >= timeInterval.getLastTime() - EPS ){
02020             x[time_index] = timeInterval.getLastTime();
02021             for( run1 = 0; run1 < m; run1++ ){
02022                                 if ( acadoIsNaN( nablaY(0,run1) ) == BT_TRUE )
02023                                         return ACADOERROR( RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF );
02024                 x[diff_index[run1]] = nablaY(0,run1);
02025             }
02026 
02027             if( soa == SOA_FREEZING_MESH ){
02028                 soa = SOA_MESH_FROZEN;
02029             }
02030             if( soa == SOA_FREEZING_ALL ){
02031                 soa = SOA_EVERYTHING_FROZEN;
02032             }
02033 
02034             return SUCCESSFUL_RETURN;
02035         }
02036      }
02037 
02038 
02039      if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
02040 
02041 
02042      // recompute the scaling based on the actual states:
02043      // -------------------------------------------------
02044 
02045         double atol;
02046         get( ABSOLUTE_TOLERANCE, atol );
02047 
02048         for( run1 = 0; run1 < m; run1++ )
02049             diff_scale(run1) = fabs(nablaY(0,run1)) + atol/TOL;
02050 
02051 
02052      // apply a numeric stabilization of the step size control:
02053      // -------------------------------------------------------
02054         double Emin = 1e-3*sqrt(TOL)*pow( hini, 2.5 );
02055 
02056         if( E < Emin     ) E = Emin    ;
02057         if( E < 10.0*EPS ) E = 10.0*EPS;
02058 
02059 
02060 
02061      // determine the new step size:
02062      // ----------------------------------------------
02063         rel_time_scale = h[0];
02064         h[0] = h[0]*pow( tune*(TOL/E) , 0.20 );
02065 
02066         if( h[0] / rel_time_scale > 2.0 ) h[0] = 2.0*rel_time_scale;
02067 
02068         if( h[0] > hmax ){
02069           h[0] = hmax;
02070         }
02071         if( h[0] < hmin ){
02072           h[0] = hmin;
02073         }
02074 
02075         if( t + h[0] >= timeInterval.getLastTime() ){
02076           h[0] = timeInterval.getLastTime()-t;
02077         }
02078     }
02079 
02080     return RET_FINAL_STEP_NOT_PERFORMED_YET;
02081 }
02082 
02083 
02084 
02085 returnValue IntegratorBDF::stop(){
02086 
02087     return ACADOERROR(RET_NOT_IMPLEMENTED_YET);
02088 }
02089 
02090 
02091 
02092 returnValue IntegratorBDF::getProtectedX( DVector *xEnd ) const{
02093 
02094     int run1;
02095 
02096     if( (int) xEnd[0].getDim() != m )
02097         return RET_INPUT_HAS_WRONG_DIMENSION;
02098 
02099     for( run1 = 0; run1 < m; run1++ )
02100         xEnd[0](run1) = nablaY(0,run1);
02101 
02102     return SUCCESSFUL_RETURN;
02103 }
02104 
02105 
02106 returnValue IntegratorBDF::getProtectedForwardSensitivities( DMatrix *Dx, int order ) const{
02107 
02108     int run1;
02109 
02110     if( Dx == NULL ){
02111         return SUCCESSFUL_RETURN;
02112     }
02113 
02114     if( order == 1 && nFDirs2 == 0 ){
02115 
02116         if( (int) Dx[0].getNumCols() != nFDirs )
02117             return RET_INPUT_HAS_WRONG_DIMENSION;
02118         if( (int) Dx[0].getNumRows() != m )
02119             return RET_INPUT_HAS_WRONG_DIMENSION;
02120 
02121         for( run1 = 0; run1 < m; run1++ ){
02122             Dx[0](run1,0) = nablaG(0,run1);
02123         }
02124     }
02125 
02126     if( order == 2 ){
02127 
02128         if( (int) Dx[0].getNumCols() != nFDirs2 )
02129             return RET_INPUT_HAS_WRONG_DIMENSION;
02130         if( (int) Dx[0].getNumRows() != m )
02131             return RET_INPUT_HAS_WRONG_DIMENSION;
02132 
02133         for( run1 = 0; run1 < m; run1++ ){
02134             Dx[0](run1,0) = nablaG3(0,run1);
02135         }
02136     }
02137 
02138     if( order == 1 && nFDirs2 > 0 ){
02139 
02140         if( (int) Dx[0].getNumCols() != nFDirs2 )
02141             return RET_INPUT_HAS_WRONG_DIMENSION;
02142         if( (int) Dx[0].getNumRows() != m )
02143             return RET_INPUT_HAS_WRONG_DIMENSION;
02144 
02145         for( run1 = 0; run1 < m; run1++ ){
02146             Dx[0](run1,0) = nablaG2(0,run1);
02147         }
02148     }
02149 
02150     if( order < 1 || order > 2 ){
02151         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
02152     }
02153 
02154     return SUCCESSFUL_RETURN;
02155 }
02156 
02157 
02158 returnValue IntegratorBDF::getProtectedBackwardSensitivities(DVector &Dx_x0,
02159                                                              DVector &Dx_p ,
02160                                                              DVector &Dx_u ,
02161                                                              DVector &Dx_w ,
02162                                                              int order      ) const{
02163 
02164     if( order == 1 && nBDirs2 == 0 )
02165         copyBackward( Dx_x0, Dx_p, Dx_u, Dx_w, nablaH );
02166 
02167     if( order == 1 && nBDirs2 > 0 )
02168         copyBackward( Dx_x0, Dx_p, Dx_u, Dx_w, nablaH2 );
02169 
02170     if( order == 2 )
02171         copyBackward( Dx_x0, Dx_p, Dx_u, Dx_w, nablaH3 );
02172 
02173     if( order < 1 || order > 2 )
02174         return ACADOERROR(RET_INPUT_OUT_OF_RANGE);
02175 
02176     return SUCCESSFUL_RETURN;
02177 }
02178 
02179 
02180 
02181 returnValue IntegratorBDF::setDxInitialization( double *dx0 ){
02182 
02183     int run1;
02184 
02185     if( dx0 != NULL ){
02186 
02187         for( run1 = 0; run1 < md; run1++ ){
02188             initial_guess[run1] = dx0[run1];
02189         }
02190     }
02191     else{
02192 
02193         for( run1 = 0; run1 < md; run1++ ){
02194             initial_guess[run1] = 0.0;
02195         }
02196     }
02197 
02198     for( run1 = 0; run1 < md; run1++ ){
02199         initial_guess[md+run1] = 0.0;
02200     }
02201 
02202     return SUCCESSFUL_RETURN;
02203 }
02204 
02205 
02206 int IntegratorBDF::getNumberOfSteps() const{
02207 
02208     return count2-2;
02209 }
02210 
02211 int IntegratorBDF::getNumberOfRejectedSteps() const{
02212 
02213     return count3;
02214 }
02215 
02216 
02217 
02218 
02219 double IntegratorBDF::getStepSize() const{
02220 
02221     return h[0];
02222 }
02223 
02224 
02225 
02226 //
02227 // PROTECTED MEMBER FUNCTIONS:
02228 //
02229 
02230 
02231 double IntegratorBDF::determinePredictor( int number_, BooleanType ini ){
02232 
02233     int run1;
02234     double pp;
02235 
02236     if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
02237 
02238         if( soa != SOA_FREEZING_ALL ){
02239 
02240             psi_[0] = psi[number_][0];
02241             psi_[1] = psi[number_][1];
02242             psi_[2] = psi[number_][2];
02243             psi_[3] = psi[number_][3];
02244 
02245             psi[number_][3] = psi[number_][2] + h[0];
02246             psi[number_][2] = psi[number_][1] + h[0];
02247             psi[number_][1] = psi[number_][0] + h[0];
02248             psi[number_][0] = h[0];
02249 
02250             gamma[number_][4] =   1.0/psi[number_][0] + 1.0/psi[number_][1]
02251                                 + 1.0/psi[number_][2] + 1.0/psi[number_][3];
02252         }
02253         else{
02254 
02255             psi[number_][3] = psi[number_-1][2] + h[0];
02256             psi[number_][2] = psi[number_-1][1] + h[0];
02257             psi[number_][1] = psi[number_-1][0] + h[0];
02258             psi[number_][0] = h[0];
02259 
02260             gamma[number_][4] =   1.0/psi[number_][0] + 1.0/psi[number_][1]
02261                                 + 1.0/psi[number_][2] + 1.0/psi[number_][3];
02262         }
02263     }
02264 
02265     const double scalT = h[0]/rel_time_scale;
02266 
02267     pp = psi[number_][0]*scalT/h[0];
02268     for( run1 = 0; run1 < m; run1++ )
02269         phi(1,run1) = pp*nablaY(1,run1);
02270 
02271     pp *= psi[number_][1]*scalT/h[0];
02272     for( run1 = 0; run1 < m; run1++ )
02273         phi(2,run1) = pp*nablaY(2,run1);
02274 
02275     pp *= psi[number_][2]*scalT/h[0];
02276     for( run1 = 0; run1 < m; run1++ )
02277         phi(3,run1) = pp*nablaY(3,run1);
02278 
02279     pp *= psi[number_][3]*scalT/h[0];
02280     for( run1 = 0; run1 < m; run1++ )
02281         phi(4,run1) = pp*nablaY(4,run1);
02282 
02283 
02284     for( run1 = 0; run1 < m; run1++ )
02285         delta(1,run1) = nablaY(0,run1) + phi(1,run1);
02286 
02287     for( run1 = 0; run1 < m; run1++ )
02288         delta(2,run1) = delta(1,run1) + phi(2,run1);
02289 
02290     for( run1 = 0; run1 < m; run1++ )
02291         delta(3,run1) = delta(2,run1) + phi(3,run1);
02292 
02293     for( run1 = 0; run1 < m; run1++ )
02294         eta[0][run1] = delta(3,run1) + phi(4,run1);
02295 
02296 
02297     for( run1 = 0; run1 < md; run1++ ){
02298 
02299         c2(run1) = - nablaY(0,run1)/psi[number_][0]
02300                    - delta (1,run1)/psi[number_][1]
02301                    - delta (2,run1)/psi[number_][2]
02302                    - delta (3,run1)/psi[number_][3];
02303     }
02304 
02305     returnValue returnvalue;
02306 
02307     correctorTime.start();
02308 
02309     returnvalue = determineCorrector(number_, ini);
02310 
02311     correctorTime.stop();
02312 
02313     if( returnvalue == RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF ){
02314         return -1;
02315     }
02316     if( returnvalue != SUCCESSFUL_RETURN ){
02317         return INFTY;
02318     }
02319 
02320     for( run1 = 0; run1 < m; run1++ ){
02321 
02322         nablaY_(0,run1) = eta[nOfNewtonSteps[number_]][run1];
02323 
02324         nablaY_(1,run1) = nablaY_(0,run1) - nablaY(0,run1);
02325         nablaY_(2,run1) = (nablaY_(1,run1) - nablaY(1,run1)*scalT)*(h[0]/psi[number_][1]);
02326         nablaY_(3,run1) = (nablaY_(2,run1) - nablaY(2,run1)*scalT*scalT)*(h[0]/psi[number_][2]);
02327         nablaY_(4,run1) = (nablaY_(3,run1) - nablaY(3,run1)*scalT*scalT*scalT)*(h[0]/psi[number_][3]);
02328     }
02329 
02330     double EE = EPS;
02331     for( run1 = 0; run1 < md; run1++ ){
02332         if( (eta[0][run1]-nablaY_(0,run1))/diff_scale(run1) >= EE  ){
02333             EE = (eta[0][run1]-nablaY_(0,run1))/diff_scale(run1);
02334         }
02335         if( (eta[0][run1]-nablaY_(0,run1))/diff_scale(run1) <= -EE ){
02336             EE = -(eta[0][run1]-nablaY_(0,run1))/diff_scale(run1);
02337         }
02338     }
02339 
02340     return 0.05*EE;
02341 }
02342 
02343 
02344 returnValue IntegratorBDF::determineCorrector( int stepnumber, BooleanType ini ){
02345 
02346     int run1, run2;
02347     int newtonsteps;
02348 
02349     double norm1 = 0.0;
02350     double norm2 = 0.0;
02351 
02352     BooleanType COMPUTE_JACOBIAN  = BT_FALSE;
02353     BooleanType JACOBIAN_COMPUTED = BT_FALSE;
02354 
02355     if( soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN ){
02356         COMPUTE_JACOBIAN = ini;
02357     }
02358 
02359     if( soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN ){
02360        nOfNewtonSteps[stepnumber] = 0;
02361     }
02362 
02363     if( stepnumber > 0 ){
02364         M_index[stepnumber] = M_index[stepnumber-1];
02365     }
02366 
02367     newtonsteps = 0;
02368 
02369     while( newtonsteps < 3 ){
02370 
02371        // evaluate F:
02372        // -----------
02373        x[time_index] = t + h[0];
02374        for( run2 = 0; run2 < m; run2++ ){
02375             x[diff_index[run2]] = eta[newtonsteps][run2];
02376        }
02377        for( run2 = 0; run2 < md; run2++ ){
02378            x[ddiff_index[run2]] = gamma[stepnumber][4]*eta[newtonsteps][run2]+c2(run2);
02379        }
02380 
02381        functionEvaluation.start();
02382 
02383        if( rhs[0].evaluate( 3*stepnumber+newtonsteps, x, F ) != SUCCESSFUL_RETURN ){
02384            return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02385        }
02386 
02387        relaxAlgebraic( F, x[time_index] );
02388 
02389        nFcnEvaluations++;
02390        functionEvaluation.stop();
02391 
02392        if( COMPUTE_JACOBIAN == BT_TRUE ){
02393 
02394            jacComputation.start();
02395 
02396            if( PrintLevel == HIGH ){
02397                cout  << "(RE-) COMPUTE-JACOBIAN... \n";
02398            }
02399 
02400            if( soa == SOA_FREEZING_MESH || soa == SOA_FREEZING_ALL ){
02401 
02402                if( nOfM >= maxNM ){
02403                    int oldMN = maxNM;
02404                    maxNM += maxNM;
02405                    M = (DMatrix**)realloc(M, maxNM*sizeof(DMatrix*));
02406                    for( run1 = oldMN; run1 < maxNM; run1++ )
02407                        M[run1] = 0;
02408                }
02409                M_index[stepnumber] = nOfM;
02410                M[nOfM] = new DMatrix(m,m);
02411                nOfM++;
02412            }
02413            else{
02414                if( M[0] == 0 ) M[0] = new DMatrix(m,m);
02415                M_index[stepnumber] = 0;
02416                M[0]->init(m,m);
02417            }
02418 
02419            for( run1 = 0; run1 < md; run1++ ){
02420                iseed[ddiff_index[run1]] = gamma[stepnumber][4];
02421                iseed[ diff_index[run1]] = 1.0;
02422                if( rhs[0].AD_forward( 3*stepnumber+newtonsteps, iseed,
02423                                       k2[0][0] ) != SUCCESSFUL_RETURN ){
02424                   return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02425                }
02426                for( run2 = 0; run2 < m; run2++ )
02427                    M[M_index[stepnumber]]->operator()(run2,run1) = k2[0][0][run2];
02428 
02429                iseed[ddiff_index[run1]] = 0.0;
02430                iseed[ diff_index[run1]] = 0.0;
02431            }
02432 
02433            for( run1 = 0; run1 < ma; run1++ ){
02434                iseed[ diff_index[md+run1]] = 1.0;
02435                if( rhs[0].AD_forward( 3*stepnumber+newtonsteps, iseed,
02436                                         k2[0][0] ) != SUCCESSFUL_RETURN ){
02437                   return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02438                }
02439                for( run2 = 0; run2 < m; run2++ )
02440                    M[M_index[stepnumber]]->operator()(run2,md+run1) = k2[0][0][run2];
02441 
02442                iseed[diff_index[md+run1]] = 0.0;
02443            }
02444 
02445            nJacEvaluations++;
02446            jacComputation.stop();
02447            jacDecomposition.start();
02448 
02449            if( decomposeJacobian(M_index[stepnumber], *M[M_index[stepnumber]] ) != SUCCESSFUL_RETURN )
02450                return ACADOERROR(RET_THE_DAE_INDEX_IS_TOO_LARGE);
02451 
02452            jacDecomposition.stop();
02453 
02454            JACOBIAN_COMPUTED = BT_TRUE;
02455            COMPUTE_JACOBIAN  = BT_FALSE;
02456        }
02457 
02458        norm1 = applyNewtonStep( M_index[stepnumber],
02459                                                            eta[newtonsteps+1],
02460                                 eta[newtonsteps],
02461                                *M[M_index[stepnumber]],
02462                                 F );
02463 
02464        if( soa == SOA_MESH_FROZEN || soa == SOA_EVERYTHING_FROZEN ){
02465            if( newtonsteps == nOfNewtonSteps[stepnumber] ){
02466                return SUCCESSFUL_RETURN;
02467            }
02468        }
02469 
02470 
02471        double subTOL = 1e-3*TOL;
02472 
02473        if( norm1 < subTOL ){
02474            nOfNewtonSteps[stepnumber] = newtonsteps+1;
02475            return SUCCESSFUL_RETURN;
02476        }
02477 
02478        if( newtonsteps == 0 ){
02479            norm2 = norm1;
02480        }
02481 
02482        if( newtonsteps == 1 && (norm1/norm2 > 0.33 || norm1/norm2 > sqrt(0.33*TOL/norm2)) ){
02483 
02484            if( JACOBIAN_COMPUTED == BT_FALSE ){
02485                COMPUTE_JACOBIAN = BT_TRUE;
02486                newtonsteps      = -1     ;
02487            }
02488        }
02489 
02490        if( newtonsteps == 1 ){
02491            norm2 = norm1;
02492        }
02493 
02494        if( newtonsteps == 2 ){
02495 
02496            if( norm1 > 0.33*TOL ){
02497 
02498                if( JACOBIAN_COMPUTED == BT_FALSE ){
02499                    COMPUTE_JACOBIAN = BT_TRUE;
02500                    newtonsteps      = -1     ;
02501                }
02502                else{
02503                    return RET_INFO_UNDEFINED;
02504                }
02505            }
02506            else{
02507                nOfNewtonSteps[stepnumber] = newtonsteps+1;
02508                return SUCCESSFUL_RETURN;
02509            }
02510        }
02511 
02512        newtonsteps++;
02513        nOfNewtonSteps[stepnumber] = newtonsteps;
02514     }
02515 
02516     return RET_INFO_UNDEFINED;
02517 }
02518 
02519 
02520 returnValue IntegratorBDF::rk_start(){
02521 
02522     int run1, run2, run3;
02523     double E = EPS;
02524 
02525      if( soa == SOA_EVERYTHING_FROZEN || soa == SOA_MESH_FROZEN ){
02526          rel_time_scale = h[0];
02527          h[0] = h[1];
02528      }
02529 
02530      // MEMORY REALLOCATION:
02531      // --------------------
02532      if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
02533 
02534          int oldAlloc = maxAlloc;
02535 
02536          if( maxAlloc < 8 ){
02537              maxAlloc = 8;
02538 
02539              psi            = (double**)realloc(psi,maxAlloc*sizeof(double*));
02540              gamma          = (double**)realloc(gamma,maxAlloc*sizeof(double*));
02541              nOfNewtonSteps = (int*)realloc(nOfNewtonSteps,(maxAlloc)*sizeof(int));
02542 
02543              for( run1 = oldAlloc; run1 < maxAlloc; run1++ ){
02544                  nOfNewtonSteps[run1] = 0;
02545              }
02546              for( run1 = oldAlloc; run1 < maxAlloc; run1++ ){
02547 
02548                  psi  [run1] = (double*)calloc(nstep,sizeof(double));
02549                  gamma[run1] = (double*)calloc(nstep,sizeof(double));
02550              }
02551          }
02552 
02553          M_index      = (int*)realloc(M_index,maxAlloc*sizeof(int));
02554          h = (double*)realloc(h,maxAlloc*sizeof(double));
02555      }
02556 
02557 
02558     if( soa != SOA_EVERYTHING_FROZEN ){
02559 
02560     returnValue returnvalue;
02561     int step_rejections = 0;
02562 
02563     while( step_rejections < 1 ){
02564 
02565     // determine k:
02566     // ------------
02567     run1 = 0;
02568 
02569     for( run2 = 0; run2 < m; run2++ ){
02570         k[0][run1][run2] = initial_guess[run2];
02571     }
02572 
02573     while( run1 < dim ){
02574 
02575         if( run1 > 0 ){
02576             for( run2 = 0; run2 < m; run2++ ){
02577                 k[0][run1][run2] = k[nOfNewtonSteps[run1-1]][run1-1][run2];
02578             }
02579         }
02580 
02581         returnvalue = rk_start_solve(run1);
02582 
02583         if( returnvalue != SUCCESSFUL_RETURN ){
02584 
02585             if( PrintLevel == HIGH ){
02586                 cout << "RUNGE-KUTTA STARTER: \n" << scientific
02587                          << "STEP REJECTED: error estimate           = " << E << endl
02588                          << "               required local tolerance = " << TOL << endl;
02589                 count3++;
02590             }
02591 
02592             if( soa == SOA_EVERYTHING_FROZEN || soa == SOA_MESH_FROZEN ){
02593                 return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02594             }
02595 
02596             if( h[0] <= hmin + EPS ){
02597                 return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02598             }
02599             h[0] = 0.2*h[0];
02600             if( h[0] < hmin ){
02601                 h[0] = hmin;
02602             }
02603             run1 = -1;
02604         }
02605 
02606         run1++;
02607     }
02608 
02609 
02610     // save previous eta4:
02611     // ----------------------------------------------
02612 
02613        for( run1 = 0; run1 < m; run1++ ){
02614            eta4_[run1]  = eta4[run1];
02615            eta5 [run1]  = eta4[run1];
02616        }
02617 
02618     // determine eta4 and eta5:
02619     // ----------------------------------------------
02620 
02621        for( run1 = 0; run1 < dim; run1++ ){
02622            for( run2 = 0; run2 < md; run2++ ){
02623                eta4[run2] = eta4[run2] + b4[run1]*h[0]*k[nOfNewtonSteps[run1]][run1][run2];
02624                eta5[run2] = eta5[run2] + b5[run1]*h[0]*k[nOfNewtonSteps[run1]][run1][run2];
02625            }
02626        }
02627        for( run1 = 0; run1 < ma; run1++ ){
02628            eta4[md+run1] = k[nOfNewtonSteps[dim-1]][dim-1][md+run1];
02629        }
02630 
02631     // determine local error estimate E:
02632     // ----------------------------------------------
02633 
02634        E = EPS;
02635        for( run2 = 0; run2 < md; run2++ ){
02636            if( (eta4[run2]-eta5[run2])/diff_scale(run2) >= E  ){
02637                E = (eta4[run2]-eta5[run2])/diff_scale(run2);
02638            }
02639            if( (eta4[run2]-eta5[run2])/diff_scale(run2) <= -E ){
02640                E = (-eta4[run2]+eta5[run2])/diff_scale(run2);
02641            }
02642        }
02643 
02644        if( E < TOL*h[0] || soa == SOA_EVERYTHING_FROZEN || soa == SOA_MESH_FROZEN  ){
02645           break;
02646        }
02647        else{
02648 
02649             if( PrintLevel == HIGH ){
02650                 cout << "RUNGE-KUTTA STARTER: \n" << scientific
02651                          << "STEP REJECTED: error estimate           = " << E << endl
02652                          << "               required local tolerance = " << TOL*h[0] << endl;
02653                 count3++;
02654             }
02655 
02656             for( run1 = 0; run1 < m; run1++ ){
02657                 eta4[run1] = eta4_[run1];
02658             }
02659             if( h[0] <= hmin + EPS ){
02660                 return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02661             }
02662             h[0] = 0.2*h[0];
02663             if( h[0] < hmin ){
02664                 h[0] = hmin;
02665             }
02666        }
02667 
02668        step_rejections++;
02669     }
02670 
02671     }
02672 
02673 
02674     // PROCEED IF THE STEP IS ACCEPTED:
02675     // --------------------------------
02676 
02677     // plan the first step of the BDF method:
02678     // --------------------------------------
02679 
02680     if( soa != SOA_EVERYTHING_FROZEN ){
02681 
02682     for( run1 = 0; run1 < nstep-1; run1++ ){
02683         for( run2 = 0; run2 < md; run2++ ){
02684             nablaY(nstep-2-run1,run2) = eta4_[run2];
02685             for( run3 = 0; run3 <= run1+3; run3++ ){
02686                 nablaY(nstep-2-run1,run2) = nablaY(nstep-2-run1,run2) +
02687                   h[0]*A[run1+3][run3]*k[nOfNewtonSteps[run3]][run3][run2];
02688             }
02689         }
02690         for( run2 = 0; run2 < ma; run2++ ){
02691             nablaY(nstep-2-run1,md+run2) = k[nOfNewtonSteps[run1+3]][run1+3][md+run2];
02692         }
02693     }
02694 
02695     }
02696 
02697 
02698      // compute forward derivatives if requested:
02699      // ------------------------------------------
02700 
02701      if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 ){
02702 
02703          if( nBDirs != 0 ){
02704              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
02705          }
02706 
02707          determineRKEtaGForward();
02708      }
02709 
02710      if( nBDirs > 0 ){
02711 
02712          if( soa != SOA_EVERYTHING_FROZEN ){
02713              return ACADOERROR(RET_NOT_FROZEN);
02714          }
02715          if( nFDirs != 0 || nBDirs2 != 0 || nFDirs2 != 0 ){
02716              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
02717          }
02718          determineRKEtaHBackward();
02719      }
02720      if( nFDirs2 > 0 ){
02721 
02722          if( soa != SOA_EVERYTHING_FROZEN ){
02723              return ACADOERROR(RET_NOT_FROZEN);
02724          }
02725          if( nBDirs != 0 || nBDirs2 != 0 || nFDirs != 1 ){
02726              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
02727          }
02728          determineRKEtaGForward2();
02729      }
02730      if( nBDirs2 > 0 ){
02731 
02732          if( soa != SOA_EVERYTHING_FROZEN ){
02733              return ACADOERROR(RET_NOT_FROZEN);
02734          }
02735          if( nBDirs != 0 || nFDirs2 != 0 || nFDirs != 1 ){
02736              return ACADOERROR(RET_WRONG_DEFINITION_OF_SEEDS);
02737          }
02738 
02739          determineRKEtaHBackward2();
02740      }
02741 
02742     // Printing:
02743     // ---------
02744     printRKIntermediateResults();
02745 
02746     const double hhh = c[3]*h[0];
02747 
02748 
02749     // STORAGE:
02750     // --------
02751     if( soa == SOA_FREEZING_MESH || soa == SOA_FREEZING_ALL ){
02752 
02753         h[1] = h[0];
02754         for( run3 = 3; run3 >= 0; run3-- )
02755             h[2+run3] = hhh;
02756     }
02757 
02758 
02759     int i1 = timeInterval.getFloorIndex( t             );
02760     int i2 = timeInterval.getFloorIndex( t + c[6]*h[0] );
02761     int jj;
02762 
02763     for( jj = i1+1; jj <= i2; jj++ ){
02764 
02765         for( run1 = 0; run1 < m; run1++ ){
02766             if( nFDirs == 0 && nBDirs  == 0 && nFDirs2 == 0 && nBDirs == 0 )   xStore( jj, run1 ) = nablaY (3,run1);
02767             if( nFDirs  > 0 && nBDirs2 == 0 && nFDirs2 == 0                )  dxStore( jj, run1 ) = nablaG (3,run1);
02768             if( nFDirs2 > 0                                                ) ddxStore( jj, run1 ) = nablaG3(3,run1);
02769         }
02770         for( run1 = 0; run1 < mn; run1++ )
02771             iStore( jj, run1 ) = x[rhs->index( VT_INTERMEDIATE_STATE, run1 )];
02772     }
02773 
02774 
02775     // PREPARE MODIFIED DIVIDED DIFFERENCES:
02776     // -------------------------------------
02777     if( soa != SOA_EVERYTHING_FROZEN )
02778         prepareDividedDifferences( nablaY );
02779 
02780     if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 )
02781         prepareDividedDifferences( nablaG );
02782 
02783     if( nFDirs2 > 0 ){
02784         prepareDividedDifferences( nablaG2 );
02785         prepareDividedDifferences( nablaG3 );
02786     }
02787 
02788 
02789     if( soa != SOA_EVERYTHING_FROZEN ){
02790 
02791         if( soa != SOA_FREEZING_ALL ){
02792 
02793             psi[7][0] = hhh    ;
02794             psi[7][1] = 2.0*hhh;
02795             psi[7][2] = 3.0*hhh;
02796             psi[7][3] = 4.0*hhh;
02797         }
02798         else{
02799 
02800             psi[6][0] = hhh    ;
02801             psi[6][1] = 2.0*hhh;
02802             psi[6][2] = 3.0*hhh;
02803             psi[6][3] = 4.0*hhh;
02804         }
02805     }
02806 
02807     // increase the time:
02808     // ----------------------------------------------
02809 
02810     t = t + c[6]*h[0];
02811 
02812     if( soa != SOA_EVERYTHING_FROZEN && soa != SOA_MESH_FROZEN ){
02813 
02814     // determine the new step size:
02815     // ----------------------------------------------
02816 
02817         rel_time_scale = 0.2*h[0];
02818         h[0] = 0.2*h[0]*pow( tune*(TOL*h[0]/E) , 1.0/3.0 );
02819 
02820         if( h[0] > hmax ){
02821             h[0] = hmax;
02822         }
02823         if( h[0] < hmin ){
02824             h[0] = hmin;
02825         }
02826 
02827         if( t + h[0] >= timeInterval.getLastTime() ){
02828             h[0] = timeInterval.getLastTime()-t;
02829         }
02830     }
02831     else{
02832         h[0] = 0.2*h[0];
02833     }
02834 
02835     return SUCCESSFUL_RETURN;
02836 }
02837 
02838 
02839 returnValue IntegratorBDF::rk_start_solve( int stepnumber ){
02840 
02841     int run1, run2, run3;
02842     int newtonsteps;
02843 
02844     double norm1 = 0.0;
02845     double norm2 = 0.0;
02846 
02847     BooleanType COMPUTE_JACOBIAN;
02848     BooleanType JACOBIAN_COMPUTED = BT_FALSE;
02849 
02850     if( stepnumber == 0 && soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN ){
02851         COMPUTE_JACOBIAN = BT_TRUE;
02852     }
02853     else{
02854         COMPUTE_JACOBIAN = BT_FALSE;
02855     }
02856 
02857     if( soa != SOA_MESH_FROZEN && soa != SOA_EVERYTHING_FROZEN ){
02858        nOfNewtonSteps[stepnumber] = 0;
02859     }
02860 
02861     if( stepnumber > 0 ){
02862         M_index[stepnumber] = M_index[stepnumber-1];
02863     }
02864 
02865     newtonsteps = 0;
02866 
02867     while( newtonsteps < 3 ){
02868 
02869        // evaluate F:
02870        // -----------
02871        x[time_index] = t + c[stepnumber]*h[0];
02872        for( run2 = 0; run2 < md; run2++ ){
02873             x[diff_index[run2]] = eta4[run2];
02874             for( run3 = 0; run3 < stepnumber+1; run3++ ){
02875                  x[diff_index[run2]] = x[diff_index[run2]] +
02876                           A[stepnumber][run3]*h[0]*k[nOfNewtonSteps[run3]][run3][run2];
02877             }
02878        }
02879        for( run2 = 0; run2 < ma; run2++ ){
02880            x[diff_index[md+run2]] = k[newtonsteps][stepnumber][md+run2];
02881        }
02882        for( run2 = 0; run2 < md; run2++ ){
02883            x[ddiff_index[run2]] = k[newtonsteps][stepnumber][run2];
02884        }
02885 
02886        functionEvaluation.start();
02887 
02888        if( rhs[0].evaluate( 3*stepnumber+newtonsteps, x, F ) != SUCCESSFUL_RETURN ){
02889            return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02890        }
02891 
02892        relaxAlgebraic( F, x[time_index] );
02893 
02894        nFcnEvaluations++;
02895        functionEvaluation.stop();
02896 
02897 
02898        if( COMPUTE_JACOBIAN == BT_TRUE ){
02899 
02900            jacComputation.start();
02901 
02902            if( PrintLevel == HIGH ){
02903                cout << "(RE-) COMPUTE-JACOBIAN... \n";
02904            }
02905 
02906            const double ise = h[0]*A[stepnumber][stepnumber];
02907 
02908            if( soa == SOA_FREEZING_MESH || soa == SOA_FREEZING_ALL ){
02909 
02910                if( nOfM >= maxNM ){
02911                    int oldMN = maxNM;
02912                    maxNM += maxNM;
02913                    M = (DMatrix**)realloc(M, maxNM*sizeof(DMatrix*));
02914                    for( run1 = oldMN; run1 < maxNM; run1++ )
02915                        M[run1] = 0;
02916                }
02917                M_index[stepnumber] = nOfM;
02918                M[nOfM] = new DMatrix(m,m);
02919                nOfM++;
02920            }
02921            else{
02922                if( M[0] == 0 ) M[0] = new DMatrix(m,m);
02923                M_index[stepnumber] = 0;
02924                M[0]->init(m,m);
02925            }
02926 
02927            for( run1 = 0; run1 < md; run1++ ){
02928                iseed[ddiff_index[run1]] = 1.0;
02929                iseed[ diff_index[run1]] = ise;
02930                if( rhs[0].AD_forward( 3*stepnumber+newtonsteps, iseed,
02931                                       k2[0][0] ) != SUCCESSFUL_RETURN ){
02932                   return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02933                }
02934 
02935                for( run2 = 0; run2 < m; run2++ )
02936                    M[M_index[stepnumber]]->operator()(run2,run1) = k2[0][0][run2];
02937 
02938                iseed[ddiff_index[run1]] = 0.0;
02939                iseed[ diff_index[run1]] = 0.0;
02940            }
02941            for( run1 = 0; run1 < ma; run1++ ){
02942                iseed[ diff_index[md+run1]] = 1.0;
02943                if( rhs[0].AD_forward( 3*stepnumber+newtonsteps, iseed,
02944                                       k2[0][0] ) != SUCCESSFUL_RETURN ){
02945                   return ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
02946                }
02947 
02948                for( run2 = 0; run2 < m; run2++ )
02949                    M[M_index[stepnumber]]->operator()(run2,md+run1) = k2[0][0][run2];
02950 
02951                iseed[diff_index[md+run1]] = 0.0;
02952            }
02953 
02954            nJacEvaluations++;
02955            jacComputation.stop();
02956            jacDecomposition.start();
02957 
02958            if( decomposeJacobian(M_index[stepnumber], *M[M_index[stepnumber]] ) != SUCCESSFUL_RETURN )
02959                    return ACADOERROR(RET_THE_DAE_INDEX_IS_TOO_LARGE);
02960 
02961            jacDecomposition.stop();
02962 
02963            JACOBIAN_COMPUTED = BT_TRUE;
02964            COMPUTE_JACOBIAN  = BT_FALSE;
02965        }
02966 
02967        norm1 = applyNewtonStep( M_index[stepnumber],
02968                                                    k[newtonsteps+1][stepnumber],
02969                                 k[newtonsteps][stepnumber]  ,
02970                                *M[M_index[stepnumber]],
02971                                 F );
02972 
02973        if( soa == SOA_MESH_FROZEN || soa == SOA_EVERYTHING_FROZEN ){
02974            if( newtonsteps == nOfNewtonSteps[stepnumber] ){
02975                return SUCCESSFUL_RETURN;
02976            }
02977        }
02978 
02979        double subTOL = 1e-3*TOL;
02980 
02981        if( norm1 < subTOL ){
02982 
02983            nOfNewtonSteps[stepnumber] = newtonsteps+1;
02984            return SUCCESSFUL_RETURN;
02985        }
02986 
02987        if( newtonsteps == 0 ){
02988            norm2 = norm1;
02989        }
02990 
02991        if( newtonsteps == 1 && (norm1/norm2 > 0.33 || norm1/norm2 > sqrt(0.33*TOL/norm2)) ){
02992 
02993            if( JACOBIAN_COMPUTED == BT_FALSE ){
02994                COMPUTE_JACOBIAN = BT_TRUE;
02995                newtonsteps      = -1     ;
02996            }
02997        }
02998 
02999        if( newtonsteps == 1 ){
03000            norm2 = norm1;
03001        }
03002 
03003        if( newtonsteps == 2 ){
03004 
03005            if( JACOBIAN_COMPUTED == BT_FALSE ){
03006                COMPUTE_JACOBIAN = BT_TRUE;
03007                newtonsteps      = -1     ;
03008            }
03009            else{
03010 
03011                if( norm1 > 0.33*TOL ){
03012                    // printf("RK - REJECTION !!!!!!  norm1 = %.16e \n", norm1 );
03013                    return RET_INFO_UNDEFINED;
03014                }
03015 
03016                nOfNewtonSteps[stepnumber] = newtonsteps+1;
03017                return SUCCESSFUL_RETURN;
03018            }
03019        }
03020 
03021        newtonsteps++;
03022        nOfNewtonSteps[stepnumber] = newtonsteps;
03023     }
03024 
03025     return RET_INFO_UNDEFINED;
03026 }
03027 
03028 
03029 void IntegratorBDF::determineRKEtaGForward(){
03030 
03031     int run1, run2, run3, newtonsteps;
03032 
03033 
03034     for( run2 = 0; run2 < m; run2++ ){
03035         k[0][0][run2] = 0.0;
03036     }
03037 
03038 
03039     // determine k:
03040     // -----------------------------------------------
03041        for( run1 = 0; run1 < dim; run1++ ){
03042 
03043            if( run1 > 0 ){
03044                for( run2 = 0; run2 < m; run2++ ){
03045                    k[0][run1][run2] = k[nOfNewtonSteps[run1-1]][run1-1][run2];
03046                }
03047            }
03048 
03049            newtonsteps = 0;
03050            while( newtonsteps < nOfNewtonSteps[run1] ){
03051 
03052                for( run2 = 0; run2 < md; run2++ ){
03053                    G[diff_index[run2]] = etaG[run2];
03054                    for( run3 = 0; run3 < run1; run3++ ){
03055                        G[diff_index[run2]] = G[diff_index[run2]] +
03056                                    A[run1][run3]*h[0]*k[nOfNewtonSteps[run3]][run3][run2];
03057                    }
03058                    G[diff_index[run2]] = G[diff_index[run2]] +
03059                                A[run1][run1]*h[0]*k[newtonsteps][run1][run2];
03060                }
03061                for( run2 = 0; run2 < ma; run2++ ){
03062                    G[diff_index[md+run2]] = k[newtonsteps][run1][md+run2];
03063                }
03064                for( run2 = 0; run2 < md; run2++ ){
03065                    G[ddiff_index[run2]] = k[newtonsteps][run1][run2];
03066                }
03067 
03068                if( rhs[0].AD_forward( 3*run1+newtonsteps, G, F )
03069                                       != SUCCESSFUL_RETURN ){
03070                    ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03071                    return;
03072                }
03073 
03074 
03075                applyNewtonStep( M_index[run1],
03076                                                    k[newtonsteps+1][run1],
03077                                 k[newtonsteps][run1]  ,
03078                                *M[M_index[run1]],
03079                                 F );
03080 
03081                newtonsteps++;
03082            }
03083        }
03084 
03085     // determine nablaG:
03086     // ----------------------------------------------
03087 
03088        for( run1 = 0; run1 < nstep-1; run1++ ){
03089            for( run2 = 0; run2 < md; run2++ ){
03090                nablaG(nstep-2-run1,run2) = etaG[run2];
03091                for( run3 = 0; run3 <= run1+3; run3++ ){
03092                    nablaG(nstep-2-run1,run2) = nablaG(nstep-2-run1,run2) +
03093                      h[0]*A[run1+3][run3]*k[nOfNewtonSteps[run3]][run3][run2];
03094                }
03095            }
03096            for( run2 = 0; run2 < ma; run2++ ){
03097                nablaG(nstep-2-run1,md+run2) = k[nOfNewtonSteps[run1+3]][run1+3][md+run2];
03098            }
03099         }
03100 }
03101 
03102 
03103 void IntegratorBDF::determineRKEtaHBackward(){
03104 
03105     int run1, run2, run3, newtonsteps;
03106     const double hhh = c[3]*h[0];
03107     int number_;
03108     double *Hstore = new double[ndir];
03109 
03110     const double scalT = rel_time_scale/hhh;
03111 
03112         for( run3 = 0; run3 < 4; run3++ ){
03113             for( run1 = 0; run1 < dim; run1++ ){
03114                 for( run2 = 0; run2 < ndir; run2++ ){
03115                     l[run3][run1][run2] = 0.0;
03116                 }
03117             }
03118         }
03119 
03120         for( run1 = 0; run1 < ndir; run1++ ){
03121             Hstore[run1] = nablaH(0,run1);
03122         }
03123 
03124         for( run1 = 0; run1 < ndir; run1++ ){
03125 
03126             zH[5][run1] = nablaH(3,run1)*scalT*scalT;
03127             zH[4][run1] = nablaH(2,run1)*scalT + zH[5][run1]/(3.0);
03128             zH[3][run1] = - zH[5][run1]/(3.0);
03129             zH[2][run1] = nablaH(1,run1) + zH[4][run1]/(2.0);
03130             zH[1][run1] = (zH[3][run1]-zH[4][run1])/(2.0);
03131             zH[0][run1] = - zH[3][run1]/(2.0);
03132 
03133             nablaH(3,run1) = nablaH(0,run1)*h[0] + zH[2][run1]*(h[0]/hhh);
03134             nablaH(2,run1) = (zH[1][run1]-zH[2][run1])*(h[0]/hhh);
03135             nablaH(1,run1) = (zH[0][run1]-zH[1][run1])*(h[0]/hhh);
03136             nablaH(0,run1) = -zH[0][run1]*(h[0]/hhh);
03137         }
03138 
03139         for( run1 = 0; run1 < ndir; run1++ ){
03140             kH[6][nOfNewtonSteps[6]][run1] = nablaH(3,run1)/h[0];
03141         }
03142         for( run1 = 0; run1 < md; run1++ ){
03143             kH[6][nOfNewtonSteps[6]][diff_index[run1]] =
03144             kH[6][nOfNewtonSteps[6]][diff_index[run1]]*h[0]*A[6][6];
03145         }
03146 
03147         newtonsteps = nOfNewtonSteps[6];
03148         newtonsteps--;
03149 
03150         number_ = 6;
03151 
03152         while( newtonsteps >= 0 ){
03153 
03154             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03155                             *M[M_index[number_]],
03156                              H );
03157 
03158             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03159                 != SUCCESSFUL_RETURN ){
03160                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03161             }
03162 
03163 
03164             for( run1 = 0; run1 < ndir; run1++ ){
03165                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03166                                                 - l[newtonsteps][number_][run1];
03167             }
03168             for( run1 = 0; run1 < md; run1++ ){
03169                 kH[number_][newtonsteps][diff_index[run1]] = 
03170                    kH[number_][newtonsteps+1][diff_index[run1]]
03171                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03172                  - l[newtonsteps][number_][ddiff_index[run1]];
03173             }
03174 
03175             newtonsteps--;
03176         }
03177 
03178 
03179         for( run1 = 0; run1 < ndir; run1++ ){
03180             kH[5][nOfNewtonSteps[5]][run1] = kH[6][0][run1] + nablaH(2,run1)/h[0];
03181         }
03182         for( run1 = 0; run1 < md; run1++ ){
03183             kH[5][nOfNewtonSteps[5]][diff_index[run1]] = kH[6][0][diff_index[run1]]
03184             +  nablaH(3,diff_index[run1])*A[6][5]
03185             +  nablaH(2,diff_index[run1])*A[5][5];
03186             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03187                kH[5][nOfNewtonSteps[5]][diff_index[run1]] =
03188                kH[5][nOfNewtonSteps[5]][diff_index[run1]]
03189                -l[run3][number_][diff_index[run1]]*h[0]*A[6][5];
03190             }
03191         }
03192 
03193 
03194         newtonsteps = nOfNewtonSteps[5];
03195         newtonsteps--;
03196 
03197         number_ = 5;
03198 
03199         while( newtonsteps >= 0 ){
03200 
03201             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03202                             *M[M_index[number_]],
03203                              H );
03204 
03205             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03206                 != SUCCESSFUL_RETURN ){
03207                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03208             }
03209 
03210 
03211             for( run1 = 0; run1 < ndir; run1++ ){
03212                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03213                                                 - l[newtonsteps][number_][run1];
03214             }
03215             for( run1 = 0; run1 < md; run1++ ){
03216                 kH[number_][newtonsteps][diff_index[run1]] = 
03217                    kH[number_][newtonsteps+1][diff_index[run1]]
03218                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03219                  - l[newtonsteps][number_][ddiff_index[run1]];
03220             }
03221 
03222             newtonsteps--;
03223         }
03224 
03225 
03226         for( run1 = 0; run1 < ndir; run1++ ){
03227             kH[4][nOfNewtonSteps[4]][run1] = kH[5][0][run1] + nablaH(1,run1)/h[0];
03228         }
03229         for( run1 = 0; run1 < md; run1++ ){
03230             kH[4][nOfNewtonSteps[4]][diff_index[run1]] = kH[5][0][diff_index[run1]]
03231             +  nablaH(3,diff_index[run1])*A[6][4]
03232             +  nablaH(2,diff_index[run1])*A[5][4]
03233             +  nablaH(1,diff_index[run1])*A[4][4];
03234             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03235                kH[4][nOfNewtonSteps[4]][diff_index[run1]] =
03236                kH[4][nOfNewtonSteps[4]][diff_index[run1]]
03237                -l[run3][5][diff_index[run1]]*h[0]*A[5][4];
03238             }
03239             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03240                kH[4][nOfNewtonSteps[4]][diff_index[run1]] =
03241                kH[4][nOfNewtonSteps[4]][diff_index[run1]]
03242                -l[run3][6][diff_index[run1]]*h[0]*A[6][4];
03243             }
03244         }
03245 
03246 
03247         newtonsteps = nOfNewtonSteps[4];
03248         newtonsteps--;
03249 
03250         number_ = 4;
03251 
03252         while( newtonsteps >= 0 ){
03253 
03254             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03255                             *M[M_index[number_]],
03256                              H );
03257 
03258             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03259                 != SUCCESSFUL_RETURN ){
03260                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03261             }
03262 
03263 
03264             for( run1 = 0; run1 < ndir; run1++ ){
03265                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03266                                                 - l[newtonsteps][number_][run1];
03267             }
03268             for( run1 = 0; run1 < md; run1++ ){
03269                 kH[number_][newtonsteps][diff_index[run1]] = 
03270                    kH[number_][newtonsteps+1][diff_index[run1]]
03271                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03272                  - l[newtonsteps][number_][ddiff_index[run1]];
03273             }
03274 
03275             newtonsteps--;
03276         }
03277 
03278 
03279         for( run1 = 0; run1 < ndir; run1++ ){
03280             kH[3][nOfNewtonSteps[3]][run1] = kH[4][0][run1] + nablaH(0,run1)/h[0];
03281         }
03282         for( run1 = 0; run1 < md; run1++ ){
03283             kH[3][nOfNewtonSteps[3]][diff_index[run1]] = kH[4][0][diff_index[run1]]
03284             +  nablaH(3,diff_index[run1])*A[6][3]
03285             +  nablaH(2,diff_index[run1])*A[5][3]
03286             +  nablaH(1,diff_index[run1])*A[4][3]
03287             +  nablaH(0,diff_index[run1])*A[3][3];
03288             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03289                kH[3][nOfNewtonSteps[3]][diff_index[run1]] =
03290                kH[3][nOfNewtonSteps[3]][diff_index[run1]]
03291                -l[run3][4][diff_index[run1]]*h[0]*A[4][3];
03292             }
03293             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03294                kH[3][nOfNewtonSteps[3]][diff_index[run1]] =
03295                kH[3][nOfNewtonSteps[3]][diff_index[run1]]
03296                -l[run3][5][diff_index[run1]]*h[0]*A[5][3];
03297             }
03298             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03299                kH[3][nOfNewtonSteps[3]][diff_index[run1]] =
03300                kH[3][nOfNewtonSteps[3]][diff_index[run1]]
03301                -l[run3][6][diff_index[run1]]*h[0]*A[6][3];
03302             }
03303         }
03304 
03305 
03306         newtonsteps = nOfNewtonSteps[3];
03307         newtonsteps--;
03308 
03309         number_ = 3;
03310 
03311         while( newtonsteps >= 0 ){
03312 
03313             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03314                             *M[M_index[number_]],
03315                              H );
03316 
03317             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03318                 != SUCCESSFUL_RETURN ){
03319                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03320             }
03321 
03322 
03323             for( run1 = 0; run1 < ndir; run1++ ){
03324                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03325                                                 - l[newtonsteps][number_][run1];
03326             }
03327             for( run1 = 0; run1 < md; run1++ ){
03328                 kH[number_][newtonsteps][diff_index[run1]] = 
03329                    kH[number_][newtonsteps+1][diff_index[run1]]
03330                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03331                  - l[newtonsteps][number_][ddiff_index[run1]];
03332             }
03333 
03334             newtonsteps--;
03335         }
03336 
03337 
03338         for( run1 = 0; run1 < ndir; run1++ ){
03339             kH[2][nOfNewtonSteps[2]][run1] = kH[3][0][run1];
03340         }
03341         for( run1 = 0; run1 < md; run1++ ){
03342             kH[2][nOfNewtonSteps[2]][diff_index[run1]] = kH[3][0][diff_index[run1]]
03343             +  nablaH(3,diff_index[run1])*A[6][2]
03344             +  nablaH(2,diff_index[run1])*A[5][2]
03345             +  nablaH(1,diff_index[run1])*A[4][2]
03346             +  nablaH(0,diff_index[run1])*A[3][2];
03347             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
03348                kH[2][nOfNewtonSteps[2]][diff_index[run1]] =
03349                kH[2][nOfNewtonSteps[2]][diff_index[run1]]
03350                -l[run3][3][diff_index[run1]]*h[0]*A[3][2];
03351             }
03352             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03353                kH[2][nOfNewtonSteps[2]][diff_index[run1]] =
03354                kH[2][nOfNewtonSteps[2]][diff_index[run1]]
03355                -l[run3][4][diff_index[run1]]*h[0]*A[4][2];
03356             }
03357             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03358                kH[2][nOfNewtonSteps[2]][diff_index[run1]] =
03359                kH[2][nOfNewtonSteps[2]][diff_index[run1]]
03360                -l[run3][5][diff_index[run1]]*h[0]*A[5][2];
03361             }
03362             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03363                kH[2][nOfNewtonSteps[2]][diff_index[run1]] =
03364                kH[2][nOfNewtonSteps[2]][diff_index[run1]]
03365                -l[run3][6][diff_index[run1]]*h[0]*A[6][2];
03366             }
03367         }
03368 
03369 
03370         newtonsteps = nOfNewtonSteps[2];
03371         newtonsteps--;
03372 
03373         number_ = 2;
03374 
03375         while( newtonsteps >= 0 ){
03376 
03377             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03378                             *M[M_index[number_]],
03379                              H );
03380 
03381             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03382                 != SUCCESSFUL_RETURN ){
03383                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03384             }
03385 
03386 
03387             for( run1 = 0; run1 < ndir; run1++ ){
03388                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03389                                                 - l[newtonsteps][number_][run1];
03390             }
03391             for( run1 = 0; run1 < md; run1++ ){
03392                 kH[number_][newtonsteps][diff_index[run1]] = 
03393                    kH[number_][newtonsteps+1][diff_index[run1]]
03394                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03395                  - l[newtonsteps][number_][ddiff_index[run1]];
03396             }
03397 
03398             newtonsteps--;
03399         }
03400 
03401 
03402         for( run1 = 0; run1 < ndir; run1++ ){
03403             kH[1][nOfNewtonSteps[1]][run1] = kH[2][0][run1];
03404         }
03405         for( run1 = 0; run1 < md; run1++ ){
03406             kH[1][nOfNewtonSteps[1]][diff_index[run1]] = kH[2][0][diff_index[run1]]
03407             +  nablaH(3,diff_index[run1])*A[6][1]
03408             +  nablaH(2,diff_index[run1])*A[5][1]
03409             +  nablaH(1,diff_index[run1])*A[4][1]
03410             +  nablaH(0,diff_index[run1])*A[3][1];
03411             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
03412                kH[1][nOfNewtonSteps[1]][diff_index[run1]] =
03413                kH[1][nOfNewtonSteps[1]][diff_index[run1]]
03414                -l[run3][2][diff_index[run1]]*h[0]*A[2][1];
03415             }
03416             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
03417                kH[1][nOfNewtonSteps[1]][diff_index[run1]] =
03418                kH[1][nOfNewtonSteps[1]][diff_index[run1]]
03419                -l[run3][3][diff_index[run1]]*h[0]*A[3][1];
03420             }
03421             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03422                kH[1][nOfNewtonSteps[1]][diff_index[run1]] =
03423                kH[1][nOfNewtonSteps[1]][diff_index[run1]]
03424                -l[run3][4][diff_index[run1]]*h[0]*A[4][1];
03425             }
03426             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03427                kH[1][nOfNewtonSteps[1]][diff_index[run1]] =
03428                kH[1][nOfNewtonSteps[1]][diff_index[run1]]
03429                -l[run3][5][diff_index[run1]]*h[0]*A[5][1];
03430             }
03431             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03432                kH[1][nOfNewtonSteps[1]][diff_index[run1]] =
03433                kH[1][nOfNewtonSteps[1]][diff_index[run1]]
03434                -l[run3][6][diff_index[run1]]*h[0]*A[6][1];
03435             }
03436         }
03437 
03438 
03439         newtonsteps = nOfNewtonSteps[1];
03440         newtonsteps--;
03441 
03442         number_ = 1;
03443 
03444         while( newtonsteps >= 0 ){
03445 
03446             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03447                             *M[M_index[number_]],
03448                              H );
03449 
03450             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03451                 != SUCCESSFUL_RETURN ){
03452                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03453             }
03454 
03455 
03456             for( run1 = 0; run1 < ndir; run1++ ){
03457                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03458                                                 - l[newtonsteps][number_][run1];
03459             }
03460             for( run1 = 0; run1 < md; run1++ ){
03461                 kH[number_][newtonsteps][diff_index[run1]] = 
03462                    kH[number_][newtonsteps+1][diff_index[run1]]
03463                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03464                  - l[newtonsteps][number_][ddiff_index[run1]];
03465             }
03466 
03467             newtonsteps--;
03468         }
03469 
03470 
03471 
03472         for( run1 = 0; run1 < ndir; run1++ ){
03473             kH[0][nOfNewtonSteps[0]][run1] = kH[1][0][run1];
03474         }
03475         for( run1 = 0; run1 < md; run1++ ){
03476             kH[0][nOfNewtonSteps[0]][diff_index[run1]] = kH[1][0][diff_index[run1]]
03477             +  nablaH(3,diff_index[run1])*A[6][0]
03478             +  nablaH(2,diff_index[run1])*A[5][0]
03479             +  nablaH(1,diff_index[run1])*A[4][0]
03480             +  nablaH(0,diff_index[run1])*A[3][0];
03481             for( run3 = 0; run3 < nOfNewtonSteps[1]; run3++){
03482                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03483                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03484                -l[run3][1][diff_index[run1]]*h[0]*A[1][0];
03485             }
03486             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
03487                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03488                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03489                -l[run3][2][diff_index[run1]]*h[0]*A[2][0];
03490             }
03491             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
03492                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03493                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03494                -l[run3][3][diff_index[run1]]*h[0]*A[3][0];
03495             }
03496             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03497                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03498                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03499                -l[run3][4][diff_index[run1]]*h[0]*A[4][0];
03500             }
03501             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03502                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03503                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03504                -l[run3][5][diff_index[run1]]*h[0]*A[5][0];
03505             }
03506             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03507                kH[0][nOfNewtonSteps[0]][diff_index[run1]] =
03508                kH[0][nOfNewtonSteps[0]][diff_index[run1]]
03509                -l[run3][6][diff_index[run1]]*h[0]*A[6][0];
03510             }
03511         }
03512 
03513         newtonsteps = nOfNewtonSteps[0];
03514         newtonsteps--;
03515 
03516         number_ = 0;
03517 
03518         while( newtonsteps >= 0 ){
03519 
03520             applyMTranspose( M_index[number_], kH[number_][newtonsteps+1],
03521                             *M[M_index[number_]],
03522                              H );
03523 
03524             if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][number_] )
03525                 != SUCCESSFUL_RETURN ){
03526                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03527             }
03528 
03529             for( run1 = 0; run1 < ndir; run1++ ){
03530                 kH[number_][newtonsteps][run1] =   kH[number_][newtonsteps+1][run1]
03531                                                 - l[newtonsteps][number_][run1];
03532             }
03533             for( run1 = 0; run1 < md; run1++ ){
03534                 kH[number_][newtonsteps][diff_index[run1]] = 
03535                    kH[number_][newtonsteps+1][diff_index[run1]]
03536                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03537                  - l[newtonsteps][number_][ddiff_index[run1]];
03538             }
03539 
03540             newtonsteps--;
03541         }
03542 
03543         for( run1 = 0; run1 < ndir; run1++ ){
03544             nablaH(0,run1) = Hstore[run1];
03545         }
03546         for( run1 = 0; run1 < ndir; run1++ ){
03547             for( run3 = 0; run3 < nOfNewtonSteps[0]; run3++)
03548                nablaH(0,run1) -= l[run3][0][run1];
03549             for( run3 = 0; run3 < nOfNewtonSteps[1]; run3++)
03550                nablaH(0,run1) -= l[run3][1][run1];
03551             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++)
03552                nablaH(0,run1) -= l[run3][2][run1];
03553             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++)
03554                nablaH(0,run1) -= l[run3][3][run1];
03555             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++)
03556                nablaH(0,run1) -=  l[run3][4][run1];
03557             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++)
03558                nablaH(0,run1) -= l[run3][5][run1];
03559             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++)
03560                nablaH(0,run1) -= l[run3][6][run1];
03561         }
03562 
03563     delete[] Hstore;
03564 }
03565 
03566 
03567 void IntegratorBDF::determineRKEtaGForward2(){
03568 
03569     int run1, run2, run3, newtonsteps;
03570 
03571 
03572     for( run2 = 0; run2 < m; run2++ ){
03573         k [0][0][run2] = 0.0;
03574         k2[0][0][run2] = 0.0;
03575     }
03576 
03577 
03578     // determine k:
03579     // -----------------------------------------------
03580        for( run1 = 0; run1 < dim; run1++ ){
03581 
03582            if( run1 > 0 ){
03583                for( run2 = 0; run2 < m; run2++ ){
03584                    k [0][run1][run2] = k [nOfNewtonSteps[run1-1]][run1-1][run2];
03585                    k2[0][run1][run2] = k2[nOfNewtonSteps[run1-1]][run1-1][run2];
03586                }
03587            }
03588 
03589            newtonsteps = 0;
03590            while( newtonsteps < nOfNewtonSteps[run1] ){
03591 
03592                for( run2 = 0; run2 < md; run2++ ){
03593                    G2[diff_index[run2]] = etaG2[run2];
03594                    G3[diff_index[run2]] = etaG3[run2];
03595                    for( run3 = 0; run3 < run1; run3++ ){
03596                        G2[diff_index[run2]] = G2[diff_index[run2]] +
03597                                    A[run1][run3]*h[0]*k [nOfNewtonSteps[run3]][run3][run2];
03598                        G3[diff_index[run2]] = G3[diff_index[run2]] +
03599                                    A[run1][run3]*h[0]*k2[nOfNewtonSteps[run3]][run3][run2];
03600                    }
03601                    G2[diff_index[run2]] = G2[diff_index[run2]] +
03602                                A[run1][run1]*h[0]*k [newtonsteps][run1][run2];
03603                    G3[diff_index[run2]] = G3[diff_index[run2]] +
03604                                A[run1][run1]*h[0]*k2[newtonsteps][run1][run2];
03605                }
03606                for( run2 = 0; run2 < ma; run2++ ){
03607                    G2[diff_index[md+run2]] = k [newtonsteps][run1][md+run2];
03608                    G3[diff_index[md+run2]] = k2[newtonsteps][run1][md+run2];
03609                }
03610                for( run2 = 0; run2 < md; run2++ ){
03611                    G2[ddiff_index[run2]] = k [newtonsteps][run1][run2];
03612                    G3[ddiff_index[run2]] = k2[newtonsteps][run1][run2];
03613                }
03614 
03615                if( rhs[0].AD_forward2( 3*run1+newtonsteps, G2, G3, F, F2 )
03616                                       != SUCCESSFUL_RETURN ){
03617                    ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03618                    return;
03619                }
03620 
03621 
03622                applyNewtonStep( M_index[run1], k[newtonsteps+1][run1],
03623                                 k[newtonsteps][run1]  ,
03624                                *M[M_index[run1]],
03625                                 F );
03626                applyNewtonStep( M_index[run1], k2[newtonsteps+1][run1],
03627                                 k2[newtonsteps][run1]  ,
03628                                *M[M_index[run1]],
03629                                 F2 );
03630 
03631                newtonsteps++;
03632            }
03633        }
03634 
03635     // determine nablaG 2,3:
03636     // ----------------------------------------------
03637 
03638        for( run1 = 0; run1 < nstep-1; run1++ ){
03639            for( run2 = 0; run2 < md; run2++ ){
03640                nablaG2(nstep-2-run1,run2) = etaG2[run2];
03641                nablaG3(nstep-2-run1,run2) = etaG3[run2];
03642                for( run3 = 0; run3 <= run1+3; run3++ ){
03643                    nablaG2(nstep-2-run1,run2) = nablaG2(nstep-2-run1,run2) +
03644                      h[0]*A[run1+3][run3]*k [nOfNewtonSteps[run3]][run3][run2];
03645                    nablaG3(nstep-2-run1,run2) = nablaG3(nstep-2-run1,run2) +
03646                      h[0]*A[run1+3][run3]*k2[nOfNewtonSteps[run3]][run3][run2];
03647                }
03648            }
03649            for( run2 = 0; run2 < ma; run2++ ){
03650                nablaG2(nstep-2-run1,md+run2) = k [nOfNewtonSteps[run1+3]][run1+3][md+run2];
03651                nablaG3(nstep-2-run1,md+run2) = k2[nOfNewtonSteps[run1+3]][run1+3][md+run2];
03652            }
03653         }
03654 }
03655 
03656 
03657 void IntegratorBDF::determineRKEtaHBackward2(){
03658 
03659     int run1, run2, run3, run4, newtonsteps;
03660     const double hhh = c[3]*h[0];
03661     int number_;
03662 
03663     double *H1store = new double[ndir];
03664     double *H2store = new double[ndir];
03665 
03666     const double scalT = rel_time_scale/hhh;
03667 
03668     for( run4 = 0; run4 < nBDirs2; run4++ ){
03669 
03670         for( run3 = 0; run3 < 4; run3++ ){
03671             for( run1 = 0; run1 < dim; run1++ ){
03672                 for( run2 = 0; run2 < ndir; run2++ ){
03673                     l [run3][run1][run2] = 0.0;
03674                     l2[run3][run1][run2] = 0.0;
03675                 }
03676             }
03677         }
03678 
03679         for( run1 = 0; run1 < ndir; run1++ ){
03680             H1store[run1] = nablaH2(0,run1);
03681             H2store[run1] = nablaH3(0,run1);
03682         }
03683 
03684         for( run1 = 0; run1 < ndir; run1++ ){
03685 
03686             zH2[5][run1] = nablaH2(3,run1)*scalT*scalT;
03687             zH2[4][run1] = nablaH2(2,run1)*scalT + zH2[5][run1]/(3.0);
03688             zH2[3][run1] = - zH2[5][run1]/(3.0);
03689             zH2[2][run1] = nablaH2(1,run1) + zH2[4][run1]/(2.0);
03690             zH2[1][run1] = (zH2[3][run1]-zH2[4][run1])/(2.0);
03691             zH2[0][run1] = - zH2[3][run1]/(2.0);
03692 
03693             zH3[5][run1] = nablaH3(3,run1)*scalT*scalT;
03694             zH3[4][run1] = nablaH3(2,run1)*scalT + zH3[5][run1]/(3.0);
03695             zH3[3][run1] = - zH3[5][run1]/(3.0);
03696             zH3[2][run1] = nablaH3(1,run1) + zH3[4][run1]/(2.0);
03697             zH3[1][run1] = (zH3[3][run1]-zH3[4][run1])/(2.0);
03698             zH3[0][run1] = - zH3[3][run1]/(2.0);
03699 
03700             nablaH2(3,run1) = nablaH2(0,run1)*h[0] + zH2[2][run1]*(h[0]/hhh);
03701             nablaH2(2,run1) = (zH2[1][run1]-zH2[2][run1])*(h[0]/hhh);
03702             nablaH2(1,run1) = (zH2[0][run1]-zH2[1][run1])*(h[0]/hhh);
03703             nablaH2(0,run1) = -zH2[0][run1]*(h[0]/hhh);
03704 
03705             nablaH3(3,run1) = nablaH3(0,run1)*h[0] + zH3[2][run1]*(h[0]/hhh);
03706             nablaH3(2,run1) = (zH3[1][run1]-zH3[2][run1])*(h[0]/hhh);
03707             nablaH3(1,run1) = (zH3[0][run1]-zH3[1][run1])*(h[0]/hhh);
03708             nablaH3(0,run1) = -zH3[0][run1]*(h[0]/hhh);
03709         }
03710 
03711 
03712         for( run1 = 0; run1 < ndir; run1++ ){
03713             kH2[6][nOfNewtonSteps[6]][run1] = nablaH2(3,run1)/h[0];
03714             kH3[6][nOfNewtonSteps[6]][run1] = nablaH3(3,run1)/h[0];
03715         }
03716         for( run1 = 0; run1 < md; run1++ ){
03717             kH2[6][nOfNewtonSteps[6]][diff_index[run1]] =
03718             kH2[6][nOfNewtonSteps[6]][diff_index[run1]]*h[0]*A[6][6];
03719             kH3[6][nOfNewtonSteps[6]][diff_index[run1]] =
03720             kH3[6][nOfNewtonSteps[6]][diff_index[run1]]*h[0]*A[6][6];
03721         }
03722 
03723         newtonsteps = nOfNewtonSteps[6];
03724         newtonsteps--;
03725 
03726         number_ = 6;
03727 
03728         while( newtonsteps >= 0 ){
03729 
03730             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
03731                             *M[M_index[number_]],
03732                              H2 );
03733 
03734             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
03735                             *M[M_index[number_]],
03736                              H3 );
03737 
03738 
03739             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
03740                                     l[newtonsteps][number_], l2[newtonsteps][number_] )
03741                 != SUCCESSFUL_RETURN ){
03742                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03743             }
03744 
03745 
03746             for( run1 = 0; run1 < ndir; run1++ ){
03747                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
03748                                                  - l[newtonsteps][number_][run1];
03749                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
03750                                                  - l2[newtonsteps][number_][run1];
03751             }
03752             for( run1 = 0; run1 < md; run1++ ){
03753                 kH2[number_][newtonsteps][diff_index[run1]] = 
03754                    kH2[number_][newtonsteps+1][diff_index[run1]]
03755                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03756                  - l[newtonsteps][number_][ddiff_index[run1]];
03757                 kH3[number_][newtonsteps][diff_index[run1]] = 
03758                    kH3[number_][newtonsteps+1][diff_index[run1]]
03759                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03760                  - l2[newtonsteps][number_][ddiff_index[run1]];
03761             }
03762 
03763             newtonsteps--;
03764         }
03765 
03766         for( run1 = 0; run1 < ndir; run1++ ){
03767             kH2[5][nOfNewtonSteps[5]][run1] = kH2[6][0][run1] + nablaH2(2,run1)/h[0];
03768             kH3[5][nOfNewtonSteps[5]][run1] = kH3[6][0][run1] + nablaH3(2,run1)/h[0];
03769         }
03770         for( run1 = 0; run1 < md; run1++ ){
03771             kH2[5][nOfNewtonSteps[5]][diff_index[run1]] = kH2[6][0][diff_index[run1]]
03772             +  nablaH2(3,diff_index[run1])*A[6][5]
03773             +  nablaH2(2,diff_index[run1])*A[5][5];
03774             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03775                kH2[5][nOfNewtonSteps[5]][diff_index[run1]] =
03776                kH2[5][nOfNewtonSteps[5]][diff_index[run1]]
03777                -l[run3][number_][diff_index[run1]]*h[0]*A[6][5];
03778             }
03779             kH3[5][nOfNewtonSteps[5]][diff_index[run1]] = kH3[6][0][diff_index[run1]]
03780             +  nablaH3(3,diff_index[run1])*A[6][5]
03781             +  nablaH3(2,diff_index[run1])*A[5][5];
03782             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03783                kH3[5][nOfNewtonSteps[5]][diff_index[run1]] =
03784                kH3[5][nOfNewtonSteps[5]][diff_index[run1]]
03785                -l2[run3][number_][diff_index[run1]]*h[0]*A[6][5];
03786             }
03787         }
03788 
03789 
03790         newtonsteps = nOfNewtonSteps[5];
03791         newtonsteps--;
03792 
03793         number_ = 5;
03794 
03795         while( newtonsteps >= 0 ){
03796 
03797             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
03798                             *M[M_index[number_]],
03799                              H2 );
03800             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
03801                             *M[M_index[number_]],
03802                              H3 );
03803 
03804             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
03805                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
03806                 != SUCCESSFUL_RETURN ){
03807                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03808             }
03809 
03810 
03811             for( run1 = 0; run1 < ndir; run1++ ){
03812                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
03813                                                  - l[newtonsteps][number_][run1];
03814                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
03815                                                  - l2[newtonsteps][number_][run1];
03816             }
03817             for( run1 = 0; run1 < md; run1++ ){
03818                 kH2[number_][newtonsteps][diff_index[run1]] = 
03819                    kH2[number_][newtonsteps+1][diff_index[run1]]
03820                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03821                  - l[newtonsteps][number_][ddiff_index[run1]];
03822                 kH3[number_][newtonsteps][diff_index[run1]] = 
03823                    kH3[number_][newtonsteps+1][diff_index[run1]]
03824                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03825                  - l2[newtonsteps][number_][ddiff_index[run1]];
03826             }
03827 
03828             newtonsteps--;
03829         }
03830 
03831 
03832         for( run1 = 0; run1 < ndir; run1++ ){
03833             kH2[4][nOfNewtonSteps[4]][run1] = kH2[5][0][run1] + nablaH2(1,run1)/h[0];
03834             kH3[4][nOfNewtonSteps[4]][run1] = kH3[5][0][run1] + nablaH3(1,run1)/h[0];
03835         }
03836         for( run1 = 0; run1 < md; run1++ ){
03837             kH2[4][nOfNewtonSteps[4]][diff_index[run1]] = kH2[5][0][diff_index[run1]]
03838             +  nablaH2(3,diff_index[run1])*A[6][4]
03839             +  nablaH2(2,diff_index[run1])*A[5][4]
03840             +  nablaH2(1,diff_index[run1])*A[4][4];
03841             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03842                kH2[4][nOfNewtonSteps[4]][diff_index[run1]] =
03843                kH2[4][nOfNewtonSteps[4]][diff_index[run1]]
03844                -l[run3][5][diff_index[run1]]*h[0]*A[5][4];
03845             }
03846             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03847                kH2[4][nOfNewtonSteps[4]][diff_index[run1]] =
03848                kH2[4][nOfNewtonSteps[4]][diff_index[run1]]
03849                -l[run3][6][diff_index[run1]]*h[0]*A[6][4];
03850             }
03851             kH3[4][nOfNewtonSteps[4]][diff_index[run1]] = kH3[5][0][diff_index[run1]]
03852             +  nablaH3(3,diff_index[run1])*A[6][4]
03853             +  nablaH3(2,diff_index[run1])*A[5][4]
03854             +  nablaH3(1,diff_index[run1])*A[4][4];
03855             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03856                kH3[4][nOfNewtonSteps[4]][diff_index[run1]] =
03857                kH3[4][nOfNewtonSteps[4]][diff_index[run1]]
03858                -l2[run3][5][diff_index[run1]]*h[0]*A[5][4];
03859             }
03860             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03861                kH3[4][nOfNewtonSteps[4]][diff_index[run1]] =
03862                kH3[4][nOfNewtonSteps[4]][diff_index[run1]]
03863                -l2[run3][6][diff_index[run1]]*h[0]*A[6][4];
03864             }
03865         }
03866 
03867 
03868         newtonsteps = nOfNewtonSteps[4];
03869         newtonsteps--;
03870 
03871         number_ = 4;
03872 
03873         while( newtonsteps >= 0 ){
03874 
03875             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
03876                             *M[M_index[number_]],
03877                              H2 );
03878             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
03879                             *M[M_index[number_]],
03880                              H3 );
03881 
03882             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
03883                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
03884                 != SUCCESSFUL_RETURN ){
03885                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03886             }
03887 
03888 
03889             for( run1 = 0; run1 < ndir; run1++ ){
03890                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
03891                                                  - l[newtonsteps][number_][run1];
03892                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
03893                                                  - l2[newtonsteps][number_][run1];
03894             }
03895             for( run1 = 0; run1 < md; run1++ ){
03896                 kH2[number_][newtonsteps][diff_index[run1]] = 
03897                    kH2[number_][newtonsteps+1][diff_index[run1]]
03898                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03899                  - l[newtonsteps][number_][ddiff_index[run1]];
03900                 kH3[number_][newtonsteps][diff_index[run1]] = 
03901                    kH3[number_][newtonsteps+1][diff_index[run1]]
03902                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03903                  - l2[newtonsteps][number_][ddiff_index[run1]];
03904             }
03905 
03906             newtonsteps--;
03907         }
03908 
03909 
03910         for( run1 = 0; run1 < ndir; run1++ ){
03911             kH2[3][nOfNewtonSteps[3]][run1] = kH2[4][0][run1] + nablaH2(0,run1)/h[0];
03912             kH3[3][nOfNewtonSteps[3]][run1] = kH3[4][0][run1] + nablaH3(0,run1)/h[0];
03913         }
03914         for( run1 = 0; run1 < md; run1++ ){
03915             kH2[3][nOfNewtonSteps[3]][diff_index[run1]] = kH2[4][0][diff_index[run1]]
03916             +  nablaH2(3,diff_index[run1])*A[6][3]
03917             +  nablaH2(2,diff_index[run1])*A[5][3]
03918             +  nablaH2(1,diff_index[run1])*A[4][3]
03919             +  nablaH2(0,diff_index[run1])*A[3][3];
03920             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03921                kH2[3][nOfNewtonSteps[3]][diff_index[run1]] =
03922                kH2[3][nOfNewtonSteps[3]][diff_index[run1]]
03923                -l[run3][4][diff_index[run1]]*h[0]*A[4][3];
03924             }
03925             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03926                kH2[3][nOfNewtonSteps[3]][diff_index[run1]] =
03927                kH2[3][nOfNewtonSteps[3]][diff_index[run1]]
03928                -l[run3][5][diff_index[run1]]*h[0]*A[5][3];
03929             }
03930             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03931                kH2[3][nOfNewtonSteps[3]][diff_index[run1]] =
03932                kH2[3][nOfNewtonSteps[3]][diff_index[run1]]
03933                -l[run3][6][diff_index[run1]]*h[0]*A[6][3];
03934             }
03935             kH3[3][nOfNewtonSteps[3]][diff_index[run1]] = kH3[4][0][diff_index[run1]]
03936             +  nablaH3(3,diff_index[run1])*A[6][3]
03937             +  nablaH3(2,diff_index[run1])*A[5][3]
03938             +  nablaH3(1,diff_index[run1])*A[4][3]
03939             +  nablaH3(0,diff_index[run1])*A[3][3];
03940             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
03941                kH3[3][nOfNewtonSteps[3]][diff_index[run1]] =
03942                kH3[3][nOfNewtonSteps[3]][diff_index[run1]]
03943                -l2[run3][4][diff_index[run1]]*h[0]*A[4][3];
03944             }
03945             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
03946                kH3[3][nOfNewtonSteps[3]][diff_index[run1]] =
03947                kH3[3][nOfNewtonSteps[3]][diff_index[run1]]
03948                -l2[run3][5][diff_index[run1]]*h[0]*A[5][3];
03949             }
03950             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
03951                kH3[3][nOfNewtonSteps[3]][diff_index[run1]] =
03952                kH3[3][nOfNewtonSteps[3]][diff_index[run1]]
03953                -l2[run3][6][diff_index[run1]]*h[0]*A[6][3];
03954             }
03955         }
03956 
03957 
03958         newtonsteps = nOfNewtonSteps[3];
03959         newtonsteps--;
03960 
03961         number_ = 3;
03962 
03963         while( newtonsteps >= 0 ){
03964 
03965             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
03966                             *M[M_index[number_]],
03967                              H2 );
03968             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
03969                             *M[M_index[number_]],
03970                              H3 );
03971 
03972             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
03973                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
03974                 != SUCCESSFUL_RETURN ){
03975                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
03976             }
03977 
03978 
03979             for( run1 = 0; run1 < ndir; run1++ ){
03980                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
03981                                                  - l[newtonsteps][number_][run1];
03982                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
03983                                                  - l2[newtonsteps][number_][run1];
03984             }
03985             for( run1 = 0; run1 < md; run1++ ){
03986                 kH2[number_][newtonsteps][diff_index[run1]] = 
03987                    kH2[number_][newtonsteps+1][diff_index[run1]]
03988                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03989                  - l[newtonsteps][number_][ddiff_index[run1]];
03990                 kH3[number_][newtonsteps][diff_index[run1]] = 
03991                    kH3[number_][newtonsteps+1][diff_index[run1]]
03992                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
03993                  - l2[newtonsteps][number_][ddiff_index[run1]];
03994             }
03995 
03996             newtonsteps--;
03997         }
03998 
03999 
04000         for( run1 = 0; run1 < ndir; run1++ ){
04001             kH2[2][nOfNewtonSteps[2]][run1] = kH2[3][0][run1];
04002             kH3[2][nOfNewtonSteps[2]][run1] = kH3[3][0][run1];
04003         }
04004         for( run1 = 0; run1 < md; run1++ ){
04005             kH2[2][nOfNewtonSteps[2]][diff_index[run1]] = kH2[3][0][diff_index[run1]]
04006             +  nablaH2(3,diff_index[run1])*A[6][2]
04007             +  nablaH2(2,diff_index[run1])*A[5][2]
04008             +  nablaH2(1,diff_index[run1])*A[4][2]
04009             +  nablaH2(0,diff_index[run1])*A[3][2];
04010             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04011                kH2[2][nOfNewtonSteps[2]][diff_index[run1]] =
04012                kH2[2][nOfNewtonSteps[2]][diff_index[run1]]
04013                -l[run3][3][diff_index[run1]]*h[0]*A[3][2];
04014             }
04015             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04016                kH2[2][nOfNewtonSteps[2]][diff_index[run1]] =
04017                kH2[2][nOfNewtonSteps[2]][diff_index[run1]]
04018                -l[run3][4][diff_index[run1]]*h[0]*A[4][2];
04019             }
04020             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04021                kH2[2][nOfNewtonSteps[2]][diff_index[run1]] =
04022                kH2[2][nOfNewtonSteps[2]][diff_index[run1]]
04023                -l[run3][5][diff_index[run1]]*h[0]*A[5][2];
04024             }
04025             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04026                kH2[2][nOfNewtonSteps[2]][diff_index[run1]] =
04027                kH2[2][nOfNewtonSteps[2]][diff_index[run1]]
04028                -l[run3][6][diff_index[run1]]*h[0]*A[6][2];
04029             }
04030             kH3[2][nOfNewtonSteps[2]][diff_index[run1]] = kH3[3][0][diff_index[run1]]
04031             +  nablaH3(3,diff_index[run1])*A[6][2]
04032             +  nablaH3(2,diff_index[run1])*A[5][2]
04033             +  nablaH3(1,diff_index[run1])*A[4][2]
04034             +  nablaH3(0,diff_index[run1])*A[3][2];
04035             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04036                kH3[2][nOfNewtonSteps[2]][diff_index[run1]] =
04037                kH3[2][nOfNewtonSteps[2]][diff_index[run1]]
04038                -l2[run3][3][diff_index[run1]]*h[0]*A[3][2];
04039             }
04040             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04041                kH3[2][nOfNewtonSteps[2]][diff_index[run1]] =
04042                kH3[2][nOfNewtonSteps[2]][diff_index[run1]]
04043                -l2[run3][4][diff_index[run1]]*h[0]*A[4][2];
04044             }
04045             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04046                kH3[2][nOfNewtonSteps[2]][diff_index[run1]] =
04047                kH3[2][nOfNewtonSteps[2]][diff_index[run1]]
04048                -l2[run3][5][diff_index[run1]]*h[0]*A[5][2];
04049             }
04050             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04051                kH3[2][nOfNewtonSteps[2]][diff_index[run1]] =
04052                kH3[2][nOfNewtonSteps[2]][diff_index[run1]]
04053                -l2[run3][6][diff_index[run1]]*h[0]*A[6][2];
04054             }
04055         }
04056 
04057 
04058         newtonsteps = nOfNewtonSteps[2];
04059         newtonsteps--;
04060 
04061         number_ = 2;
04062 
04063         while( newtonsteps >= 0 ){
04064 
04065             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
04066                             *M[M_index[number_]],
04067                              H2 );
04068             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
04069                             *M[M_index[number_]],
04070                              H3 );
04071 
04072             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
04073                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
04074                 != SUCCESSFUL_RETURN ){
04075                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04076             }
04077 
04078 
04079             for( run1 = 0; run1 < ndir; run1++ ){
04080                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
04081                                                  - l[newtonsteps][number_][run1];
04082                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
04083                                                  - l2[newtonsteps][number_][run1];
04084             }
04085             for( run1 = 0; run1 < md; run1++ ){
04086                 kH2[number_][newtonsteps][diff_index[run1]] = 
04087                    kH2[number_][newtonsteps+1][diff_index[run1]]
04088                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04089                  - l[newtonsteps][number_][ddiff_index[run1]];
04090                 kH3[number_][newtonsteps][diff_index[run1]] = 
04091                    kH3[number_][newtonsteps+1][diff_index[run1]]
04092                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04093                  - l2[newtonsteps][number_][ddiff_index[run1]];
04094             }
04095 
04096             newtonsteps--;
04097         }
04098 
04099 
04100         for( run1 = 0; run1 < ndir; run1++ ){
04101             kH2[1][nOfNewtonSteps[1]][run1] = kH2[2][0][run1];
04102             kH3[1][nOfNewtonSteps[1]][run1] = kH3[2][0][run1];
04103         }
04104         for( run1 = 0; run1 < md; run1++ ){
04105             kH2[1][nOfNewtonSteps[1]][diff_index[run1]] = kH2[2][0][diff_index[run1]]
04106             +  nablaH2(3,diff_index[run1])*A[6][1]
04107             +  nablaH2(2,diff_index[run1])*A[5][1]
04108             +  nablaH2(1,diff_index[run1])*A[4][1]
04109             +  nablaH2(0,diff_index[run1])*A[3][1];
04110             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
04111                kH2[1][nOfNewtonSteps[1]][diff_index[run1]] =
04112                kH2[1][nOfNewtonSteps[1]][diff_index[run1]]
04113                -l[run3][2][diff_index[run1]]*h[0]*A[2][1];
04114             }
04115             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04116                kH2[1][nOfNewtonSteps[1]][diff_index[run1]] =
04117                kH2[1][nOfNewtonSteps[1]][diff_index[run1]]
04118                -l[run3][3][diff_index[run1]]*h[0]*A[3][1];
04119             }
04120             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04121                kH2[1][nOfNewtonSteps[1]][diff_index[run1]] =
04122                kH2[1][nOfNewtonSteps[1]][diff_index[run1]]
04123                -l[run3][4][diff_index[run1]]*h[0]*A[4][1];
04124             }
04125             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04126                kH2[1][nOfNewtonSteps[1]][diff_index[run1]] =
04127                kH2[1][nOfNewtonSteps[1]][diff_index[run1]]
04128                -l[run3][5][diff_index[run1]]*h[0]*A[5][1];
04129             }
04130             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04131                kH2[1][nOfNewtonSteps[1]][diff_index[run1]] =
04132                kH2[1][nOfNewtonSteps[1]][diff_index[run1]]
04133                -l[run3][6][diff_index[run1]]*h[0]*A[6][1];
04134             }
04135             kH3[1][nOfNewtonSteps[1]][diff_index[run1]] = kH3[2][0][diff_index[run1]]
04136             +  nablaH3(3,diff_index[run1])*A[6][1]
04137             +  nablaH3(2,diff_index[run1])*A[5][1]
04138             +  nablaH3(1,diff_index[run1])*A[4][1]
04139             +  nablaH3(0,diff_index[run1])*A[3][1];
04140             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
04141                kH3[1][nOfNewtonSteps[1]][diff_index[run1]] =
04142                kH3[1][nOfNewtonSteps[1]][diff_index[run1]]
04143                -l2[run3][2][diff_index[run1]]*h[0]*A[2][1];
04144             }
04145             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04146                kH3[1][nOfNewtonSteps[1]][diff_index[run1]] =
04147                kH3[1][nOfNewtonSteps[1]][diff_index[run1]]
04148                -l2[run3][3][diff_index[run1]]*h[0]*A[3][1];
04149             }
04150             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04151                kH3[1][nOfNewtonSteps[1]][diff_index[run1]] =
04152                kH3[1][nOfNewtonSteps[1]][diff_index[run1]]
04153                -l2[run3][4][diff_index[run1]]*h[0]*A[4][1];
04154             }
04155             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04156                kH3[1][nOfNewtonSteps[1]][diff_index[run1]] =
04157                kH3[1][nOfNewtonSteps[1]][diff_index[run1]]
04158                -l2[run3][5][diff_index[run1]]*h[0]*A[5][1];
04159             }
04160             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04161                kH3[1][nOfNewtonSteps[1]][diff_index[run1]] =
04162                kH3[1][nOfNewtonSteps[1]][diff_index[run1]]
04163                -l2[run3][6][diff_index[run1]]*h[0]*A[6][1];
04164             }
04165         }
04166 
04167 
04168         newtonsteps = nOfNewtonSteps[1];
04169         newtonsteps--;
04170 
04171         number_ = 1;
04172 
04173         while( newtonsteps >= 0 ){
04174 
04175             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
04176                             *M[M_index[number_]],
04177                              H2 );
04178             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
04179                             *M[M_index[number_]],
04180                              H3 );
04181 
04182             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
04183                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
04184                 != SUCCESSFUL_RETURN ){
04185                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04186             }
04187 
04188 
04189             for( run1 = 0; run1 < ndir; run1++ ){
04190                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
04191                                                  - l[newtonsteps][number_][run1];
04192                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
04193                                                  - l2[newtonsteps][number_][run1];
04194             }
04195             for( run1 = 0; run1 < md; run1++ ){
04196                 kH2[number_][newtonsteps][diff_index[run1]] = 
04197                    kH2[number_][newtonsteps+1][diff_index[run1]]
04198                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04199                  - l[newtonsteps][number_][ddiff_index[run1]];
04200                 kH3[number_][newtonsteps][diff_index[run1]] = 
04201                    kH3[number_][newtonsteps+1][diff_index[run1]]
04202                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04203                  - l2[newtonsteps][number_][ddiff_index[run1]];
04204             }
04205 
04206             newtonsteps--;
04207         }
04208 
04209 
04210 
04211         for( run1 = 0; run1 < ndir; run1++ ){
04212             kH2[0][nOfNewtonSteps[0]][run1] = kH2[1][0][run1];
04213             kH3[0][nOfNewtonSteps[0]][run1] = kH3[1][0][run1];
04214         }
04215         for( run1 = 0; run1 < md; run1++ ){
04216             kH2[0][nOfNewtonSteps[0]][diff_index[run1]] = kH2[1][0][diff_index[run1]]
04217             +  nablaH2(3,diff_index[run1])*A[6][0]
04218             +  nablaH2(2,diff_index[run1])*A[5][0]
04219             +  nablaH2(1,diff_index[run1])*A[4][0]
04220             +  nablaH2(0,diff_index[run1])*A[3][0];
04221             for( run3 = 0; run3 < nOfNewtonSteps[1]; run3++){
04222                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04223                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04224                -l[run3][1][diff_index[run1]]*h[0]*A[1][0];
04225             }
04226             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
04227                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04228                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04229                -l[run3][2][diff_index[run1]]*h[0]*A[2][0];
04230             }
04231             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04232                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04233                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04234                -l[run3][3][diff_index[run1]]*h[0]*A[3][0];
04235             }
04236             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04237                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04238                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04239                -l[run3][4][diff_index[run1]]*h[0]*A[4][0];
04240             }
04241             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04242                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04243                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04244                -l[run3][5][diff_index[run1]]*h[0]*A[5][0];
04245             }
04246             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04247                kH2[0][nOfNewtonSteps[0]][diff_index[run1]] =
04248                kH2[0][nOfNewtonSteps[0]][diff_index[run1]]
04249                -l[run3][6][diff_index[run1]]*h[0]*A[6][0];
04250             }
04251             kH3[0][nOfNewtonSteps[0]][diff_index[run1]] = kH3[1][0][diff_index[run1]]
04252             +  nablaH3(3,diff_index[run1])*A[6][0]
04253             +  nablaH3(2,diff_index[run1])*A[5][0]
04254             +  nablaH3(1,diff_index[run1])*A[4][0]
04255             +  nablaH3(0,diff_index[run1])*A[3][0];
04256             for( run3 = 0; run3 < nOfNewtonSteps[1]; run3++){
04257                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04258                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04259                -l2[run3][1][diff_index[run1]]*h[0]*A[1][0];
04260             }
04261             for( run3 = 0; run3 < nOfNewtonSteps[2]; run3++){
04262                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04263                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04264                -l2[run3][2][diff_index[run1]]*h[0]*A[2][0];
04265             }
04266             for( run3 = 0; run3 < nOfNewtonSteps[3]; run3++){
04267                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04268                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04269                -l2[run3][3][diff_index[run1]]*h[0]*A[3][0];
04270             }
04271             for( run3 = 0; run3 < nOfNewtonSteps[4]; run3++){
04272                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04273                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04274                -l2[run3][4][diff_index[run1]]*h[0]*A[4][0];
04275             }
04276             for( run3 = 0; run3 < nOfNewtonSteps[5]; run3++){
04277                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04278                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04279                -l2[run3][5][diff_index[run1]]*h[0]*A[5][0];
04280             }
04281             for( run3 = 0; run3 < nOfNewtonSteps[6]; run3++){
04282                kH3[0][nOfNewtonSteps[0]][diff_index[run1]] =
04283                kH3[0][nOfNewtonSteps[0]][diff_index[run1]]
04284                -l2[run3][6][diff_index[run1]]*h[0]*A[6][0];
04285             }
04286         }
04287 
04288         newtonsteps = nOfNewtonSteps[0];
04289         newtonsteps--;
04290 
04291         number_ = 0;
04292 
04293         while( newtonsteps >= 0 ){
04294 
04295             applyMTranspose( M_index[number_], kH2[number_][newtonsteps+1],
04296                             *M[M_index[number_]],
04297                              H2 );
04298             applyMTranspose( M_index[number_], kH3[number_][newtonsteps+1],
04299                             *M[M_index[number_]],
04300                              H3 );
04301 
04302             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
04303                                      l[newtonsteps][number_], l2[newtonsteps][number_] )
04304                 != SUCCESSFUL_RETURN ){
04305                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04306             }
04307 
04308 
04309             for( run1 = 0; run1 < ndir; run1++ ){
04310                 kH2[number_][newtonsteps][run1] =   kH2[number_][newtonsteps+1][run1]
04311                                                  - l[newtonsteps][number_][run1];
04312                 kH3[number_][newtonsteps][run1] =   kH3[number_][newtonsteps+1][run1]
04313                                                  - l2[newtonsteps][number_][run1];
04314             }
04315             for( run1 = 0; run1 < md; run1++ ){
04316                 kH2[number_][newtonsteps][diff_index[run1]] = 
04317                    kH2[number_][newtonsteps+1][diff_index[run1]]
04318                  - l[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04319                  - l[newtonsteps][number_][ddiff_index[run1]];
04320                 kH3[number_][newtonsteps][diff_index[run1]] = 
04321                    kH3[number_][newtonsteps+1][diff_index[run1]]
04322                  - l2[newtonsteps][number_][diff_index[run1]]*h[0]*A[number_][number_]
04323                  - l2[newtonsteps][number_][ddiff_index[run1]];
04324             }
04325 
04326             newtonsteps--;
04327         }
04328 
04329 
04330         for( run1 = 0; run1 < ndir; run1++ ){
04331             nablaH2(0,run1) = H1store[run1];
04332             nablaH3(0,run1) = H2store[run1];
04333         }
04334 
04335         for( run1 = 0; run1 < ndir; run1++ ){
04336             for( run3 = 0; run3 < nOfNewtonSteps[0]; run3++){
04337                for( run2 = 0; run2 < dim; run2++ ){
04338                    nablaH2(0,run1) = nablaH2(0,run1)
04339                    -l[run3][run2][run1];
04340                    nablaH3(0,run1) = nablaH3(0,run1)
04341                   -l2[run3][run2][run1];
04342                }
04343             }
04344         }
04345 
04346     }
04347     delete[] H1store;
04348     delete[] H2store;
04349 }
04350 
04351 
04352 
04353 void IntegratorBDF::determineBDFEtaGForward( int number_ ){
04354 
04355     int run1, run2, run4;
04356     int newtonsteps;
04357     double pp;
04358 
04359     const double scalT = psi[number_][0]/rel_time_scale;
04360 
04361     for( run4 = 0; run4 < nFDirs; run4++ ){
04362 
04363         for( run1 = 0; run1 < m; run1++ ){
04364             phiG(1,run1) = nablaG(1,run1)*scalT;
04365         }
04366 
04367         pp = (psi[number_][1]/psi[number_][0]);
04368         for( run1 = 0; run1 < m; run1++ ){
04369             phiG(2,run1) = pp*nablaG(2,run1)*scalT*scalT;
04370         }
04371 
04372         pp = pp*(psi[number_][2]/psi[number_][0]);
04373         for( run1 = 0; run1 < m; run1++ ){
04374             phiG(3,run1) = pp*nablaG(3,run1)*scalT*scalT*scalT;
04375         }
04376 
04377         pp = pp*(psi[number_][3]/psi[number_][0]);
04378         for( run1 = 0; run1 < m; run1++ ){
04379             phiG(4,run1) = pp*nablaG(4,run1)*scalT*scalT*scalT*scalT;
04380         }
04381 
04382         for( run1 = 0; run1 < m; run1++ ){
04383             deltaG(1,run1) = nablaG(0,run1) + phiG(1,run1);
04384         }
04385 
04386         for( run1 = 0; run1 < m; run1++ ){
04387             deltaG(2,run1) = deltaG(1,run1) + phiG(2,run1);
04388         }
04389 
04390         for( run1 = 0; run1 < m; run1++ ){
04391             deltaG(3,run1) = deltaG(2,run1) + phiG(3,run1);
04392         }
04393 
04394         for( run1 = 0; run1 < m; run1++ ){
04395             eta[0][run1] = deltaG(3,run1) + phiG(4,run1);
04396         }
04397 
04398 
04399         for( run1 = 0; run1 < md; run1++ ){
04400 
04401             c2G(run1) = - nablaG (0,run1)/psi[number_][0]
04402                         - deltaG (1,run1)/psi[number_][1]
04403                         - deltaG (2,run1)/psi[number_][2]
04404                         - deltaG (3,run1)/psi[number_][3];
04405         }
04406 
04407 
04408         newtonsteps = 0;
04409         while( newtonsteps < nOfNewtonSteps[number_] ){
04410             for( run2 = 0; run2 < m; run2++ ){
04411                  G[diff_index[run2]] = eta[newtonsteps][run2];
04412             }
04413             for( run2 = 0; run2 < md; run2++ ){
04414                  G[ddiff_index[run2]] =
04415                  gamma[number_][4]*eta[newtonsteps][run2]+c2G(run2);
04416             }
04417 
04418             if( rhs[0].AD_forward( 3*number_+newtonsteps, G, F )
04419                 != SUCCESSFUL_RETURN ){
04420                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04421             }
04422 
04423             applyNewtonStep( M_index[number_], eta[newtonsteps+1],
04424                              eta[newtonsteps],
04425                             *M[M_index[number_]],
04426                              F );
04427 
04428             newtonsteps++;
04429         }
04430 
04431         for( run1 = 0; run1 < m; run1++ ){
04432             nablaY_(0,run1) = eta[nOfNewtonSteps[number_]][run1];
04433             nablaY_(1,run1) = nablaY_(0,run1) - nablaG(0,run1);
04434             nablaY_(2,run1) = (nablaY_(1,run1) - nablaG(1,run1)*scalT)*(psi[number_][0]/psi[number_][1]);
04435             nablaY_(3,run1) = (nablaY_(2,run1) - nablaG(2,run1)*scalT*scalT)*(psi[number_][0]/psi[number_][2]);
04436             nablaY_(4,run1) = (nablaY_(3,run1) - nablaG(3,run1)*scalT*scalT*scalT)*(psi[number_][0]/psi[number_][3]);
04437         }
04438 
04439         nablaG = nablaY_;
04440     }
04441 }
04442 
04443 
04444 void IntegratorBDF::determineBDFEtaHBackward( int number_ ){
04445 
04446     int run1, run2;
04447     int newtonsteps;
04448     double pp;
04449 
04450     const double scalT = rel_time_scale/psi[number_][0];
04451 
04452     newtonsteps = nOfNewtonSteps[number_];
04453 
04454     for( run1 = 0; run1 < newtonsteps; run1++ ){
04455         for( run2 = 0; run2 < ndir; run2++ ){
04456             l[run1][0][run2] = 0.0;
04457         }
04458     }
04459 
04460     for( run1 = 0; run1 < ndir; run1++ ){
04461 
04462         nablaH_(4,run1) = nablaH(4,run1)*scalT*scalT*scalT;
04463         nablaH_(3,run1) = nablaH(3,run1)*scalT*scalT + nablaH_(4,run1)*(psi[number_][0]/psi[number_][3]);
04464         nablaH_(2,run1) = nablaH(2,run1)*scalT + nablaH_(3,run1)*(psi[number_][0]/psi[number_][2]);
04465         nablaH_(1,run1) = nablaH(1,run1) + nablaH_(2,run1)*(psi[number_][0]/psi[number_][1]);
04466         nablaH_(0,run1) = nablaH(0,run1) + nablaH_(1,run1)/psi[number_][0];
04467 
04468         etaH[newtonsteps][run1] = nablaH_(0,run1);
04469         c2H(run1) = 0.0;
04470     }
04471 
04472     newtonsteps--;
04473     while( newtonsteps >= 0 ){
04474 
04475         applyMTranspose( M_index[number_], etaH[newtonsteps+1], M[M_index[number_]][0], H );
04476 
04477         if( rhs[0].AD_backward( 3*number_+newtonsteps, H, l[newtonsteps][0] ) != SUCCESSFUL_RETURN )
04478             ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04479 
04480         for( run1 = 0; run1 < ndir; run1++ ){
04481             etaH[newtonsteps][run1] = etaH[newtonsteps+1][run1] - l[newtonsteps][0][run1];
04482         }
04483 
04484         for( run1 = 0; run1 < md; run1++ )
04485             etaH[newtonsteps][diff_index[run1]] -= l[newtonsteps][0][ddiff_index[run1]]*gamma[number_][4];
04486 
04487         for( run1 = 0; run1 < md; run1++ )
04488             c2H(diff_index[run1]) -= l[newtonsteps][0][ddiff_index[run1]];
04489 
04490         newtonsteps--;
04491     }
04492 
04493 
04494     for( run1 = 0; run1 < ndir; run1++ ){
04495 
04496         deltaH(3,run1) = etaH[0][run1]   - c2H(run1)/psi[number_][3];
04497         deltaH(2,run1) = deltaH(3,run1) - c2H(run1)/psi[number_][2];
04498         deltaH(1,run1) = deltaH(2,run1) - c2H(run1)/psi[number_][1];
04499     }
04500 
04501     for( run1 = 0; run1 < ndir; run1++ ){
04502 
04503         nablaH(0,run1) =  - nablaH_(1,run1) /psi[number_][0]
04504                                  - c2H(run1) /psi[number_][0]
04505                                  + deltaH(1,run1);
04506     }
04507     pp =  psi[number_][0];
04508     for( run1 = 0; run1 < ndir; run1++ ){
04509         nablaH(1,run1) =  - nablaH_(2,run1)*(psi[number_][0]/psi[number_][1])
04510                                      + pp*deltaH(1,run1);
04511     }
04512     pp = pp*(psi[number_][1]/psi[number_][0]);
04513     for( run1 = 0; run1 < ndir; run1++ ){
04514         nablaH(2,run1) =  - nablaH_(3,run1)*(psi[number_][0]/psi[number_][2])
04515                                  + pp*deltaH(2,run1);
04516     }
04517     pp = pp*(psi[number_][2]/psi[number_][0]);
04518 
04519     for( run1 = 0; run1 < ndir; run1++ ){
04520         nablaH(3,run1) =  - nablaH_(4,run1)*(psi[number_][0]/psi[number_][3])
04521                                  + pp*deltaH(3,run1);
04522     }
04523     pp = pp*(psi[number_][3]/psi[number_][0]);
04524 
04525     for( run1 = 0; run1 < ndir; run1++ ){
04526         nablaH(4,run1) =  pp*etaH[0][run1];
04527     }
04528 }
04529 
04530 
04531 
04532 void IntegratorBDF::determineBDFEtaGForward2( int number_ ){
04533 
04534     int run1, run2, run4;
04535     int newtonsteps;
04536     double pp;
04537 
04538     const double scalT = psi[number_][0]/rel_time_scale;
04539 
04540     for( run4 = 0; run4 < nFDirs2; run4++ ){
04541 
04542         for( run1 = 0; run1 < m; run1++ ){
04543             phiG2(1,run1) = nablaG2(1,run1)*scalT;
04544         }
04545         for( run1 = 0; run1 < m; run1++ ){
04546             phiG3(1,run1) = nablaG3(1,run1)*scalT;
04547         }
04548 
04549         pp = (psi[number_][1]/psi[number_][0]);
04550         for( run1 = 0; run1 < m; run1++ ){
04551             phiG2(2,run1) = pp*nablaG2(2,run1)*scalT*scalT;
04552         }
04553         for( run1 = 0; run1 < m; run1++ ){
04554             phiG3(2,run1) = pp*nablaG3(2,run1)*scalT*scalT;
04555         }
04556 
04557         pp = pp*(psi[number_][2]/psi[number_][0]);
04558         for( run1 = 0; run1 < m; run1++ ){
04559             phiG2(3,run1) = pp*nablaG2(3,run1)*scalT*scalT*scalT;
04560         }
04561         for( run1 = 0; run1 < m; run1++ ){
04562             phiG3(3,run1) = pp*nablaG3(3,run1)*scalT*scalT*scalT;
04563         }
04564 
04565         pp = pp*(psi[number_][3]/psi[number_][0]);
04566         for( run1 = 0; run1 < m; run1++ ){
04567             phiG2(4,run1) = pp*nablaG2(4,run1)*scalT*scalT*scalT*scalT;
04568         }
04569         for( run1 = 0; run1 < m; run1++ ){
04570             phiG3(4,run1) = pp*nablaG3(4,run1)*scalT*scalT*scalT*scalT;
04571         }
04572 
04573         for( run1 = 0; run1 < m; run1++ ){
04574             deltaG2(1,run1) = nablaG2(0,run1) + phiG2(1,run1);
04575         }
04576         for( run1 = 0; run1 < m; run1++ ){
04577             deltaG3(1,run1) = nablaG3(0,run1) + phiG3(1,run1);
04578         }
04579 
04580         for( run1 = 0; run1 < m; run1++ ){
04581             deltaG2(2,run1) = deltaG2(1,run1) + phiG2(2,run1);
04582         }
04583         for( run1 = 0; run1 < m; run1++ ){
04584             deltaG3(2,run1) = deltaG3(1,run1) + phiG3(2,run1);
04585         }
04586 
04587         for( run1 = 0; run1 < m; run1++ ){
04588             deltaG2(3,run1) = deltaG2(2,run1) + phiG2(3,run1);
04589         }
04590         for( run1 = 0; run1 < m; run1++ ){
04591             deltaG3(3,run1) = deltaG3(2,run1) + phiG3(3,run1);
04592         }
04593 
04594         for( run1 = 0; run1 < m; run1++ ){
04595             eta [0][run1] = deltaG2(3,run1) + phiG2(4,run1);
04596         }
04597         for( run1 = 0; run1 < m; run1++ ){
04598             eta2[0][run1] = deltaG3(3,run1) + phiG3(4,run1);
04599         }
04600 
04601         for( run1 = 0; run1 < md; run1++ ){
04602 
04603             c2G2(run1) = - nablaG2 (0,run1)/psi[number_][0]
04604                          - deltaG2 (1,run1)/psi[number_][1]
04605                          - deltaG2 (2,run1)/psi[number_][2]
04606                          - deltaG2 (3,run1)/psi[number_][3];
04607         }
04608         for( run1 = 0; run1 < md; run1++ ){
04609 
04610             c2G3(run1) = - nablaG3 (0,run1)/psi[number_][0]
04611                          - deltaG3 (1,run1)/psi[number_][1]
04612                          - deltaG3 (2,run1)/psi[number_][2]
04613                          - deltaG3 (3,run1)/psi[number_][3];
04614         }
04615 
04616 
04617         newtonsteps = 0;
04618         while( newtonsteps < nOfNewtonSteps[number_] ){
04619             for( run2 = 0; run2 < m; run2++ ){
04620                  G2[diff_index[run2]] = eta [newtonsteps][run2];
04621                  G3[diff_index[run2]] = eta2[newtonsteps][run2];
04622             }
04623             for( run2 = 0; run2 < md; run2++ ){
04624                  G2[ddiff_index[run2]] =
04625                  gamma[number_][4]*eta [newtonsteps][run2]+c2G2(run2);
04626                  G3[ddiff_index[run2]] =
04627                  gamma[number_][4]*eta2[newtonsteps][run2]+c2G3(run2);
04628             }
04629 
04630             if( rhs[0].AD_forward2( 3*number_+newtonsteps, G2, G3, F, F2 )
04631                 != SUCCESSFUL_RETURN ){
04632                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04633             }
04634 
04635             applyNewtonStep( M_index[number_], eta[newtonsteps+1],
04636                              eta[newtonsteps],
04637                             *M[M_index[number_]],
04638                              F );
04639 
04640             applyNewtonStep( M_index[number_], eta2[newtonsteps+1],
04641                              eta2[newtonsteps],
04642                             *M[M_index[number_]],
04643                              F2 );
04644 
04645             newtonsteps++;
04646         }
04647 
04648         for( run1 = 0; run1 < m; run1++ ){
04649             nablaY_(0,run1) = eta[nOfNewtonSteps[number_]][run1];
04650             nablaY_(1,run1) = nablaY_(0,run1) - nablaG2(0,run1);
04651             nablaY_(2,run1) = (nablaY_(1,run1) - nablaG2(1,run1)*scalT)*(psi[number_][0]/psi[number_][1]);
04652             nablaY_(3,run1) = (nablaY_(2,run1) - nablaG2(2,run1)*scalT*scalT)*(psi[number_][0]/psi[number_][2]);
04653             nablaY_(4,run1) = (nablaY_(3,run1) - nablaG2(3,run1)*scalT*scalT*scalT)*(psi[number_][0]/psi[number_][3]);
04654         }
04655         nablaG2 = nablaY_;
04656 
04657 
04658         for( run1 = 0; run1 < m; run1++ ){
04659             nablaY_(0,run1) = eta2[nOfNewtonSteps[number_]][run1];
04660             nablaY_(1,run1) = nablaY_(0,run1) - nablaG3(0,run1);
04661             nablaY_(2,run1) = (nablaY_(1,run1) - nablaG3(1,run1)*scalT)*(psi[number_][0]/psi[number_][1]);
04662             nablaY_(3,run1) = (nablaY_(2,run1) - nablaG3(2,run1)*scalT*scalT)*(psi[number_][0]/psi[number_][2]);
04663             nablaY_(4,run1) = (nablaY_(3,run1) - nablaG3(3,run1)*scalT*scalT*scalT)*(psi[number_][0]/psi[number_][3]);
04664         }
04665         nablaG3 = nablaY_;
04666     }
04667 }
04668 
04669 
04670 
04671 void IntegratorBDF::determineBDFEtaHBackward2( int number_ ){
04672 
04673     int run1, run2, run4;
04674     int newtonsteps;
04675     double pp;
04676 
04677     const double scalT = rel_time_scale/psi[number_][0];
04678 
04679     for( run4 = 0; run4 < nBDirs2; run4++ ){
04680 
04681         newtonsteps = nOfNewtonSteps[number_];
04682 
04683         for( run1 = 0; run1 < newtonsteps; run1++ ){
04684             for( run2 = 0; run2 < ndir; run2++ ){
04685                 l [run1][0][run2] = 0.0;
04686                 l2[run1][0][run2] = 0.0;
04687             }
04688         }
04689 
04690         for( run1 = 0; run1 < ndir; run1++ ){
04691 
04692             nablaH2_(4,run1) = nablaH2(4,run1)*scalT*scalT*scalT;
04693             nablaH2_(3,run1) = nablaH2(3,run1)*scalT*scalT + nablaH2_(4,run1)*(psi[number_][0]/psi[number_][3]);
04694             nablaH2_(2,run1) = nablaH2(2,run1)*scalT + nablaH2_(3,run1)*(psi[number_][0]/psi[number_][2]);
04695             nablaH2_(1,run1) = nablaH2(1,run1) + nablaH2_(2,run1)*(psi[number_][0]/psi[number_][1]);
04696             nablaH2_(0,run1) = nablaH2(0,run1) + nablaH2_(1,run1)/psi[number_][0];
04697 
04698             nablaH3_(4,run1) = nablaH3(4,run1)*scalT*scalT*scalT;
04699             nablaH3_(3,run1) = nablaH3(3,run1)*scalT*scalT + nablaH3_(4,run1)*(psi[number_][0]/psi[number_][3]);
04700             nablaH3_(2,run1) = nablaH3(2,run1)*scalT + nablaH3_(3,run1)*(psi[number_][0]/psi[number_][2]);
04701             nablaH3_(1,run1) = nablaH3(1,run1) + nablaH3_(2,run1)*(psi[number_][0]/psi[number_][1]);
04702             nablaH3_(0,run1) = nablaH3(0,run1) + nablaH3_(1,run1)/psi[number_][0];
04703 
04704             etaH2[newtonsteps][run1] = nablaH2_(0,run1);
04705             etaH3[newtonsteps][run1] = nablaH3_(0,run1);
04706             c2H2              (run1) = 0.0;
04707             c2H3              (run1) = 0.0;
04708         }
04709 
04710         newtonsteps--;
04711         while( newtonsteps >= 0 ){
04712 
04713             applyMTranspose( M_index[number_], etaH2[newtonsteps+1],
04714                             *M[M_index[number_]],
04715                              H2 );
04716 
04717             applyMTranspose( M_index[number_], etaH3[newtonsteps+1],
04718                             *M[M_index[number_]],
04719                              H3 );
04720 
04721             if( rhs[0].AD_backward2( 3*number_+newtonsteps, H2, H3,
04722                                      l[newtonsteps][0], l2[newtonsteps][0] )
04723                 != SUCCESSFUL_RETURN ){
04724                 ACADOERROR(RET_UNSUCCESSFUL_RETURN_FROM_INTEGRATOR_BDF);
04725             }
04726 
04727             for( run1 = 0; run1 < ndir; run1++ ){
04728                 etaH2[newtonsteps][run1] = etaH2[newtonsteps+1][run1] -
04729                                                  l [newtonsteps][0][run1];
04730                 etaH3[newtonsteps][run1] = etaH3[newtonsteps+1][run1] -
04731                                                  l2[newtonsteps][0][run1];
04732             }
04733 
04734             for( run1 = 0; run1 < md; run1++ ){
04735                 etaH2[newtonsteps][diff_index[run1]] =
04736                 etaH2[newtonsteps][diff_index[run1]] -
04737                 l[newtonsteps][0][ddiff_index[run1]]*gamma[number_][4];
04738                 etaH3[newtonsteps][diff_index[run1]] =
04739                 etaH3[newtonsteps][diff_index[run1]] -
04740                 l2[newtonsteps][0][ddiff_index[run1]]*gamma[number_][4];
04741             }
04742 
04743             for( run1 = 0; run1 < md; run1++ ){
04744                 c2H2(diff_index[run1]) -= l[newtonsteps][0][ddiff_index[run1]];
04745                 c2H3(diff_index[run1]) -= l2[newtonsteps][0][ddiff_index[run1]];
04746             }
04747 
04748             newtonsteps--;
04749         }
04750 
04751         for( run1 = 0; run1 < ndir; run1++ ){
04752 
04753             deltaH2(3,run1) = etaH2[0][run1]   - c2H2(run1)/psi[number_][3];
04754             deltaH2(2,run1) = deltaH2(3,run1) - c2H2(run1)/psi[number_][2];
04755             deltaH2(1,run1) = deltaH2(2,run1) - c2H2(run1)/psi[number_][1];
04756 
04757             deltaH3(3,run1) = etaH3[0][run1]   - c2H3(run1)/psi[number_][3];
04758             deltaH3(2,run1) = deltaH3(3,run1) - c2H3(run1)/psi[number_][2];
04759             deltaH3(1,run1) = deltaH3(2,run1) - c2H3(run1)/psi[number_][1];
04760         }
04761 
04762         for( run1 = 0; run1 < ndir; run1++ ){
04763             nablaH2(0,run1) =  - nablaH2_(1,run1) /psi[number_][0]
04764                                       - c2H2(run1) /psi[number_][0]
04765                                       + deltaH2(1,run1);
04766             nablaH3(0,run1) =  - nablaH3_(1,run1) /psi[number_][0]
04767                                       - c2H3(run1) /psi[number_][0]
04768                                       + deltaH3(1,run1);
04769         }
04770 
04771         pp = psi[number_][0];
04772         for( run1 = 0; run1 < ndir; run1++ ){
04773             nablaH2(1,run1) =  - nablaH2_(2,run1)*(psi[number_][0]/psi[number_][1])
04774                                       + pp*deltaH2(1,run1);
04775             nablaH3(1,run1) =  - nablaH3_(2,run1)*(psi[number_][0]/psi[number_][1])
04776                                       + pp*deltaH3(1,run1);
04777         }
04778         pp = pp*(psi[number_][1]/psi[number_][0]);
04779         for( run1 = 0; run1 < ndir; run1++ ){
04780             nablaH2(2,run1) =  - nablaH2_(3,run1)*(psi[number_][0]/psi[number_][2])
04781                                       + pp*deltaH2(2,run1);
04782             nablaH3(2,run1) =  - nablaH3_(3,run1)*(psi[number_][0]/psi[number_][2])
04783                                       + pp*deltaH3(2,run1);
04784         }
04785         pp = pp*(psi[number_][2]/psi[number_][0]);
04786         for( run1 = 0; run1 < ndir; run1++ ){
04787             nablaH2(3,run1) =  - nablaH2_(4,run1)*(psi[number_][0]/psi[number_][3])
04788                                       + pp*deltaH2(3,run1);
04789             nablaH3(3,run1) =  - nablaH3_(4,run1)*(psi[number_][0]/psi[number_][3])
04790                                       + pp*deltaH3(3,run1);
04791         }
04792         pp = pp*(psi[number_][3]/psi[number_][0]);
04793         for( run1 = 0; run1 < ndir; run1++ ){
04794             nablaH2(4,run1) =  pp*etaH2[0][run1];
04795             nablaH3(4,run1) =  pp*etaH3[0][run1];
04796         }
04797     }
04798 }
04799 
04800 
04801 
04802 void IntegratorBDF::initializeButcherTableau(){
04803 
04804 
04805     A[0][0] =  0.096   ;
04806     A[0][1] =  0.0     ;
04807     A[0][2] =  0.0     ;
04808     A[0][3] =  0.0     ;
04809     A[0][4] =  0.0     ;
04810     A[0][5] =  0.0     ;
04811     A[0][6] =  0.0     ;
04812 
04813     A[1][0] =  0.104   ;
04814     A[1][1] =  0.096   ;
04815     A[1][2] =  0.0     ;
04816     A[1][3] =  0.0     ;
04817     A[1][4] =  0.0     ;
04818     A[1][5] =  0.0     ;
04819     A[1][6] =  0.0     ;
04820 
04821     A[2][0] =  0.1310450349650284    ;
04822     A[2][1] =  0.07295496503497158102;
04823     A[2][2] =  0.096   ;
04824     A[2][3] =  0.0     ;
04825     A[2][4] =  0.0     ;
04826     A[2][5] =  0.0     ;
04827     A[2][6] =  0.0     ;
04828 
04829     A[3][0] =  0.2199597787833081951;
04830     A[3][1] = -0.1447179487179487179;
04831     A[3][2] =  0.02875816993464052288;
04832     A[3][3] =  0.096   ;
04833     A[3][4] =  0.0     ;
04834     A[3][5] =  0.0     ;
04835     A[3][6] =  0.0     ;
04836 
04837     A[4][0] =  0.1608848667672197084;
04838     A[4][1] = -0.0512400094214750;
04839     A[4][2] = -0.02467973856209150327;
04840     A[4][3] =  0.2190348812163467684;
04841     A[4][4] =  0.096   ;
04842     A[4][5] =  0.0     ;
04843     A[4][6] =  0.0     ;
04844 
04845     A[5][0] =  0.1761704161863276;
04846     A[5][1] = -0.2376951929172075;
04847     A[5][2] =  0.1249803878146932;
04848     A[5][3] =  0.3034259664066430296;
04849     A[5][4] =  0.1371184225095437181;
04850     A[5][5] =  0.096   ;
04851     A[5][6] =  0.0     ;
04852 
04853     A[6][0] =  0.1822523881347410759;
04854     A[6][1] = -0.3465441147470548;
04855     A[6][2] =  0.213542483660130719;
04856     A[6][3] =  0.3547492429521829614;
04857     A[6][4] =  0.1     ;
04858     A[6][5] =  0.2     ;
04859     A[6][6] =  0.096   ;
04860 
04861 
04862     b4[0] =  0.0;
04863     b4[1] =  0.0;
04864     b4[2] =  0.0;
04865     b4[3] =  0.4583333333333333333;
04866     b4[4] =  0.04166666666666666667;
04867     b4[5] =  0.04166666666666666667;
04868     b4[6] =  0.4583333333333333333;
04869 
04870     b5[0] = -0.3413924474339141;
04871     b5[1] =  0.8554210048117954;
04872     b5[2] = -0.5726796951403962;
04873     b5[3] =  0.4498353628579520;
04874     b5[4] =  0.0888157749045627;
04875     b5[5] =  0.0400000000000000;
04876     b5[6] =  0.4800000000000000;
04877 
04878     c[0] = 0.096;
04879     c[1] = 0.2;
04880     c[2] = 0.3;
04881     c[3] = 0.2;
04882     c[4] = 0.4;
04883     c[5] = 0.6;
04884     c[6] = 0.8;
04885 }
04886 
04887 
04888 
04889 void IntegratorBDF::printBDFfinalResults2( DMatrix &div ){
04890 
04891     int run2;
04892 
04893     cout <<"w.r.t. the states:\n" << scientific;
04894     for( run2 = 0; run2 < md; run2++ ){
04895         cout << div(0,diff_index[run2]) << "  ";
04896     }
04897     cout <<"\n";
04898 
04899     if( mu > 0 ){
04900         cout <<"w.r.t. the controls:\n" << scientific;
04901         for( run2 = 0; run2 < mu; run2++ ){
04902                 cout << div(0,control_index[run2]) << "  ";
04903         }
04904         cout <<"\n";
04905     }
04906     if( mp > 0 ){
04907         cout <<"w.r.t. the parameters:\n" << scientific;
04908         for( run2 = 0; run2 < mp; run2++ ){
04909                 cout << div(0,parameter_index[run2]) << "  ";
04910         }
04911         cout <<"\n";
04912     }
04913     if( mw > 0 ){
04914         cout <<"w.r.t. the diturbances:\n" << scientific;
04915         for( run2 = 0; run2 < mw; run2++ ){
04916                 cout << div(0,disturbance_index[run2]) << "  ";
04917         }
04918         cout <<"\n";
04919     }
04920 }
04921 
04922 
04923 
04924 void IntegratorBDF::printBDFfinalResults(){
04925 
04926         cout << scientific;
04927 
04928     int run1;
04929 
04930         // Forward Sensitivities:
04931         // ----------------------
04932 
04933         if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 ){
04934             cout <<"BDF: Forward Sensitivities:\n";
04935             for( run1 = 0; run1 < m; run1++ ){
04936 
04937                         cout << nablaG(0, run1);
04938             }
04939             cout <<"\n";
04940         }
04941 
04942         if( nFDirs2 > 0 ){
04943 
04944             cout <<"BDF: First Order Forward Sensitivities:\n";
04945             for( run1 = 0; run1 < m; run1++ ){
04946                 cout << nablaG2(0, run1);
04947             }
04948             cout <<"\n";
04949 
04950             cout <<"BDF: Second Order Forward Sensitivities:\n";
04951             for( run1 = 0; run1 < m; run1++ ){
04952                 cout << nablaG3(0, run1);
04953             }
04954             cout <<"\n";
04955         }
04956 
04957         // Backward Sensitivities:
04958         // -----------------------
04959 
04960         if( nBDirs > 0 ){
04961             cout <<"BDF: t = " << t-c[6]*h[0] << "  h = " <<  h[0] << endl;
04962             cout <<"BDF: Backward Sensitivities:\n";
04963             printBDFfinalResults2( nablaH );
04964         }
04965 
04966         // 2nd order Backward Sensitivities:
04967         // ---------------------------------
04968 
04969         if( nBDirs2 > 0 ){
04970                 cout <<"BDF: t = " << t-c[6]*h[0] << "  h = " <<  h[0] << endl;
04971 
04972             cout << "BDF: Backward Sensitivities:\n";
04973             printBDFfinalResults2( nablaH2 );
04974             cout <<"BDF: 2nd order Backward Sensitivities:\n";
04975             printBDFfinalResults2( nablaH3 );
04976         }
04977 }
04978 
04979 
04980 
04981 void IntegratorBDF::printBDFIntermediateResults(){
04982 
04983         cout << scientific;
04984 
04985     int run1;
04986 
04987      if( PrintLevel == HIGH ){
04988 
04989        if( nBDirs == 0 && nBDirs2 == 0 )
04990            cout <<"BDF: t = " << t-c[6]*h[0] << "  h = " <<  h[0] << endl;
04991 
04992        if( soa != SOA_EVERYTHING_FROZEN ){
04993 
04994          for( run1 = 0; run1 < md; run1++ ){
04995                  cout << "x[" << run1 << "] = " << nablaY(0,run1) << "  ";
04996          }
04997          for( run1 = 0; run1 < ma; run1++ ){
04998                  cout << "xa[" << run1 << "] = " << nablaY(0,md + run1) << "  ";
04999          }
05000          cout <<"\n";
05001        }
05002        else{
05003          if( nBDirs == 0 && nBDirs2 == 0 )
05004             cout <<"\n";
05005        }
05006 
05007         // Forward Sensitivities:
05008         // ----------------------
05009 
05010         if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 ){
05011             cout <<"BDF: Forward Sensitivities:\n";
05012             for( run1 = 0; run1 < m; run1++ ){
05013                 cout << nablaG(0,run1) << "  ";
05014             }
05015             cout <<"\n";
05016         }
05017 
05018         if( nFDirs2 > 0 ){
05019 
05020             cout <<"BDF: First Order Forward Sensitivities:\n";
05021             for( run1 = 0; run1 < m; run1++ ){
05022                 cout << nablaG2(0,run1) << "  ";
05023             }
05024             cout <<"\n";
05025 
05026             cout <<"BDF: Second Order Forward Sensitivities:\n";
05027             for( run1 = 0; run1 < m; run1++ ){
05028                 cout << nablaG3(0,run1) << "  ";
05029             }
05030             cout <<"\n";
05031         }
05032 
05033      }
05034 }
05035 
05036 
05037 void IntegratorBDF::printRKIntermediateResults(){
05038 
05039     int run1, run3;
05040 
05041     cout << scientific;
05042 
05043     if( PrintLevel == HIGH ){
05044 
05045         for( run3 = 0; run3 < 4; run3++ ){
05046             if( nBDirs == 0 && nBDirs2 == 0 )
05047             cout <<"BDF: t = " << t+h[0]*c[3+run3] << "  h = " << h[0]*c[3];
05048             if( soa != SOA_EVERYTHING_FROZEN ){
05049                 for( run1 = 0; run1 < md; run1++ ){
05050                                         cout << "x[" << run1 << "] = " << nablaY(3 - run3, run1) << "  ";
05051                 }
05052                 for( run1 = 0; run1 < ma; run1++ ){
05053                         cout << "xa[" << run1 << "] = " << nablaY(3-run3,run1+md) << "  ";
05054                 }
05055                 cout <<"\n";
05056             }
05057             else{
05058                 if( nBDirs == 0 && nBDirs2 == 0 )
05059                     cout <<"\n";
05060             }
05061 
05062             // Forward Sensitivities:
05063             // ----------------------
05064 
05065             if( nFDirs > 0 && nBDirs2 == 0 && nFDirs2 == 0 ){
05066                 cout <<"BDF: Forward Sensitivities:\n";
05067                 for( run1 = 0; run1 < m; run1++ ){
05068                         cout << nablaG(3-run3,run1) << "  ";
05069                 }
05070                 cout <<"\n";
05071             }
05072 
05073             if( nFDirs2 > 0 ){
05074 
05075                 cout <<"BDF: First Order Forward Sensitivities:\n";
05076                 for( run1 = 0; run1 < m; run1++ ){
05077                     cout << nablaG2(3-run3,run1) << "  ";
05078                 }
05079                 cout <<"\n";
05080 
05081                 cout <<"BDF: Second Order Forward Sensitivities:\n";
05082                 for( run1 = 0; run1 < m; run1++ ){
05083                     cout << nablaG3(3-run3,run1) << "  ";
05084                 }
05085                 cout <<"\n";
05086             }
05087         }
05088 
05089         // Backward Sensitivities:
05090         // -----------------------
05091 
05092         if( nBDirs > 0 ){
05093             cout << scientific
05094                  << "BDF: t = " << t-c[6]*h[0] << "  h_RK = " << h[0] << endl
05095                  << "BDF: Backward Sensitivities:\n";
05096             printBDFfinalResults2( nablaH );
05097         }
05098 
05099         // 2nd order Backward Sensitivities:
05100         // ---------------------------------
05101 
05102         if( nBDirs2 > 0 ){
05103             cout << scientific
05104                  << "BDF: t = " << t-c[6]*h[0] << "  " << "h_RK = " << h[0] << endl
05105                  << "BDF: Backward Sensitivities:\n";
05106             printBDFfinalResults2( nablaH2 );
05107             cout << "BDF: 2nd order Backward Sensitivities:\n";
05108             printBDFfinalResults2( nablaH3 );
05109         }
05110     }
05111 }
05112 
05113 
05114 returnValue IntegratorBDF::decomposeJacobian(int index, DMatrix &J){
05115 
05116 //      ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05117 
05118     switch( las ){
05119 
05120         case HOUSEHOLDER_METHOD:
05121 //              ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05122                 break;
05123 //              ASSERT( index < qr.size() );
05124 //              qr[ index ] = Eigen::HouseholderQR< DMatrix::Base >( J );
05125 //              return SUCCESSFUL_RETURN;
05126 
05127         case SPARSE_LU:
05128                         ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05129                         break;
05130 //             return J.computeSparseLUdecomposition();
05131 
05132         default:
05133              return ACADOERROR( RET_NOT_IMPLEMENTED_YET );
05134     }
05135 
05136     return SUCCESSFUL_RETURN;
05137 }
05138 
05139 
05140 double IntegratorBDF::applyNewtonStep( int index, double *etakplus1, const double *etak, const DMatrix &J, const double *FFF ){
05141 
05142     int run1;
05143     DVector bb(m,FFF);
05144     DVector deltaX;
05145 
05146         switch (las)
05147         {
05148         case HOUSEHOLDER_METHOD:
05149                 deltaX = J.householderQr().solve( bb );
05150 //              deltaX = qr[ index ].solve(bb);
05151                 break;
05152         case SPARSE_LU:
05153                 ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05154 //              deltaX = J.solveSparseLU(bb);
05155                 break;
05156         default:
05157                 deltaX.setZero();
05158                 break;
05159         }
05160 
05161     for( run1 = 0; run1 < m; run1++ )
05162         etakplus1[run1] = etak[run1] - deltaX(run1);
05163 
05164     return deltaX.getNorm( VN_LINF, diff_scale );
05165 }
05166 
05167 
05168 void IntegratorBDF::applyMTranspose( int index, double *seed1, DMatrix &J, double *seed2 ){
05169 
05170     int run1;
05171     DVector bb(m);
05172 
05173     for( run1 = 0; run1 < m; run1++ )
05174         bb(run1) = seed1[diff_index[run1]];
05175 
05176     DVector deltaX;
05177 
05178         switch (las)
05179         {
05180         case HOUSEHOLDER_METHOD:
05181                 // Dodgy
05182                 deltaX = J.transpose().householderQr().solve( bb );
05183 //              deltaX = J.solveTransposeQR(bb);
05184 //              deltaX = qr[ index ].matrixQR().block(0, 0, J.cols(), J.cols()).
05185 //                      triangularView<Eigen::Upper>().transpose().solve( bb );
05186                 break;
05187         case SPARSE_LU:
05188                 ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05189 //              deltaX = J.solveTransposeSparseLU(bb);
05190                 break;
05191         default:
05192                 ACADOFATAL(  RET_NOT_IMPLEMENTED_YET );
05193                 deltaX.setZero();
05194                 break;
05195         }
05196 
05197     for( run1 = 0; run1 < m; run1++ )
05198         seed2[run1] = deltaX(run1);
05199 }
05200 
05201 
05202 
05203 void IntegratorBDF::relaxAlgebraic( double *residuum, double timePoint ){
05204 
05205     int           relaxationType  ;
05206     double        relaxationPar   ;
05207 
05208     get( RELAXATION_PARAMETER, relaxationPar );
05209 
05210     const double  a = relaxationPar*(timeInterval.getIntervalLength());
05211     const double  b = 1.0;
05212     double        damping = 1.0;
05213     int           run1   ;
05214     double        normRES;
05215 
05216     get( ALGEBRAIC_RELAXATION, relaxationType );
05217 
05218     switch( relaxationType ){
05219 
05220         case ART_EXPONENTIAL:
05221              damping = exp(relaxationConstant*(timeInterval.getFirstTime()-timePoint)/( timeInterval.getIntervalLength() ));
05222              for( run1 = 0; run1 < ma; run1++ )
05223                  residuum[md+run1] -= initialAlgebraicResiduum[md+run1]*damping;
05224              break;
05225 
05226         case ART_ADAPTIVE_POLYNOMIAL:
05227              normRES = 0.0;
05228              for( run1 = 0; run1 < ma; run1++ )
05229                  if( fabs(initialAlgebraicResiduum[md+run1]) > normRES )
05230                      normRES = fabs(initialAlgebraicResiduum[md+run1]);
05231 
05232              damping = pow( 1.0 - ( timePoint-timeInterval.getFirstTime())/( timeInterval.getIntervalLength() ) ,
05233                                     b + a/(normRES+100.0*EPS) );
05234 
05235              break;
05236 
05237         case ART_UNKNOWN:
05238              damping = 1.0;
05239              break;
05240     }
05241 
05242     for( run1 = 0; run1 < ma; run1++ )
05243         residuum[md+run1] -= initialAlgebraicResiduum[md+run1]*damping;
05244 
05245     return;
05246 }
05247 
05248 
05249 
05250 int IntegratorBDF::getDim() const{
05251 
05252     return md+ma;
05253 }
05254 
05255 int IntegratorBDF::getDimX() const{
05256 
05257     return md;
05258 }
05259 
05260 
05261 void IntegratorBDF::prepareDividedDifferences( DMatrix &div ){
05262 
05263     int run1;
05264 
05265     for( run1 = 0; run1 < m; run1++ ){
05266 
05267         div(3,run1) = div(2,run1) - div(3,run1);
05268 
05269         div(2,run1) = div(1,run1) - div(2,run1);
05270         div(3,run1) = ( div(2,run1) - div(3,run1) )/(2.0);
05271 
05272         div(1,run1) = div(0,run1) - div(1,run1);
05273         div(2,run1) = ( div(1,run1) - div(2,run1) )/(2.0);
05274         div(3,run1) = ( div(2,run1) - div(3,run1) )/(3.0);
05275 
05276         div(4,run1) = 0.0; // the RK-starter has order 3 only.
05277     }
05278     return;
05279 }
05280 
05281 
05282 void IntegratorBDF::copyBackward( DVector       &Dx_x0,
05283                                   DVector       &Dx_p ,
05284                                   DVector       &Dx_u ,
05285                                   DVector       &Dx_w ,
05286                                   const DMatrix &div    ) const{
05287 
05288     int run1;
05289 
05290     if( Dx_x0.getDim() != 0 )
05291         for( run1 = 0; run1 < m; run1++ )
05292             Dx_x0(run1)= div(0,diff_index[run1]);
05293 
05294     if( Dx_u.getDim() != 0 )
05295         for( run1 = 0; run1 < mu; run1++ )
05296             Dx_u(run1) = div(0,control_index[run1]);
05297 
05298     if( Dx_p.getDim() != 0 )
05299         for( run1 = 0; run1 < mp; run1++ )
05300             Dx_p(run1) = div(0,parameter_index[run1]);
05301 
05302     if( Dx_w.getDim() != 0 )
05303         for( run1 = 0; run1 < mw; run1++ )
05304             Dx_w(run1) = div(0,disturbance_index[run1]);
05305 
05306     return;
05307 }
05308 
05309 
05310 void IntegratorBDF::interpolate( int number_, DMatrix &div, VariablesGrid &poly ){
05311 
05312     int i1 = timeInterval.getFloorIndex( t-h[0] );
05313     int i2 = timeInterval.getFloorIndex( t      );
05314     int jj, run1;
05315 
05316     for( jj = i1+1; jj <= i2; jj++ ){
05317 
05318         for( run1 = 0; run1 < m; run1++ )
05319             poly( jj, run1 ) = div(0,run1);
05320 
05321         double pp1 = (timeInterval.getTime(jj) - t)/h[0];
05322         double pp2 = pp1*(pp1*h[0] + h[0])/h[0];
05323         double pp3 = pp2*(pp1*h[0] + psi[number_][1])/h[0];
05324         double pp4 = pp3*(pp1*h[0] + psi[number_][2])/h[0];
05325 
05326         for( run1 = 0; run1 < m; run1++ )
05327             poly( jj, run1 ) += pp1*div(1,run1);
05328 
05329         for( run1 = 0; run1 < m; run1++ )
05330             poly( jj, run1 ) += pp2*div(2,run1);
05331 
05332         for( run1 = 0; run1 < m; run1++ )
05333             poly( jj, run1 ) += pp3*div(3,run1);
05334 
05335         for( run1 = 0; run1 < m; run1++ )
05336             poly( jj, run1 ) += pp4*div(4,run1);
05337     }
05338 }
05339 
05340 
05341 void IntegratorBDF::logCurrentIntegratorStep(   const DVector& currentX,
05342                                                                                                 const DVector& currentXA
05343                                                                                                 )
05344 {
05345         return;
05346         
05347         int run1;
05348 
05349         // log differential states
05350         if ( currentX.isEmpty( ) == BT_TRUE )
05351         {
05352                 DVector currentDiffStates(md);
05353                 for( run1 = 0; run1 < md; run1++ )
05354                         currentDiffStates( run1 ) = nablaY( 0,run1 );
05355                 setLast( LOG_DIFFERENTIAL_STATES,currentDiffStates,t );
05356         }
05357         else
05358         {
05359                 setLast( LOG_DIFFERENTIAL_STATES,currentX,t );
05360         }
05361 
05362         // log algebraic states
05363         if ( currentX.isEmpty( ) == BT_TRUE )
05364         {
05365                 DVector currentAlgStates(ma);
05366                 for( run1 = 0; run1 < ma; run1++ )
05367                         currentAlgStates( run1 ) = nablaY( 0,md+run1 );
05368                 setLast( LOG_ALGEBRAIC_STATES,currentAlgStates,t );
05369         }
05370         else
05371         {
05372                 setLast( LOG_ALGEBRAIC_STATES,currentXA,t );
05373         }
05374 }
05375 
05376 
05377 CLOSE_NAMESPACE_ACADO
05378 
05379 // end of file.


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:27