[2696] | 1 | #define MAX(a,b) ( ((a) >= (b)) ? (a):(b) ) |
---|
| 2 | #define MIN(b,c) ( ((b) < (c)) ? (b):(c) ) |
---|
| 3 | #define ABS(x) ( ((x) >= 0 ) ? (x):(-x) ) |
---|
| 4 | #define SQRT(d) ( pow((d),0.5) ) |
---|
| 5 | |
---|
| 6 | /* Numerical Constants */ |
---|
| 7 | #define ZERO (KPP_REAL)0.0 |
---|
| 8 | #define ONE (KPP_REAL)1.0 |
---|
| 9 | #define HALF (KPP_REAL)0.5 |
---|
| 10 | #define DeltaMin (KPP_REAL)1.0e-5 |
---|
| 11 | enum boolean { FALSE=0, TRUE=1 }; |
---|
| 12 | |
---|
| 13 | /* Statistics on the work performed by the Rosenbrock method */ |
---|
| 14 | enum statistics { Nfun=1, Njac=2, Nstp=3, Nacc=4, Nrej=5, Ndec=6, Nsol=7, |
---|
| 15 | Nsng=8, Ntexit=1, Nhexit=2, Nhnew=3 }; |
---|
| 16 | |
---|
| 17 | /*~~~> Parameters of the Rosenbrock method, up to 6 stages */ |
---|
| 18 | int ros_S, rosMethod; |
---|
| 19 | enum ros_Params { RS2=1, RS3=2, RS4=3, RD3=4, RD4=5 }; |
---|
| 20 | KPP_REAL ros_A[15], ros_C[15], ros_M[6], ros_E[6], ros_Alpha[6], ros_Gamma[6], |
---|
| 21 | ros_ELO; |
---|
| 22 | int ros_NewF[6]; /* Holds Boolean values */ |
---|
| 23 | char ros_Name[12]; /* Length 12 */ |
---|
| 24 | |
---|
| 25 | /*~~~> Types of Adjoints Implemented */ |
---|
| 26 | enum adjoint { Adj_none=1, Adj_discrete=2, Adj_continuous=3, |
---|
| 27 | Adj_simple_continuous=4 }; |
---|
| 28 | |
---|
| 29 | /*~~~> Checkpoints in memory */ |
---|
| 30 | int bufsize = 200000; |
---|
| 31 | int stack_ptr; /* last written entry */ |
---|
| 32 | KPP_REAL *chk_H, *chk_T; |
---|
| 33 | KPP_REAL **chk_Y, **chk_K, **chk_J; /* 2D arrays */ |
---|
| 34 | KPP_REAL **chk_dY, **chk_d2Y; /* 2D arrays */ |
---|
| 35 | |
---|
| 36 | /* Function Headers */ |
---|
| 37 | void INTEGRATE_ADJ(int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], |
---|
| 38 | KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], |
---|
| 39 | KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], |
---|
| 40 | KPP_REAL RCNTRL_U[], int ISTATUS_U[], KPP_REAL RSTATUS_U[]); |
---|
| 41 | int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR], |
---|
| 42 | KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], |
---|
| 43 | KPP_REAL RelTol[], KPP_REAL AbsTol_adj[][NVAR], |
---|
| 44 | KPP_REAL RelTol_adj[][NVAR], KPP_REAL RCNTRL[], |
---|
| 45 | int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[] ); |
---|
| 46 | void ros_AllocateDBuffers( int S, int SaveLU ); |
---|
| 47 | void ros_FreeDBuffers( int SaveLU ); |
---|
| 48 | void ros_AllocateCBuffers(); |
---|
| 49 | void ros_FreeCBuffers(); |
---|
| 50 | void ros_DPush( int S, KPP_REAL T, KPP_REAL H, KPP_REAL Ystage[], |
---|
| 51 | KPP_REAL K[], KPP_REAL E[], int P[], int SaveLU ); |
---|
| 52 | void ros_DPop( int S, KPP_REAL* T, KPP_REAL* H, KPP_REAL* Ystage, |
---|
| 53 | KPP_REAL* K, KPP_REAL* E, int* P, int SaveLU ); |
---|
| 54 | void ros_CPush( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], |
---|
| 55 | KPP_REAL d2Y[] ); |
---|
| 56 | void ros_CPop( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], |
---|
| 57 | KPP_REAL d2Y[] ); |
---|
| 58 | int ros_ErrorMsg( int Code, KPP_REAL T, KPP_REAL H); |
---|
| 59 | int ros_FwdInt (KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL T, |
---|
| 60 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int AdjointType, |
---|
| 61 | KPP_REAL Hmin, KPP_REAL Hstart, KPP_REAL Hmax, |
---|
| 62 | KPP_REAL Roundoff, int ISTATUS[], int Max_no_steps, |
---|
| 63 | KPP_REAL RSTATUS[], int Autonomous, int VectorTol, |
---|
| 64 | KPP_REAL FacMax, KPP_REAL FacMin, KPP_REAL FacSafe, |
---|
| 65 | KPP_REAL FacRej, int SaveLU); |
---|
| 66 | int ros_DadjInt ( int NADJ, KPP_REAL Lambda[][NVAR], KPP_REAL Tstart, |
---|
| 67 | KPP_REAL Tend, KPP_REAL T, int SaveLU, int ISTATUS[], |
---|
| 68 | KPP_REAL Roundoff, int Autonomous); |
---|
| 69 | int ros_CadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, KPP_REAL Tend, |
---|
| 70 | KPP_REAL T, KPP_REAL AbsTol_adj[][NVAR], |
---|
| 71 | KPP_REAL RelTol_adj[][NVAR], KPP_REAL RSTATUS[], |
---|
| 72 | KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, |
---|
| 73 | KPP_REAL Roundoff, int Max_no_steps, int Autonomous, |
---|
| 74 | int VectorTol, KPP_REAL FacMax, KPP_REAL FacMin, |
---|
| 75 | KPP_REAL FacSafe, KPP_REAL FacRej, int ISTATUS[] ); |
---|
| 76 | int ros_SimpleCadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, |
---|
| 77 | KPP_REAL Tend, KPP_REAL T, int ISTATUS[], |
---|
| 78 | int Autonomous, KPP_REAL Roundoff ); |
---|
| 79 | KPP_REAL ros_ErrorNorm ( KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], |
---|
| 80 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int VectorTol ); |
---|
| 81 | void ros_FunTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], |
---|
| 82 | KPP_REAL Fcn0[], KPP_REAL dFdT[], int ISTATUS[] ); |
---|
| 83 | void ros_JacTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], |
---|
| 84 | KPP_REAL Jac0[], KPP_REAL dJdT[], int ISTATUS[] ); |
---|
| 85 | int ros_PrepareMatrix ( KPP_REAL H, int Direction, KPP_REAL gam, |
---|
| 86 | KPP_REAL Jac0[], KPP_REAL Ghimj[], int Pivot[], |
---|
| 87 | int ISTATUS[] ); |
---|
| 88 | void ros_Decomp( KPP_REAL A[], int Pivot[], int* ising, int ISTATUS[] ); |
---|
| 89 | void ros_Solve( char How, KPP_REAL A[], int Pivot[], KPP_REAL b[], |
---|
| 90 | int ISTATUS[] ); |
---|
| 91 | void ros_cadj_Y( KPP_REAL T, KPP_REAL Y[] ); |
---|
| 92 | void ros_Hermite3( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], |
---|
| 93 | KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Y[] ); |
---|
| 94 | void ros_Hermite5( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], |
---|
| 95 | KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Ha[], |
---|
| 96 | KPP_REAL Hb[], KPP_REAL Y[] ); |
---|
| 97 | void Ros2(); |
---|
| 98 | void Ros3(); |
---|
| 99 | void Ros4(); |
---|
| 100 | void Rodas3(); |
---|
| 101 | void Rodas4(); |
---|
| 102 | void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ); |
---|
| 103 | void HessTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Hes[] ); |
---|
| 104 | void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Fun [] ); |
---|
| 105 | void WSCAL( int N, KPP_REAL Alpha, KPP_REAL X[], int incX ); |
---|
| 106 | void WAXPY( int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], |
---|
| 107 | int incY ); |
---|
| 108 | void WCOPY( int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY ); |
---|
| 109 | KPP_REAL WLAMCH( char C ); |
---|
| 110 | void Update_SUN(); |
---|
| 111 | void Update_RCONST(); |
---|
| 112 | void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] ); |
---|
| 113 | void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]); |
---|
| 114 | void Jac_SP_Vec( KPP_REAL Jac[], KPP_REAL Fcn[], KPP_REAL K[] ); |
---|
| 115 | void JacTR_SP_Vec( KPP_REAL Jac[], KPP_REAL Fcn[], KPP_REAL K[] ); |
---|
| 116 | void HessTR_Vec(KPP_REAL Hess[], KPP_REAL U1[], KPP_REAL U2[], KPP_REAL HTU[]); |
---|
| 117 | void KppSolve( KPP_REAL A[], KPP_REAL b[] ); |
---|
| 118 | int KppDecomp( KPP_REAL A[] ); |
---|
| 119 | void KppSolveTR( KPP_REAL JVS[], KPP_REAL X[], KPP_REAL XX[] ); |
---|
| 120 | void Hessian( KPP_REAL V[], KPP_REAL F[], KPP_REAL RCT[], KPP_REAL Hess[] ); |
---|
| 121 | |
---|
| 122 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 123 | void INTEGRATE_ADJ( int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], |
---|
| 124 | KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], |
---|
| 125 | KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], |
---|
| 126 | KPP_REAL RCNTRL_U[], int ISTATUS_U[], |
---|
| 127 | KPP_REAL RSTATUS_U[] ) { |
---|
| 128 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 129 | |
---|
| 130 | /* Local Variables */ |
---|
| 131 | KPP_REAL RCNTRL[20], RSTATUS[20]; |
---|
| 132 | int ICNTRL[20], ISTATUS[20], IERR, i; |
---|
| 133 | |
---|
| 134 | for( i = 0; i < 20; i++ ) { |
---|
| 135 | ICNTRL[i] = 0; |
---|
| 136 | RCNTRL[i] = ZERO; |
---|
| 137 | ISTATUS[i] = 0; |
---|
| 138 | RSTATUS[i] = ZERO; |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | /*~~~> fine-tune the integrator: |
---|
| 142 | ICNTRL(1) = 0 ! 0 = non-autonomous, 1 = autonomous |
---|
| 143 | ICNTRL(2) = 1 ! 0 = scalar, 1 = vector tolerances |
---|
| 144 | RCNTRL(3) = STEPMIN ! starting step |
---|
| 145 | ICNTRL(3) = 5 ! choice of the method for forward integration |
---|
| 146 | ICNTRL(6) = 1 ! choice of the method for continuous adjoint |
---|
| 147 | ICNTRL(7) = 2 ! 1=none, 2=discrete, 3=full continuous, |
---|
| 148 | 4=simplified continuous adjoint |
---|
| 149 | ICNTRL(8) = 1 ! Save fwd LU factorization: 0=*don't* save, 1=save */ |
---|
| 150 | |
---|
| 151 | /* if optional parameters are given, and if they are >=0, then they overwrite |
---|
| 152 | default settings */ |
---|
| 153 | // if(ICNTRL_U != NULL) { |
---|
| 154 | // for(i=0; i<20; i++) |
---|
| 155 | // if(ICNTRL_U[i] > 0) |
---|
| 156 | // ICNTRL[i] = ICNTRL_U[i]; |
---|
| 157 | // } /* end for */ |
---|
| 158 | // } /* end if */ |
---|
| 159 | |
---|
| 160 | // if(RCNTRL_U != NULL) { |
---|
| 161 | // for(i=0; i<20; i++) |
---|
| 162 | // if(RCNTRL_U[i] > 0) |
---|
| 163 | // RCNTRL[i] = RCNTRL_U[i]; |
---|
| 164 | // } /* end for */ |
---|
| 165 | // } /* end if */ |
---|
| 166 | |
---|
| 167 | IERR = RosenbrockADJ( Y, NADJ, Lambda, TIN, TOUT, ATOL, RTOL, ATOL_adj, |
---|
| 168 | RTOL_adj, RCNTRL, ICNTRL, RSTATUS, ISTATUS ); |
---|
| 169 | |
---|
| 170 | if (IERR < 0) |
---|
| 171 | printf( "RosenbrockADJ: Unsucessful step at T=%f (IERR=%d)", TIN, IERR ); |
---|
| 172 | |
---|
| 173 | STEPMIN = RSTATUS[Nhexit]; |
---|
| 174 | |
---|
| 175 | /* if optional parameters are given for output |
---|
| 176 | copy to them to return information */ |
---|
| 177 | // if(ISTATUS_U != NULL) |
---|
| 178 | // for(i=0; i<20; i++) |
---|
| 179 | // ISTATUS_U[i] = ISTATUS[i]; |
---|
| 180 | // } |
---|
| 181 | |
---|
| 182 | // if(RSTATUS_U != NULL) |
---|
| 183 | // for(i=0; i<20; i++) |
---|
| 184 | // RSTATUS_U[i] = RSTATUS[i]; |
---|
| 185 | // } |
---|
| 186 | |
---|
| 187 | } /* End of INTEGRATE_ADJ */ |
---|
| 188 | |
---|
| 189 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 190 | int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR], |
---|
| 191 | KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], |
---|
| 192 | KPP_REAL RelTol[], KPP_REAL AbsTol_adj[][NVAR], |
---|
| 193 | KPP_REAL RelTol_adj[][NVAR], KPP_REAL RCNTRL[], |
---|
| 194 | int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[] ) { |
---|
| 195 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 196 | |
---|
| 197 | ADJ = Adjoint of the Tangent Linear Model of a Rosenbrock Method |
---|
| 198 | |
---|
| 199 | Solves the system y'=F(t,y) using a RosenbrockADJ method defined by: |
---|
| 200 | |
---|
| 201 | G = 1/(H*gamma(1)) - Jac(t0,Y0) |
---|
| 202 | T_i = t0 + Alpha(i)*H |
---|
| 203 | Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j |
---|
| 204 | G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j + |
---|
| 205 | gamma(i)*dF/dT(t0, Y0) |
---|
| 206 | Y1 = Y0 + \sum_{j=1}^S M(j)*K_j |
---|
| 207 | |
---|
| 208 | For details on RosenbrockADJ methods and their implementation consult: |
---|
| 209 | E. Hairer and G. Wanner |
---|
| 210 | "Solving ODEs II. Stiff and differential-algebraic problems". |
---|
| 211 | Springer series in computational mathematics, Springer-Verlag, 1996. |
---|
| 212 | The codes contained in the book inspired this implementation. |
---|
| 213 | |
---|
| 214 | (C) Adrian Sandu, August 2004 |
---|
| 215 | Virginia Polytechnic Institute and State University |
---|
| 216 | Contact: sandu@cs.vt.edu |
---|
| 217 | Revised by Philipp Miehe and Adrian Sandu, May 2006 |
---|
| 218 | Translation F90 to C by Paul Eller, April 2007 |
---|
| 219 | This implementation is part of KPP - the Kinetic PreProcessor |
---|
| 220 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 221 | |
---|
| 222 | ~~~> INPUT ARGUMENTS: |
---|
| 223 | |
---|
| 224 | - Y[NVAR] = vector of initial conditions (at T=Tstart) |
---|
| 225 | NADJ -> dimension of linearized system, |
---|
| 226 | i.e. the number of sensitivity coefficients |
---|
| 227 | - Lambda[NVAR][NADJ] -> vector of initial sensitivity conditions |
---|
| 228 | (at T=Tstart) |
---|
| 229 | - [Tstart,Tend] = time range of integration |
---|
| 230 | (if Tstart>Tend the integration is performed backwards in time) |
---|
| 231 | - RelTol, AbsTol = user precribed accuracy |
---|
| 232 | - void Fun( T, Y, Ydot ) = ODE function, |
---|
| 233 | returns Ydot = Y' = F(T,Y) |
---|
| 234 | - void Jac( T, Y, Jcb ) = Jacobian of the ODE function, |
---|
| 235 | returns Jcb = dF/dY |
---|
| 236 | - ICNTRL[0:9] = integer inputs parameters |
---|
| 237 | - RCNTRL[0:9] = real inputs parameters |
---|
| 238 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 239 | |
---|
| 240 | ~~~> OUTPUT ARGUMENTS: |
---|
| 241 | |
---|
| 242 | - Y[NVAR] -> vector of final states (at T->Tend) |
---|
| 243 | - Lambda[NVAR][NADJ] -> vector of final sensitivities (at T=Tend) |
---|
| 244 | - ICNTRL[9:18] -> integer output parameters |
---|
| 245 | - RCNTRL[9:18] -> real output parameters |
---|
| 246 | - IERR -> job status upon return |
---|
| 247 | - succes (positive value) or failure (negative value) - |
---|
| 248 | = 1 : Success |
---|
| 249 | = -1 : Improper value for maximal no of steps |
---|
| 250 | = -2 : Selected RosenbrockADJ method not implemented |
---|
| 251 | = -3 : Hmin/Hmax/Hstart must be positive |
---|
| 252 | = -4 : FacMin/FacMax/FacRej must be positive |
---|
| 253 | = -5 : Improper tolerance values |
---|
| 254 | = -6 : No of steps exceeds maximum bound |
---|
| 255 | = -7 : Step size too small |
---|
| 256 | = -8 : Matrix is repeatedly singular |
---|
| 257 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 258 | |
---|
| 259 | ~~~> INPUT PARAMETERS: |
---|
| 260 | |
---|
| 261 | Note: For input parameters equal to zero the default values of the |
---|
| 262 | corresponding variables are used. |
---|
| 263 | |
---|
| 264 | ICNTRL[0] = 1: F = F(y) Independent of T (AUTONOMOUS) |
---|
| 265 | = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS) |
---|
| 266 | |
---|
| 267 | ICNTRL[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors |
---|
| 268 | = 1: AbsTol, RelTol are scalars |
---|
| 269 | |
---|
| 270 | ICNTRL[2] -> selection of a particular Rosenbrock method |
---|
| 271 | = 0 : default method is Rodas3 |
---|
| 272 | = 1 : method is Ros2 |
---|
| 273 | = 2 : method is Ros3 |
---|
| 274 | = 3 : method is Ros4 |
---|
| 275 | = 4 : method is Rodas3 |
---|
| 276 | = 5: method is Rodas4 |
---|
| 277 | |
---|
| 278 | ICNTRL[3] -> maximum number of integration steps |
---|
| 279 | For ICNTRL[4]=0) the default value of BUFSIZE is used |
---|
| 280 | |
---|
| 281 | ICNTRL[5] -> selection of a particular Rosenbrock method for the |
---|
| 282 | continuous adjoint integration - for cts adjoint it |
---|
| 283 | can be different than the forward method ICNTRL(3) |
---|
| 284 | Note 1: to avoid interpolation errors (which can be huge!) |
---|
| 285 | it is recommended to use only ICNTRL[6] = 2 or 4 |
---|
| 286 | Note 2: the performance of the full continuous adjoint |
---|
| 287 | strongly depends on the forward solution accuracy Abs/RelTol |
---|
| 288 | |
---|
| 289 | ICNTRL[6] -> Type of adjoint algorithm |
---|
| 290 | = 0 : default is discrete adjoint ( of method ICNTRL[3] ) |
---|
| 291 | = 1 : no adjoint |
---|
| 292 | = 2 : discrete adjoint ( of method ICNTRL[3] ) |
---|
| 293 | = 3 : fully adaptive continuous adjoint ( with method ICNTRL[6] ) |
---|
| 294 | = 4 : simplified continuous adjoint ( with method ICNTRL[6] ) |
---|
| 295 | |
---|
| 296 | ICNTRL[7] -> checkpointing the LU factorization at each step: |
---|
| 297 | ICNTRL[7]=0 : do *not* save LU factorization (the default) |
---|
| 298 | ICNTRL[7]=1 : save LU factorization |
---|
| 299 | Note: if ICNTRL[7]=1 the LU factorization is *not* saved |
---|
| 300 | |
---|
| 301 | ~~~> Real input parameters: |
---|
| 302 | |
---|
| 303 | RCNTRL[0] -> Hmin, lower bound for the integration step size |
---|
| 304 | It is strongly recommended to keep Hmin = ZERO |
---|
| 305 | |
---|
| 306 | RCNTRL[1] -> Hmax, upper bound for the integration step size |
---|
| 307 | |
---|
| 308 | RCNTRL[2] -> Hstart, starting value for the integration step size |
---|
| 309 | |
---|
| 310 | RCNTRL[3] -> FacMin, lower bound on step decrease factor (default=0.2) |
---|
| 311 | |
---|
| 312 | RCNTRL[4] -> FacMax, upper bound on step increase factor (default=6) |
---|
| 313 | |
---|
| 314 | RCNTRL[5] -> FacRej, step decrease factor after multiple rejections |
---|
| 315 | (default=0.1) |
---|
| 316 | |
---|
| 317 | RCNTRL[6] -> FacSafe, by which the new step is slightly smaller |
---|
| 318 | than the predicted value (default=0.9) |
---|
| 319 | |
---|
| 320 | RCNTRL[7] -> ThetaMin. If Newton convergence rate smaller |
---|
| 321 | than ThetaMin the Jacobian is not recomputed; |
---|
| 322 | (default=0.001) |
---|
| 323 | |
---|
| 324 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 325 | |
---|
| 326 | ~~~> OUTPUT PARAMETERS: |
---|
| 327 | |
---|
| 328 | Note: each call to RosenbrockADJ adds the corrent no. of fcn calls |
---|
| 329 | to previous value of ISTATUS(1), and similar for the other params. |
---|
| 330 | Set ISTATUS[0:9] = 0 before call to avoid this accumulation. |
---|
| 331 | |
---|
| 332 | ISTATUS[0] = No. of function calls |
---|
| 333 | ISTATUS[1] = No. of jacobian calls |
---|
| 334 | ISTATUS[2] = No. of steps |
---|
| 335 | ISTATUS[3] = No. of accepted steps |
---|
| 336 | ISTATUS[4] = No. of rejected steps (except at the beginning) |
---|
| 337 | ISTATUS[5] = No. of LU decompositions |
---|
| 338 | ISTATUS[6] = No. of forward/backward substitutions |
---|
| 339 | ISTATUS[7] = No. of singular matrix decompositions |
---|
| 340 | |
---|
| 341 | RSTATUS[0] -> Texit, the time corresponding to the |
---|
| 342 | computed Y upon return |
---|
| 343 | RSTATUS[1] -> Hexit, last accepted step before exit |
---|
| 344 | For multiple restarts, use Hexit as Hstart in the following run |
---|
| 345 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 346 | |
---|
| 347 | /*~~~> Local variables */ |
---|
| 348 | KPP_REAL Roundoff, FacMin, FacMax, FacRej, FacSafe; |
---|
| 349 | KPP_REAL Hmin, Hmax, Hstart; |
---|
| 350 | KPP_REAL Texit=0.0; |
---|
| 351 | int i, UplimTol, Max_no_steps=0, IERR; |
---|
| 352 | int AdjointType=0, CadjMethod=0; |
---|
| 353 | int Autonomous, VectorTol, SaveLU; /* Holds boolean values */ |
---|
| 354 | |
---|
| 355 | stack_ptr = -1; |
---|
| 356 | |
---|
| 357 | /*~~~> Initialize statistics */ |
---|
| 358 | for(i=0; i<20; i++) { |
---|
| 359 | ISTATUS[i] = 0; |
---|
| 360 | RSTATUS[i] = ZERO; |
---|
| 361 | } |
---|
| 362 | |
---|
| 363 | /*~~~> Autonomous or time dependent ODE. Default is time dependent. */ |
---|
| 364 | Autonomous = !(ICNTRL[0] == 0); |
---|
| 365 | |
---|
| 366 | /*~~~> For Scalar tolerances(ICNTRL[1] != 0) the code uses AbsTol[0] and |
---|
| 367 | RelTol[0] |
---|
| 368 | For Vector tolerances(ICNTRL[1] == 0) the code uses AbsTol[1:NVAR] and |
---|
| 369 | RelTol[1:NVAR] */ |
---|
| 370 | |
---|
| 371 | if (ICNTRL[1] == 0) { |
---|
| 372 | VectorTol = TRUE; |
---|
| 373 | UplimTol = NVAR; |
---|
| 374 | } |
---|
| 375 | else { |
---|
| 376 | VectorTol = FALSE; |
---|
| 377 | UplimTol = 1; |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | /*~~~> Initialize the particular Rosenbrock method selected */ |
---|
| 381 | switch( ICNTRL[2] ) { |
---|
| 382 | case 0: |
---|
| 383 | case 4: |
---|
| 384 | Rodas3(); |
---|
| 385 | break; |
---|
| 386 | case 1: |
---|
| 387 | Ros2(); |
---|
| 388 | break; |
---|
| 389 | case 2: |
---|
| 390 | Ros3(); |
---|
| 391 | break; |
---|
| 392 | case 3: |
---|
| 393 | Ros4(); |
---|
| 394 | break; |
---|
| 395 | case 5: |
---|
| 396 | Rodas4(); |
---|
| 397 | break; |
---|
| 398 | default: |
---|
| 399 | printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] ); |
---|
| 400 | return ros_ErrorMsg(-2, Tstart, ZERO); |
---|
| 401 | } /* End switch */ |
---|
| 402 | |
---|
| 403 | /*~~~> The maximum number of steps admitted */ |
---|
| 404 | if (ICNTRL[3] == 0) |
---|
| 405 | Max_no_steps = bufsize - 1; |
---|
| 406 | else if (Max_no_steps > 0) |
---|
| 407 | Max_no_steps = ICNTRL[3]; |
---|
| 408 | else { |
---|
| 409 | printf("User-selected max no. of steps: ICNTRL[3]=%d",ICNTRL[3] ); |
---|
| 410 | return ros_ErrorMsg(-1,Tstart,ZERO); |
---|
| 411 | } |
---|
| 412 | |
---|
| 413 | /*~~~>The particular Rosenbrock method chosen for integrating the cts adjoint*/ |
---|
| 414 | if (ICNTRL[5] == 0) |
---|
| 415 | CadjMethod = 4; |
---|
| 416 | else if ( (ICNTRL[5] >= 1) && (ICNTRL[5] <= 5) ) |
---|
| 417 | CadjMethod = ICNTRL[5]; |
---|
| 418 | else { |
---|
| 419 | printf( "Unknown CADJ Rosenbrock method: ICNTRL[5]=%d", CadjMethod ); |
---|
| 420 | return ros_ErrorMsg(-2,Tstart,ZERO); |
---|
| 421 | } |
---|
| 422 | |
---|
| 423 | /*~~~> Discrete or continuous adjoint formulation */ |
---|
| 424 | if ( ICNTRL[6] == 0 ) |
---|
| 425 | AdjointType = Adj_discrete; |
---|
| 426 | else if ( (ICNTRL[6] >= 1) && (ICNTRL[6] <= 4) ) |
---|
| 427 | AdjointType = ICNTRL[6]; |
---|
| 428 | else { |
---|
| 429 | printf( "User-selected adjoint type: ICNTRL[6]=%d", AdjointType ); |
---|
| 430 | return ros_ErrorMsg(-9,Tstart,ZERO); |
---|
| 431 | } |
---|
| 432 | |
---|
| 433 | /*~~~> Save or not the forward LU factorization */ |
---|
| 434 | SaveLU = (ICNTRL[7] != 0); |
---|
| 435 | |
---|
| 436 | /*~~~> Unit roundoff (1+Roundoff>1) */ |
---|
| 437 | Roundoff = WLAMCH('E'); |
---|
| 438 | |
---|
| 439 | /*~~~> Lower bound on the step size: (positive value) */ |
---|
| 440 | if (RCNTRL[0] == ZERO) |
---|
| 441 | Hmin = ZERO; |
---|
| 442 | else if (RCNTRL[0] > ZERO) |
---|
| 443 | Hmin = RCNTRL[0]; |
---|
| 444 | else { |
---|
| 445 | printf( "User-selected Hmin: RCNTRL[0]=%f", RCNTRL[0] ); |
---|
| 446 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | /*~~~> Upper bound on the step size: (positive value) */ |
---|
| 450 | if (RCNTRL[1] == ZERO) |
---|
| 451 | Hmax = ABS(Tend-Tstart); |
---|
| 452 | else if (RCNTRL[1] > ZERO) |
---|
| 453 | Hmax = MIN(ABS(RCNTRL[1]),ABS(Tend-Tstart)); |
---|
| 454 | else { |
---|
| 455 | printf( "User-selected Hmax: RCNTRL[1]=%f", RCNTRL[1] ); |
---|
| 456 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 457 | } |
---|
| 458 | |
---|
| 459 | /*~~~> Starting step size: (positive value) */ |
---|
| 460 | if (RCNTRL[2] == ZERO) { |
---|
| 461 | Hstart = MAX(Hmin,DeltaMin); |
---|
| 462 | } |
---|
| 463 | else if (RCNTRL[2] > ZERO) |
---|
| 464 | Hstart = MIN(ABS(RCNTRL[2]),ABS(Tend-Tstart)); |
---|
| 465 | else { |
---|
| 466 | printf( "User-selected Hstart: RCNTRL[2]=%f", RCNTRL[2] ); |
---|
| 467 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 468 | } |
---|
| 469 | |
---|
| 470 | /*~~~> Step size can be changed s.t. FacMin < Hnew/Hold < FacMax */ |
---|
| 471 | if (RCNTRL[3] == ZERO) |
---|
| 472 | FacMin = (KPP_REAL)0.2; |
---|
| 473 | else if (RCNTRL[3] > ZERO) |
---|
| 474 | FacMin = RCNTRL[3]; |
---|
| 475 | else { |
---|
| 476 | printf( "User-selected FacMin: RCNTRL[3]=%f", RCNTRL[3] ); |
---|
| 477 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 478 | } |
---|
| 479 | if (RCNTRL[4] == ZERO) |
---|
| 480 | FacMax = (KPP_REAL)6.0; |
---|
| 481 | else if (RCNTRL[4] > ZERO) |
---|
| 482 | FacMax = RCNTRL[4]; |
---|
| 483 | else { |
---|
| 484 | printf( "User-selected FacMax: RCNTRL[4]=%f", RCNTRL[4] ); |
---|
| 485 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 486 | } |
---|
| 487 | |
---|
| 488 | /*~~~> FacRej: Factor to decrease step after 2 succesive rejections */ |
---|
| 489 | if (RCNTRL[5] == ZERO) |
---|
| 490 | FacRej = (KPP_REAL)0.1; |
---|
| 491 | else if (RCNTRL[5] > ZERO) |
---|
| 492 | FacRej = RCNTRL[5]; |
---|
| 493 | else { |
---|
| 494 | printf( "User-selected FacRej: RCNTRL[5]=%f", RCNTRL[5] ); |
---|
| 495 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 496 | } |
---|
| 497 | |
---|
| 498 | /*~~~> FacSafe: Safety Factor in the computation of new step size */ |
---|
| 499 | if (RCNTRL[6] == ZERO) |
---|
| 500 | FacSafe = (KPP_REAL)0.9; |
---|
| 501 | else if (RCNTRL[6] > ZERO) |
---|
| 502 | FacSafe = RCNTRL[6]; |
---|
| 503 | else { |
---|
| 504 | printf( "User-selected FacSafe: RCNTRL[6]=%f", RCNTRL[6] ); |
---|
| 505 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 506 | } |
---|
| 507 | |
---|
| 508 | /*~~~> Check if tolerances are reasonable */ |
---|
| 509 | for(i=0; i < UplimTol; i++) { |
---|
| 510 | if ( (AbsTol[i] <= ZERO) || (RelTol[i] <= (KPP_REAL)10.0*Roundoff) |
---|
| 511 | || (RelTol[i] >= (KPP_REAL)1.0) ) { |
---|
| 512 | printf( " AbsTol[%d] = %f", i, AbsTol[i] ); |
---|
| 513 | printf( " RelTol[%d] = %f", i, RelTol[i] ); |
---|
| 514 | return ros_ErrorMsg(-5,Tstart,ZERO); |
---|
| 515 | } |
---|
| 516 | } |
---|
| 517 | |
---|
| 518 | /*~~~> Allocate checkpoint space or open checkpoint files */ |
---|
| 519 | if (AdjointType == Adj_discrete) { |
---|
| 520 | ros_AllocateDBuffers( ros_S, SaveLU ); |
---|
| 521 | } |
---|
| 522 | else if ( (AdjointType == Adj_continuous) || |
---|
| 523 | (AdjointType == Adj_simple_continuous) ) { |
---|
| 524 | ros_AllocateCBuffers(); |
---|
| 525 | } |
---|
| 526 | |
---|
| 527 | /*~~~> CALL Forward Rosenbrock method */ |
---|
| 528 | IERR = ros_FwdInt(Y, Tstart, Tend, Texit, AbsTol, RelTol, AdjointType, Hmin, |
---|
| 529 | Hstart, Hmax, Roundoff, ISTATUS, Max_no_steps, |
---|
| 530 | RSTATUS, Autonomous, VectorTol, FacMax, FacMin, |
---|
| 531 | FacSafe, FacRej, SaveLU); |
---|
| 532 | |
---|
| 533 | printf( "\n\nFORWARD STATISTICS\n" ); |
---|
| 534 | printf( "Step=%d Acc=%d Rej=%d Singular=%d\n\n", Nstp, Nacc, Nrej, Nsng ); |
---|
| 535 | |
---|
| 536 | /*~~~> If Forward integration failed return */ |
---|
| 537 | if (IERR<0) |
---|
| 538 | return IERR; |
---|
| 539 | |
---|
| 540 | /*~~~> Initialize the particular Rosenbrock method for continuous adjoint */ |
---|
| 541 | if ( (AdjointType == Adj_continuous) || |
---|
| 542 | (AdjointType == Adj_simple_continuous) ) { |
---|
| 543 | switch (CadjMethod) { |
---|
| 544 | case 1: |
---|
| 545 | Ros2(); |
---|
| 546 | break; |
---|
| 547 | case 2: |
---|
| 548 | Ros3(); |
---|
| 549 | break; |
---|
| 550 | case 3: |
---|
| 551 | Ros4(); |
---|
| 552 | break; |
---|
| 553 | case 4: |
---|
| 554 | Rodas3(); |
---|
| 555 | break; |
---|
| 556 | case 5: |
---|
| 557 | Rodas4(); |
---|
| 558 | break; |
---|
| 559 | default: |
---|
| 560 | printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] ); |
---|
| 561 | return ros_ErrorMsg(-2,Tstart,ZERO); |
---|
| 562 | } |
---|
| 563 | } /* End switch */ |
---|
| 564 | |
---|
| 565 | switch( AdjointType ) { |
---|
| 566 | case Adj_discrete: |
---|
| 567 | IERR = ros_DadjInt (NADJ, Lambda, Tstart, Tend, Texit, SaveLU, ISTATUS, |
---|
| 568 | Roundoff, Autonomous ); |
---|
| 569 | break; |
---|
| 570 | case Adj_continuous: |
---|
| 571 | IERR = ros_CadjInt (NADJ, Lambda, Tend, Tstart, Texit, AbsTol_adj, |
---|
| 572 | RelTol_adj, RSTATUS, Hmin, Hmax, Hstart, Roundoff, |
---|
| 573 | Max_no_steps, Autonomous, VectorTol, FacMax, FacMin, |
---|
| 574 | FacSafe, FacRej, ISTATUS); |
---|
| 575 | break; |
---|
| 576 | case Adj_simple_continuous: |
---|
| 577 | IERR = ros_SimpleCadjInt (NADJ, Lambda, Tstart, Tend, Texit, ISTATUS, |
---|
| 578 | Autonomous, Roundoff); |
---|
| 579 | } /* End switch for AdjointType */ |
---|
| 580 | |
---|
| 581 | printf( "ADJOINT STATISTICS\n" ); |
---|
| 582 | printf( "Step=%d Acc=%d Rej=%d Singular=%d\n",Nstp,Nacc,Nrej,Nsng ); |
---|
| 583 | |
---|
| 584 | /*~~~> Free checkpoint space or close checkpoint files */ |
---|
| 585 | if (AdjointType == Adj_discrete) |
---|
| 586 | ros_FreeDBuffers( SaveLU ); |
---|
| 587 | else if ( (AdjointType == Adj_continuous) || |
---|
| 588 | (AdjointType == Adj_simple_continuous) ) |
---|
| 589 | ros_FreeCBuffers(); |
---|
| 590 | |
---|
| 591 | return IERR; |
---|
| 592 | } |
---|
| 593 | |
---|
| 594 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 595 | void ros_AllocateDBuffers( int S, int SaveLU ) { |
---|
| 596 | /*~~~> Allocate buffer space for discrete adjoint |
---|
| 597 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 598 | |
---|
| 599 | int i; |
---|
| 600 | |
---|
| 601 | chk_H = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL)); |
---|
| 602 | if (chk_H == NULL) { |
---|
| 603 | printf("Failed allocation of buffer H"); |
---|
| 604 | exit(0); |
---|
| 605 | } |
---|
| 606 | |
---|
| 607 | chk_T = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL)); |
---|
| 608 | if (chk_T == NULL) { |
---|
| 609 | printf("Failed allocation of buffer T"); |
---|
| 610 | exit(0); |
---|
| 611 | } |
---|
| 612 | |
---|
| 613 | chk_Y = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*)); |
---|
| 614 | if (chk_Y == NULL) { |
---|
| 615 | printf("Failed allocation of buffer Y"); |
---|
| 616 | exit(0); |
---|
| 617 | } |
---|
| 618 | for(i=0; i<bufsize; i++) { |
---|
| 619 | chk_Y[i] = (KPP_REAL*) malloc(NVAR * S * sizeof(KPP_REAL)); |
---|
| 620 | if (chk_Y[i] == NULL) { |
---|
| 621 | printf("Failed allocation of buffer Y"); |
---|
| 622 | exit(0); |
---|
| 623 | } |
---|
| 624 | } |
---|
| 625 | |
---|
| 626 | chk_K = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*)); |
---|
| 627 | if (chk_K == NULL) { |
---|
| 628 | printf("Failed allocation of buffer K"); |
---|
| 629 | exit(0); |
---|
| 630 | } |
---|
| 631 | for(i=0; i<bufsize; i++) { |
---|
| 632 | chk_K[i] = (KPP_REAL*) malloc(NVAR * S * sizeof(KPP_REAL)); |
---|
| 633 | if (chk_K == NULL) { |
---|
| 634 | printf("Failed allocation of buffer K"); |
---|
| 635 | exit(0); |
---|
| 636 | } |
---|
| 637 | } |
---|
| 638 | |
---|
| 639 | if (SaveLU) { |
---|
| 640 | chk_J = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*)); |
---|
| 641 | if (chk_J == NULL) { |
---|
| 642 | printf( "Failed allocation of buffer J"); |
---|
| 643 | exit(0); |
---|
| 644 | } |
---|
| 645 | #ifdef FULL_ALGEBRA |
---|
| 646 | for(i=0; i<bufsize; i++) { |
---|
| 647 | chk_J[i] = (KPP_REAL*) malloc(NVAR * NVAR * sizeof(KPP_REAL)); |
---|
| 648 | if (chk_J == NULL) { |
---|
| 649 | printf( "Failed allocation of buffer J"); |
---|
| 650 | exit(0); |
---|
| 651 | } |
---|
| 652 | } |
---|
| 653 | #else |
---|
| 654 | for(i=0; i<bufsize; i++) { |
---|
| 655 | chk_J[i] = (KPP_REAL*) malloc(LU_NONZERO * sizeof(KPP_REAL)); |
---|
| 656 | if (chk_J == NULL) { |
---|
| 657 | printf( "Failed allocation of buffer J"); |
---|
| 658 | exit(0); |
---|
| 659 | } |
---|
| 660 | } |
---|
| 661 | #endif |
---|
| 662 | } |
---|
| 663 | |
---|
| 664 | } /* End of ros_AllocateDBuffers */ |
---|
| 665 | |
---|
| 666 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 667 | void ros_FreeDBuffers( int SaveLU ) { |
---|
| 668 | /*~~~> Deallocate buffer space for discrete adjoint |
---|
| 669 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 670 | |
---|
| 671 | int i; |
---|
| 672 | |
---|
| 673 | free(chk_H); |
---|
| 674 | free(chk_T); |
---|
| 675 | |
---|
| 676 | for(i=0; i<bufsize; i++) { |
---|
| 677 | free(chk_Y[i]); |
---|
| 678 | free(chk_K[i]); |
---|
| 679 | } |
---|
| 680 | free(chk_Y); |
---|
| 681 | free(chk_K); |
---|
| 682 | |
---|
| 683 | if (SaveLU) { |
---|
| 684 | for(i=0; i<bufsize; i++) |
---|
| 685 | free(chk_J[i]); |
---|
| 686 | free(chk_J); |
---|
| 687 | } |
---|
| 688 | |
---|
| 689 | } /* End of ros_FreeDBuffers */ |
---|
| 690 | |
---|
| 691 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 692 | void ros_AllocateCBuffers() { |
---|
| 693 | /*~~~> Allocate buffer space for continuous adjoint |
---|
| 694 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 695 | |
---|
| 696 | int i; |
---|
| 697 | |
---|
| 698 | chk_H = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL)); |
---|
| 699 | if (chk_H == NULL) { |
---|
| 700 | printf( "Failed allocation of buffer H"); |
---|
| 701 | exit(0); |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | chk_T = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL)); |
---|
| 705 | if (chk_T == NULL) { |
---|
| 706 | printf( "Failed allocation of buffer T"); |
---|
| 707 | exit(0); |
---|
| 708 | } |
---|
| 709 | |
---|
| 710 | chk_Y = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize); |
---|
| 711 | if (chk_Y == NULL) { |
---|
| 712 | printf( "Failed allocation of buffer Y"); |
---|
| 713 | exit(0); |
---|
| 714 | } |
---|
| 715 | for(i=0; i<bufsize; i++) { |
---|
| 716 | chk_Y[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL)* NVAR ); |
---|
| 717 | if (chk_Y == NULL) { |
---|
| 718 | printf( "Failed allocation of buffer Y"); |
---|
| 719 | exit(0); |
---|
| 720 | } |
---|
| 721 | } |
---|
| 722 | |
---|
| 723 | chk_dY = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize); |
---|
| 724 | if (chk_dY == NULL) { |
---|
| 725 | printf( "Failed allocation of buffer dY"); |
---|
| 726 | exit(0); |
---|
| 727 | } |
---|
| 728 | for(i=0; i<bufsize; i++) { |
---|
| 729 | chk_dY[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL) * NVAR); |
---|
| 730 | if (chk_dY == NULL) { |
---|
| 731 | printf( "Failed allocation of buffer dY"); |
---|
| 732 | exit(0); |
---|
| 733 | } |
---|
| 734 | } |
---|
| 735 | |
---|
| 736 | chk_d2Y = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize); |
---|
| 737 | if (chk_d2Y == NULL) { |
---|
| 738 | printf( "Failed allocation of buffer d2Y"); |
---|
| 739 | exit(0); |
---|
| 740 | } |
---|
| 741 | for(i=0; i<bufsize; i++) { |
---|
| 742 | chk_d2Y[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL) * NVAR); |
---|
| 743 | if (chk_d2Y == NULL) { |
---|
| 744 | printf( "Failed allocation of buffer d2Y"); |
---|
| 745 | exit(0); |
---|
| 746 | } |
---|
| 747 | } |
---|
| 748 | } /* End of ros_AllocateCBuffers */ |
---|
| 749 | |
---|
| 750 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 751 | void ros_FreeCBuffers() { |
---|
| 752 | /*~~~> Dallocate buffer space for continuous adjoint |
---|
| 753 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 754 | |
---|
| 755 | int i; |
---|
| 756 | |
---|
| 757 | free(chk_H); |
---|
| 758 | free(chk_T); |
---|
| 759 | |
---|
| 760 | for(i=0; i<bufsize; i++) { |
---|
| 761 | free(chk_Y[i]); |
---|
| 762 | free(chk_dY[i]); |
---|
| 763 | free(chk_d2Y[i]); |
---|
| 764 | } |
---|
| 765 | |
---|
| 766 | free(chk_Y); |
---|
| 767 | free(chk_dY); |
---|
| 768 | free(chk_d2Y); |
---|
| 769 | |
---|
| 770 | } /* End of ros_FreeCBuffers */ |
---|
| 771 | |
---|
| 772 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 773 | void ros_DPush( int S, KPP_REAL T, KPP_REAL H, KPP_REAL Ystage[], |
---|
| 774 | KPP_REAL K[], KPP_REAL E[], int P[], int SaveLU ) { |
---|
| 775 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 776 | ~~~> Saves the next trajectory snapshot for discrete adjoints |
---|
| 777 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 778 | |
---|
| 779 | int i; |
---|
| 780 | |
---|
| 781 | stack_ptr = stack_ptr + 1; |
---|
| 782 | if ( stack_ptr >= bufsize ) { |
---|
| 783 | printf( "Push failed: buffer overflow" ); |
---|
| 784 | exit(0); |
---|
| 785 | } |
---|
| 786 | |
---|
| 787 | chk_H[ stack_ptr ] = H; |
---|
| 788 | chk_T[ stack_ptr ] = T; |
---|
| 789 | for(i=0; i<NVAR*S; i++) { |
---|
| 790 | chk_Y[stack_ptr][i] = Ystage[i]; |
---|
| 791 | chk_K[stack_ptr][i] = K[i]; |
---|
| 792 | } |
---|
| 793 | |
---|
| 794 | if (SaveLU) { |
---|
| 795 | #ifdef FULL_ALGEBRA |
---|
| 796 | int j; |
---|
| 797 | for(j=0; j<NVAR; j++) { |
---|
| 798 | for(i=0; i<NVAR; i++) |
---|
| 799 | chk_J[stack_ptr][i][j] = E[i][j]; |
---|
| 800 | chk_P[stack_ptr][j] = P[j]; |
---|
| 801 | } |
---|
| 802 | #else |
---|
| 803 | for(i=0; i<LU_NONZERO; i++) |
---|
| 804 | chk_J[stack_ptr][i] = E[i]; |
---|
| 805 | #endif |
---|
| 806 | } |
---|
| 807 | |
---|
| 808 | } /* End of ros_DPush */ |
---|
| 809 | |
---|
| 810 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 811 | void ros_DPop( int S, KPP_REAL* T, KPP_REAL* H, KPP_REAL* Ystage, |
---|
| 812 | KPP_REAL* K, KPP_REAL* E, int* P, int SaveLU ) { |
---|
| 813 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 814 | ~~~> Retrieves the next trajectory snapshot for discrete adjoints |
---|
| 815 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 816 | |
---|
| 817 | int i; |
---|
| 818 | |
---|
| 819 | if ( stack_ptr < 0 ) { |
---|
| 820 | printf( "Pop failed: empty buffer" ); |
---|
| 821 | exit(0); |
---|
| 822 | } |
---|
| 823 | |
---|
| 824 | *H = chk_H[ stack_ptr ]; |
---|
| 825 | *T = chk_T[ stack_ptr ]; |
---|
| 826 | |
---|
| 827 | for(i=0; i<NVAR*S; i++) { |
---|
| 828 | Ystage[i] = chk_Y[stack_ptr][i]; |
---|
| 829 | K[i] = chk_K[stack_ptr][i]; |
---|
| 830 | } |
---|
| 831 | |
---|
| 832 | if (SaveLU) { |
---|
| 833 | #ifdef FULL_ALGEBRA |
---|
| 834 | int j; |
---|
| 835 | for(i=0; i<NVAR; i++) { |
---|
| 836 | for(j=0; j<NVAR; j++) |
---|
| 837 | E[(j*NVAR)+i] = chk_J[stack_ptr][j][i]; |
---|
| 838 | P[i] = chk_P[stack_ptr][i]; |
---|
| 839 | } |
---|
| 840 | #else |
---|
| 841 | for(i=0; i<LU_NONZERO; i++) |
---|
| 842 | E[i] = chk_J[stack_ptr][i]; |
---|
| 843 | #endif |
---|
| 844 | } |
---|
| 845 | |
---|
| 846 | stack_ptr--; |
---|
| 847 | |
---|
| 848 | } /* End of ros_DPop */ |
---|
| 849 | |
---|
| 850 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 851 | void ros_CPush( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], |
---|
| 852 | KPP_REAL d2Y[] ) { |
---|
| 853 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 854 | ~~~> Saves the next trajectory snapshot for discrete adjoints |
---|
| 855 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 856 | |
---|
| 857 | int i; |
---|
| 858 | |
---|
| 859 | stack_ptr++; |
---|
| 860 | if ( stack_ptr > bufsize ) { |
---|
| 861 | printf( "Push failed: buffer overflow" ); |
---|
| 862 | exit(0); |
---|
| 863 | } |
---|
| 864 | chk_H[ stack_ptr ] = H; |
---|
| 865 | chk_T[ stack_ptr ] = T; |
---|
| 866 | |
---|
| 867 | for(i = 0; i< NVAR; i++ ) { |
---|
| 868 | chk_Y[stack_ptr][i] = Y[i]; |
---|
| 869 | chk_dY[stack_ptr][i] = dY[i]; |
---|
| 870 | chk_d2Y[stack_ptr][i] = d2Y[i]; |
---|
| 871 | } |
---|
| 872 | } /* End of ros_CPush */ |
---|
| 873 | |
---|
| 874 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 875 | void ros_CPop( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], |
---|
| 876 | KPP_REAL d2Y[] ) { |
---|
| 877 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 878 | ~~~> Retrieves the next trajectory snapshot for discrete adjoints |
---|
| 879 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 880 | |
---|
| 881 | int i; |
---|
| 882 | |
---|
| 883 | if ( stack_ptr <= 0 ) { |
---|
| 884 | printf( "Pop failed: empty buffer" ); |
---|
| 885 | exit(0); |
---|
| 886 | } |
---|
| 887 | H = chk_H[ stack_ptr ]; |
---|
| 888 | T = chk_T[ stack_ptr ]; |
---|
| 889 | |
---|
| 890 | for(i=0; i<NVAR; i++) { |
---|
| 891 | Y[i] = chk_Y[stack_ptr][i]; |
---|
| 892 | dY[i] = chk_dY[stack_ptr][i]; |
---|
| 893 | d2Y[i] = chk_d2Y[stack_ptr][i]; |
---|
| 894 | } |
---|
| 895 | stack_ptr--; |
---|
| 896 | } /* End of ros_CPop */ |
---|
| 897 | |
---|
| 898 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 899 | int ros_ErrorMsg( int Code, KPP_REAL T, KPP_REAL H) { |
---|
| 900 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 901 | Handles all error messages |
---|
| 902 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 903 | |
---|
| 904 | int IERR = Code; |
---|
| 905 | printf( "Forced exit from RosenbrockADJ due to the following error:"); |
---|
| 906 | |
---|
| 907 | switch (Code) { |
---|
| 908 | case -1: |
---|
| 909 | printf( "--> Improper value for maximal no of steps" ); |
---|
| 910 | break; |
---|
| 911 | case -2: |
---|
| 912 | printf( "--> Selected RosenbrockADJ method not implemented" ); |
---|
| 913 | break; |
---|
| 914 | case -3: |
---|
| 915 | printf( "--> Hmin/Hmax/Hstart must be positive" ); |
---|
| 916 | break; |
---|
| 917 | case -4: |
---|
| 918 | printf( "--> FacMin/FacMax/FacRej must be positive" ); |
---|
| 919 | break; |
---|
| 920 | case -5: |
---|
| 921 | printf( "--> Improper tolerance values" ); |
---|
| 922 | break; |
---|
| 923 | case -6: |
---|
| 924 | printf( "--> No of steps exceeds maximum buffer bound" ); |
---|
| 925 | break; |
---|
| 926 | case -7: |
---|
| 927 | printf( "--> Step size too small: T + 10*H = T or H < Roundoff" ); |
---|
| 928 | break; |
---|
| 929 | case -8: |
---|
| 930 | printf( "--> Matrix is repeatedly singular" ); |
---|
| 931 | break; |
---|
| 932 | case -9: |
---|
| 933 | printf( "--> Improper type of adjoint selected" ); |
---|
| 934 | break; |
---|
| 935 | default: |
---|
| 936 | printf( "Unknown Error code: %d", Code ); |
---|
| 937 | } /* End of switch */ |
---|
| 938 | |
---|
| 939 | return IERR; |
---|
| 940 | } /* End of ros_ErrorMsg */ |
---|
| 941 | |
---|
| 942 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 943 | int ros_FwdInt ( KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL T, |
---|
| 944 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int AdjointType, |
---|
| 945 | KPP_REAL Hmin, KPP_REAL Hstart, KPP_REAL Hmax, |
---|
| 946 | KPP_REAL Roundoff, int ISTATUS[], int Max_no_steps, |
---|
| 947 | KPP_REAL RSTATUS[], int Autonomous, int VectorTol, |
---|
| 948 | KPP_REAL FacMax, KPP_REAL FacMin, KPP_REAL FacSafe, |
---|
| 949 | KPP_REAL FacRej, int SaveLU ) { |
---|
| 950 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 951 | Template for the implementation of a generic RosenbrockADJ method |
---|
| 952 | defined by ros_S (no of stages) |
---|
| 953 | and its coefficients ros_{A,C,M,E,Alpha,Gamma} |
---|
| 954 | ~~~> Y - Input: the initial condition at Tstart; Output: the solution at T |
---|
| 955 | ~~~> Tstart, Tend - Input: integration interval |
---|
| 956 | ~~~> T - Output: time at which the solution is returned (T=Tend if success) |
---|
| 957 | ~~~> AbsTol, RelTol - Input: tolerances |
---|
| 958 | ~~~> IERR - Output: Error indicator |
---|
| 959 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 960 | |
---|
| 961 | /* ~~~~ Local variables */ |
---|
| 962 | KPP_REAL Ynew[NVAR], Fcn0[NVAR], Fcn[NVAR]; |
---|
| 963 | KPP_REAL K[NVAR*ros_S], dFdT[NVAR]; |
---|
| 964 | KPP_REAL *Ystage = NULL; /* Array pointer */ |
---|
| 965 | #ifdef FULL_ALGEBRA |
---|
| 966 | KPP_REAL Jac0[NVAR][NVAR], Ghimj[NVAR][NVAR]; |
---|
| 967 | #else |
---|
| 968 | KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO]; |
---|
| 969 | #endif |
---|
| 970 | KPP_REAL H, Hnew, HC, HG, Fac, Tau; |
---|
| 971 | KPP_REAL Err, Yerr[NVAR]; |
---|
| 972 | int Pivot[NVAR], Direction, ioffset, i, j=0, istage; |
---|
| 973 | int RejectLastH, RejectMoreH, Singular; /* Boolean Values */ |
---|
| 974 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 975 | |
---|
| 976 | /*~~~> Allocate stage vector buffer if needed */ |
---|
| 977 | if (AdjointType == Adj_discrete) { |
---|
| 978 | Ystage = (KPP_REAL*) malloc(NVAR*ros_S*sizeof(KPP_REAL)); |
---|
| 979 | /* Uninitialized Ystage may lead to NaN on some compilers */ |
---|
| 980 | if (Ystage == NULL) { |
---|
| 981 | printf( "Allocation of Ystage failed" ); |
---|
| 982 | exit(0); |
---|
| 983 | } |
---|
| 984 | } |
---|
| 985 | |
---|
| 986 | /*~~~> Initial preparations */ |
---|
| 987 | T = Tstart; |
---|
| 988 | RSTATUS[Nhexit] = ZERO; |
---|
| 989 | H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) ); |
---|
| 990 | if (ABS(H) <= ((KPP_REAL)10.0)*Roundoff) |
---|
| 991 | H = DeltaMin; |
---|
| 992 | |
---|
| 993 | if (Tend >= Tstart) |
---|
| 994 | Direction = 1; |
---|
| 995 | else |
---|
| 996 | Direction = -1; |
---|
| 997 | |
---|
| 998 | H = Direction*H; |
---|
| 999 | |
---|
| 1000 | RejectLastH = FALSE; |
---|
| 1001 | RejectMoreH = FALSE; |
---|
| 1002 | |
---|
| 1003 | /*~~~> Time loop begins below */ |
---|
| 1004 | while ( ((Direction > 0) && ((T-Tend)+Roundoff <= ZERO)) || |
---|
| 1005 | ((Direction < 0) && ((Tend-T)+Roundoff <= ZERO)) ) { /* TimeLoop */ |
---|
| 1006 | |
---|
| 1007 | if ( ISTATUS[Nstp] > Max_no_steps ) /* Too many steps */ |
---|
| 1008 | return ros_ErrorMsg(-6,T,H); |
---|
| 1009 | |
---|
| 1010 | if ( ((T+((KPP_REAL)0.1)*H) == T) || (H <= Roundoff) ) /* Step size |
---|
| 1011 | too small */ |
---|
| 1012 | return ros_ErrorMsg(-7,T,H); |
---|
| 1013 | |
---|
| 1014 | /*~~~> Limit H if necessary to avoid going beyond Tend */ |
---|
| 1015 | RSTATUS[Nhexit] = H; |
---|
| 1016 | H = MIN(H,ABS(Tend-T)); |
---|
| 1017 | |
---|
| 1018 | /*~~~> Compute the function at current time */ |
---|
| 1019 | FunTemplate(T,Y,Fcn0); |
---|
| 1020 | ISTATUS[Nfun] = ISTATUS[Nfun] + 1; |
---|
| 1021 | |
---|
| 1022 | /*~~~> Compute the function derivative with respect to T */ |
---|
| 1023 | if (!Autonomous) |
---|
| 1024 | ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT, ISTATUS ); |
---|
| 1025 | |
---|
| 1026 | /*~~~> Compute the Jacobian at current time */ |
---|
| 1027 | JacTemplate(T,Y,Jac0); |
---|
| 1028 | ISTATUS[Njac] = ISTATUS[Njac] + 1; |
---|
| 1029 | |
---|
| 1030 | /*~~~> Repeat step calculation until current step accepted */ |
---|
| 1031 | do { /* UntilAccepted */ |
---|
| 1032 | |
---|
| 1033 | Singular = ros_PrepareMatrix ( H,Direction,ros_Gamma[0],Jac0,Ghimj,Pivot, |
---|
| 1034 | ISTATUS ); |
---|
| 1035 | |
---|
| 1036 | if (Singular) /* More than 5 consecutive failed decompositions */ |
---|
| 1037 | return ros_ErrorMsg(-8,T,H); |
---|
| 1038 | |
---|
| 1039 | /*~~~> Compute the stages */ |
---|
| 1040 | for( istage = 0; istage < ros_S; istage++ ) { /* Stage */ |
---|
| 1041 | |
---|
| 1042 | /* Current istage offset. Current istage vector is |
---|
| 1043 | K(ioffset+1:ioffset+NVAR) */ |
---|
| 1044 | ioffset = NVAR*istage; |
---|
| 1045 | |
---|
| 1046 | /*For the 1st istage the function has been computed previously*/ |
---|
| 1047 | if ( istage == 0 ) { |
---|
| 1048 | WCOPY(NVAR,Fcn0,1,Fcn,1); |
---|
| 1049 | if (AdjointType == Adj_discrete) { /* Save stage solution */ |
---|
| 1050 | for(i=0; i<NVAR; i++) |
---|
| 1051 | Ystage[i] = Y[i]; |
---|
| 1052 | WCOPY(NVAR,Y,1,Ynew,1); |
---|
| 1053 | } |
---|
| 1054 | } |
---|
| 1055 | /* istage>0 and a new function evaluation is needed at the |
---|
| 1056 | current istage */ |
---|
| 1057 | else if ( ros_NewF[istage] ) { |
---|
| 1058 | WCOPY(NVAR,Y,1,Ynew,1); |
---|
| 1059 | for ( j = 0; j < istage; j++ ) { |
---|
| 1060 | WAXPY( NVAR,ros_A[(istage)*(istage-1)/2+j], |
---|
| 1061 | &K[NVAR*j],1,Ynew,1 ); |
---|
| 1062 | } |
---|
| 1063 | Tau = T + ros_Alpha[istage]*Direction*H; |
---|
| 1064 | FunTemplate(Tau,Ynew,Fcn); |
---|
| 1065 | ISTATUS[Nfun] = ISTATUS[Nfun] + 1; |
---|
| 1066 | } /* if istage == 1 elseif ros_NewF[istage] */ |
---|
| 1067 | |
---|
| 1068 | /* Save stage solution every time even if ynew is not updated */ |
---|
| 1069 | if ( ( istage > 0 ) && (AdjointType == Adj_discrete) ) { |
---|
| 1070 | for(i=0; i<NVAR; i++) |
---|
| 1071 | Ystage[ioffset+i] = Ynew[i]; |
---|
| 1072 | } |
---|
| 1073 | WCOPY(NVAR,Fcn,1,&K[ioffset],1); |
---|
| 1074 | for( j = 0; j < istage; j++ ) { |
---|
| 1075 | HC = ros_C[(istage)*(istage-1)/2+j]/(Direction*H); |
---|
| 1076 | WAXPY(NVAR,HC,&K[NVAR*j],1,&K[ioffset],1); |
---|
| 1077 | } |
---|
| 1078 | if (( !Autonomous) && (ros_Gamma[istage] != ZERO)) { |
---|
| 1079 | HG = Direction*H*ros_Gamma[istage]; |
---|
| 1080 | WAXPY(NVAR,HG,dFdT,1,&K[ioffset],1); |
---|
| 1081 | } |
---|
| 1082 | ros_Solve('N', Ghimj, Pivot, &K[ioffset], ISTATUS); |
---|
| 1083 | } /* End of Stage loop */ |
---|
| 1084 | |
---|
| 1085 | /*~~~> Compute the new solution */ |
---|
| 1086 | WCOPY(NVAR,Y,1,Ynew,1); |
---|
| 1087 | for( j=0; j<ros_S; j++ ) |
---|
| 1088 | WAXPY(NVAR,ros_M[j],&K[NVAR*j],1,Ynew,1); |
---|
| 1089 | |
---|
| 1090 | /*~~~> Compute the error estimation */ |
---|
| 1091 | WSCAL(NVAR,ZERO,Yerr,1); |
---|
| 1092 | for( j=0; j<ros_S; j++ ) |
---|
| 1093 | WAXPY(NVAR,ros_E[j],&K[NVAR*j],1,Yerr,1); |
---|
| 1094 | Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ); |
---|
| 1095 | |
---|
| 1096 | /*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */ |
---|
| 1097 | Fac = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,(ONE/ros_ELO)))); |
---|
| 1098 | Hnew = H*Fac; |
---|
| 1099 | |
---|
| 1100 | /*~~~> Check the error magnitude and adjust step size */ |
---|
| 1101 | ISTATUS[Nstp] = ISTATUS[Nstp] + 1; |
---|
| 1102 | if ( (Err <= ONE) || (H <= Hmin) ) { /*~~~> Accept step */ |
---|
| 1103 | ISTATUS[Nacc]++; |
---|
| 1104 | if (AdjointType == Adj_discrete) { /* Save current state */ |
---|
| 1105 | ros_DPush( ros_S, T, H, Ystage, K, Ghimj, Pivot, SaveLU ); |
---|
| 1106 | } |
---|
| 1107 | else if ( (AdjointType == Adj_continuous) || |
---|
| 1108 | (AdjointType == Adj_simple_continuous) ) { |
---|
| 1109 | #ifdef FULL_ALGEBRA |
---|
| 1110 | K = MATMUL(Jac0,Fcn0); |
---|
| 1111 | #else |
---|
| 1112 | Jac_SP_Vec( Jac0, Fcn0, &K[0] ); |
---|
| 1113 | #endif |
---|
| 1114 | if ( !Autonomous) |
---|
| 1115 | WAXPY(NVAR,ONE,dFdT,1,&K[0],1); |
---|
| 1116 | ros_CPush( T, H, Y, Fcn0, &K[0] ); |
---|
| 1117 | } |
---|
| 1118 | WCOPY(NVAR,Ynew,1,Y,1); |
---|
| 1119 | T = T + Direction*H; |
---|
| 1120 | Hnew = MAX(Hmin,MIN(Hnew,Hmax)); |
---|
| 1121 | if (RejectLastH) { /* No step size increase after a |
---|
| 1122 | rejected step */ |
---|
| 1123 | Hnew = MIN(Hnew,H); |
---|
| 1124 | } |
---|
| 1125 | RSTATUS[Nhexit] = H; |
---|
| 1126 | RSTATUS[Nhnew] = Hnew; |
---|
| 1127 | RSTATUS[Ntexit] = T; |
---|
| 1128 | RejectLastH = FALSE; |
---|
| 1129 | RejectMoreH = FALSE; |
---|
| 1130 | H = Hnew; |
---|
| 1131 | break; /* UntilAccepted - EXIT THE LOOP: WHILE STEP NOT ACCEPTED */ |
---|
| 1132 | } |
---|
| 1133 | |
---|
| 1134 | else { /*~~~> Reject step */ |
---|
| 1135 | if (RejectMoreH) |
---|
| 1136 | Hnew = H*FacRej; |
---|
| 1137 | RejectMoreH = RejectLastH; |
---|
| 1138 | RejectLastH = TRUE; |
---|
| 1139 | H = Hnew; |
---|
| 1140 | if (ISTATUS[Nacc] >= 1) |
---|
| 1141 | ISTATUS[Nrej]++; |
---|
| 1142 | } /* End if else - Err <= 1 */ |
---|
| 1143 | |
---|
| 1144 | } while(1); /* End of UntilAccepted do loop */ |
---|
| 1145 | } /* End of TimeLoop */ |
---|
| 1146 | |
---|
| 1147 | /*~~~> Save last state: only needed for continuous adjoint */ |
---|
| 1148 | if ( (AdjointType == Adj_continuous) || |
---|
| 1149 | (AdjointType == Adj_simple_continuous) ) { |
---|
| 1150 | FunTemplate(T,Y,Fcn0); |
---|
| 1151 | ISTATUS[Nfun]++; |
---|
| 1152 | JacTemplate(T,Y,Jac0); |
---|
| 1153 | ISTATUS[Njac]++; |
---|
| 1154 | #ifdef FULL_ALGEBRA |
---|
| 1155 | K = MATMUL(Jac0,Fcn0); |
---|
| 1156 | #else |
---|
| 1157 | Jac_SP_Vec( Jac0, Fcn0, &K[0] ); |
---|
| 1158 | #endif |
---|
| 1159 | if (!Autonomous) { |
---|
| 1160 | ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT, ISTATUS ); |
---|
| 1161 | WAXPY(NVAR,ONE,dFdT,1,&K[0],1); |
---|
| 1162 | } |
---|
| 1163 | ros_CPush( T, H, Y, Fcn0, &K[0] ); |
---|
| 1164 | /*~~~> Deallocate stage buffer: only needed for discrete adjoint */ |
---|
| 1165 | } |
---|
| 1166 | else if (AdjointType == Adj_discrete) { |
---|
| 1167 | free(Ystage); |
---|
| 1168 | } |
---|
| 1169 | |
---|
| 1170 | /*~~~> Succesful exit */ |
---|
| 1171 | return 1; /*~~~> The integration was successful */ |
---|
| 1172 | } /* End of ros_FwdInt */ |
---|
| 1173 | |
---|
| 1174 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1175 | int ros_DadjInt ( int NADJ, KPP_REAL Lambda[][NVAR], KPP_REAL Tstart, |
---|
| 1176 | KPP_REAL Tend, KPP_REAL T, int SaveLU, int ISTATUS[], |
---|
| 1177 | KPP_REAL Roundoff, int Autonomous) { |
---|
| 1178 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1179 | Template for the implementation of a generic RosenbrockSOA method |
---|
| 1180 | defined by ros_S (no of stages) |
---|
| 1181 | and its coefficients ros_{A,C,M,E,Alpha,Gamma} |
---|
| 1182 | !~~~> NADJ - Input: the initial condition at Tstart; Output: the solution at T |
---|
| 1183 | !~~~> Lambda[NADJ][NVAR] - First order adjoint |
---|
| 1184 | !~~~> Tstart, Tend - Input: integration interval |
---|
| 1185 | !~~~> T - Output: time at which the solution is returned (T=Tend if success) |
---|
| 1186 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1187 | |
---|
| 1188 | /*~~~~ Local variables */ |
---|
| 1189 | KPP_REAL Ystage[NVAR*ros_S], K[NVAR*ros_S]; |
---|
| 1190 | KPP_REAL U[NADJ][NVAR*ros_S], V[NADJ][NVAR*ros_S]; |
---|
| 1191 | #ifdef FULL_ALGEBRA |
---|
| 1192 | KPP_REAL Jac[NVAR][NVAR], dJdT[NVAR][NVAR], Ghimj[NVAR][NVAR]; |
---|
| 1193 | #else |
---|
| 1194 | KPP_REAL Jac[LU_NONZERO], dJdT[LU_NONZERO], Ghimj[LU_NONZERO]; |
---|
| 1195 | #endif |
---|
| 1196 | KPP_REAL Hes0[NHESS]; |
---|
| 1197 | KPP_REAL Tmp[NVAR], Tmp2[NVAR]; |
---|
| 1198 | KPP_REAL H=0.0, HC, HA, Tau; |
---|
| 1199 | int Pivot[NVAR], Direction; |
---|
| 1200 | int i, j, m, istage, istart, jstart; |
---|
| 1201 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1202 | |
---|
| 1203 | if (Tend >= Tstart) |
---|
| 1204 | Direction = 1; |
---|
| 1205 | else |
---|
| 1206 | Direction = -1; |
---|
| 1207 | |
---|
| 1208 | /*~~~> Time loop begins below */ |
---|
| 1209 | while ( stack_ptr > 0 ) { /* TimeLoop */ |
---|
| 1210 | |
---|
| 1211 | /*~~~> Recover checkpoints for stage values and vectors */ |
---|
| 1212 | ros_DPop( ros_S, &T, &H, &Ystage[0], &K[0], &Ghimj[0], &Pivot[0], SaveLU ); |
---|
| 1213 | |
---|
| 1214 | /*~~~> Compute LU decomposition */ |
---|
| 1215 | if (!SaveLU) { |
---|
| 1216 | JacTemplate(T,&Ystage[0],Ghimj); |
---|
| 1217 | ISTATUS[Njac] = ISTATUS[Njac] + 1; |
---|
| 1218 | Tau = ONE/(Direction*H*ros_Gamma[0]); |
---|
| 1219 | #ifdef FULL_ALGEBRA |
---|
| 1220 | for(j=0; j<NVAR; j++) { |
---|
| 1221 | for(i=0; i<NVAR; i++) |
---|
| 1222 | Ghimj[i][j] = -Ghimj[i][j]; |
---|
| 1223 | } |
---|
| 1224 | for(i=0; i<NVAR; i++) |
---|
| 1225 | Ghimj[i][i] = Ghimj[i][i]+Tau; |
---|
| 1226 | #else |
---|
| 1227 | WSCAL(LU_NONZERO,(-ONE),Ghimj,1); |
---|
| 1228 | for (i=0; i<NVAR; i++) |
---|
| 1229 | Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+Tau; |
---|
| 1230 | #endif |
---|
| 1231 | ros_Decomp(Ghimj, Pivot, &j, ISTATUS); |
---|
| 1232 | } |
---|
| 1233 | |
---|
| 1234 | /*~~~> Compute Hessian at the beginning of the interval */ |
---|
| 1235 | HessTemplate(T,&Ystage[0],Hes0); |
---|
| 1236 | |
---|
| 1237 | /*~~~> Compute the stages */ |
---|
| 1238 | for (istage = ros_S - 1; istage >= 0; istage--) { /* Stage loop */ |
---|
| 1239 | |
---|
| 1240 | /*~~~> Current istage first entry */ |
---|
| 1241 | istart = NVAR*istage; |
---|
| 1242 | |
---|
| 1243 | /*~~~> Compute U */ |
---|
| 1244 | for (m = 0; m<NADJ; m++) { |
---|
| 1245 | WCOPY(NVAR,&Lambda[m][0],1,&U[m][istart],1); |
---|
| 1246 | WSCAL(NVAR,ros_M[istage],&U[m][istart],1); |
---|
| 1247 | } /* m=0:NADJ-1 */ |
---|
| 1248 | for (j = istage+1; j < ros_S; j++) { |
---|
| 1249 | jstart = NVAR*j; |
---|
| 1250 | HA = ros_A[j*(j-1)/2+istage]; |
---|
| 1251 | HC = ros_C[j*(j-1)/2+istage]/(Direction*H); |
---|
| 1252 | for ( m = 0; m < NADJ; m++ ) { |
---|
| 1253 | WAXPY(NVAR,HA,&V[m][jstart],1,&U[m][istart],1); |
---|
| 1254 | WAXPY(NVAR,HC,&U[m][jstart],1,&U[m][istart],1); |
---|
| 1255 | } /* m=0:NADJ-1 */ |
---|
| 1256 | } |
---|
| 1257 | for ( m = 0; m < NADJ; m++ ) |
---|
| 1258 | ros_Solve('T', Ghimj, Pivot, &U[m][istart], ISTATUS); /* m=1:NADJ-1 */ |
---|
| 1259 | |
---|
| 1260 | /*~~~> Compute V */ |
---|
| 1261 | Tau = T + ros_Alpha[istage]*Direction*H; |
---|
| 1262 | JacTemplate(Tau,&Ystage[istart],Jac); |
---|
| 1263 | ISTATUS[Njac]++; |
---|
| 1264 | for ( m = 0; m < NADJ; m++ ) { |
---|
| 1265 | #ifdef FULL_ALGEBRA |
---|
| 1266 | for (i=istart; i < istart+NVAR-1; i++ ) |
---|
| 1267 | V[[m][i] = MATMUL(TRANSPOSE(Jac),U[m][i]]; |
---|
| 1268 | #else |
---|
| 1269 | JacTR_SP_Vec(Jac,&U[m][istart],&V[m][istart]); |
---|
| 1270 | #endif |
---|
| 1271 | } /* m=0:NADJ-1 */ |
---|
| 1272 | } /*End of Stage loop */ |
---|
| 1273 | |
---|
| 1274 | if (!Autonomous) |
---|
| 1275 | /*~~~> Compute the Jacobian derivative with respect to T. |
---|
| 1276 | Last "Jac" computed for stage 1 */ |
---|
| 1277 | ros_JacTimeDerivative ( T, Roundoff, &Ystage[0], Jac, dJdT, ISTATUS ); |
---|
| 1278 | |
---|
| 1279 | /*~~~> Compute the new solution */ |
---|
| 1280 | /*~~~> Compute Lambda */ |
---|
| 1281 | for( istage = 0; istage < ros_S; istage++ ) { |
---|
| 1282 | istart = NVAR*istage; |
---|
| 1283 | for (m = 0; m < NADJ; m++) { |
---|
| 1284 | /* Add V_i */ |
---|
| 1285 | WAXPY(NVAR,ONE,&V[m][istart],1,&Lambda[m][0],1); |
---|
| 1286 | /* Add (H0xK_i)^T * U_i */ |
---|
| 1287 | HessTR_Vec ( Hes0, &U[m][istart], &K[istart], Tmp ); |
---|
| 1288 | WAXPY(NVAR,ONE,Tmp,1,&Lambda[m][0],1); |
---|
| 1289 | } /* m=0:NADJ-1 */ |
---|
| 1290 | } |
---|
| 1291 | |
---|
| 1292 | /* Add H * dJac_dT_0^T * \sum(gamma_i U_i) */ |
---|
| 1293 | /* Tmp holds sum gamma_i U_i */ |
---|
| 1294 | if (!Autonomous) { |
---|
| 1295 | for( m = 0; m < NADJ; m++ ) { |
---|
| 1296 | for(i=0; i<NVAR; i++) |
---|
| 1297 | Tmp[i] = ZERO; |
---|
| 1298 | for( istage = 0; istage < ros_S; istage++ ) { |
---|
| 1299 | istart = NVAR*istage; |
---|
| 1300 | WAXPY(NVAR,ros_Gamma[istage],&U[m][istart],1,Tmp,1); |
---|
| 1301 | } |
---|
| 1302 | #ifdef FULL_ALGEBRA |
---|
| 1303 | Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp); |
---|
| 1304 | #else |
---|
| 1305 | JacTR_SP_Vec(dJdT,Tmp,Tmp2); |
---|
| 1306 | #endif |
---|
| 1307 | WAXPY(NVAR,H,Tmp2,1,&Lambda[m][0],1); |
---|
| 1308 | } /* m=0:NADJ-1 */ |
---|
| 1309 | } /* .NOT.Autonomous */ |
---|
| 1310 | } /* End of TimeLoop */ |
---|
| 1311 | |
---|
| 1312 | /*~~~> Save last state */ |
---|
| 1313 | /*~~~> Succesful exit */ |
---|
| 1314 | return 1; /*~~~> The integration was successful */ |
---|
| 1315 | |
---|
| 1316 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1317 | } /* End of ros_DadjInt */ |
---|
| 1318 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1319 | |
---|
| 1320 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1321 | int ros_CadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, KPP_REAL Tend, |
---|
| 1322 | KPP_REAL T, KPP_REAL AbsTol_adj[][NVAR], |
---|
| 1323 | KPP_REAL RelTol_adj[][NVAR], KPP_REAL RSTATUS[], |
---|
| 1324 | KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, |
---|
| 1325 | KPP_REAL Roundoff, int Max_no_steps, int Autonomous, |
---|
| 1326 | int VectorTol, KPP_REAL FacMax, KPP_REAL FacMin, |
---|
| 1327 | KPP_REAL FacSafe, KPP_REAL FacRej, int ISTATUS[] ) { |
---|
| 1328 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1329 | Template for the implementation of a generic RosenbrockADJ method |
---|
| 1330 | defined by ros_S (no of stages) |
---|
| 1331 | and its coefficients ros_{A,C,M,E,Alpha,Gamma} |
---|
| 1332 | ~~~> NADJ, Y[NADJ][NVAR] - Input: the initial condition at Tstart; |
---|
| 1333 | Output: the solution at T |
---|
| 1334 | ~~~> Tstart, Tend - Input: integration interval |
---|
| 1335 | ~~~> AbsTol_adj[NADJ][NVAR], RelTol_adj[NADJ][NVAR] - Input: adjoint tolerances |
---|
| 1336 | ~~~> T - Output: time at which the solution is returned (T=Tend if success) |
---|
| 1337 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1338 | |
---|
| 1339 | /*~~~~ Local variables */ |
---|
| 1340 | KPP_REAL Y0[NVAR]; |
---|
| 1341 | KPP_REAL Ynew[NADJ][NVAR], Fcn0[NADJ][NVAR], Fcn[NADJ][NVAR]; |
---|
| 1342 | KPP_REAL K[NADJ][NVAR*ros_S], dFdT[NADJ][NVAR]; |
---|
| 1343 | #ifdef FULL_ALGEBRA |
---|
| 1344 | KPP_REAL Jac0[NVAR][NVAR], Ghimj[NVAR][NVAR], Jac[NVAR][NVAR], |
---|
| 1345 | dJdT[NVAR][NVAR]; |
---|
| 1346 | #else |
---|
| 1347 | KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO], Jac[LU_NONZERO], |
---|
| 1348 | dJdT[LU_NONZERO]; |
---|
| 1349 | #endif |
---|
| 1350 | KPP_REAL H, Hnew, HC, HG, Fac, Tau; |
---|
| 1351 | KPP_REAL Err, Yerr[NADJ][NVAR]; |
---|
| 1352 | int Pivot[NVAR], Direction, ioffset, j, istage, iadj; |
---|
| 1353 | int RejectLastH, RejectMoreH, Singular; /* Boolean values */ |
---|
| 1354 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1355 | |
---|
| 1356 | /*~~~> Initial preparations */ |
---|
| 1357 | T = Tstart; |
---|
| 1358 | RSTATUS[Nhexit] = (KPP_REAL)0.0; |
---|
| 1359 | H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) ); |
---|
| 1360 | if (ABS(H) <= ((KPP_REAL)10.0)*Roundoff) |
---|
| 1361 | H = DeltaMin; |
---|
| 1362 | |
---|
| 1363 | if (Tend >= Tstart) |
---|
| 1364 | Direction = 1; |
---|
| 1365 | else |
---|
| 1366 | Direction = -1; |
---|
| 1367 | |
---|
| 1368 | H = Direction*H; |
---|
| 1369 | RejectLastH = FALSE; |
---|
| 1370 | RejectMoreH = FALSE; |
---|
| 1371 | |
---|
| 1372 | /*~~~> Time loop begins below */ |
---|
| 1373 | while ( ((Direction > 0) && ((T-Tend)+Roundoff <= ZERO)) |
---|
| 1374 | || ((Direction < 0) && ((Tend-T)+Roundoff <= ZERO)) ) { /*TimeLoop*/ |
---|
| 1375 | |
---|
| 1376 | if ( ISTATUS[Nstp] > Max_no_steps ) /* Too many steps */ |
---|
| 1377 | return ros_ErrorMsg(-6,T,H); |
---|
| 1378 | |
---|
| 1379 | /* Step size too small */ |
---|
| 1380 | if ( ((T+((KPP_REAL)0.1)*H) == T) || (H <= Roundoff) ) |
---|
| 1381 | return ros_ErrorMsg(-7,T,H); |
---|
| 1382 | |
---|
| 1383 | /*~~~> Limit H if necessary to avoid going beyond Tend */ |
---|
| 1384 | RSTATUS[Nhexit] = H; |
---|
| 1385 | H = MIN(H,ABS(Tend-T)); |
---|
| 1386 | |
---|
| 1387 | /*~~~> Interpolate forward solution */ |
---|
| 1388 | ros_cadj_Y( T, Y0 ); |
---|
| 1389 | /*~~~> Compute the Jacobian at current time */ |
---|
| 1390 | JacTemplate(T, Y0, Jac0); |
---|
| 1391 | ISTATUS[Njac]++; |
---|
| 1392 | |
---|
| 1393 | /*~~~> Compute the function derivative with respect to T */ |
---|
| 1394 | if (!Autonomous) { |
---|
| 1395 | ros_JacTimeDerivative ( T, Roundoff, Y0, Jac0, dJdT, ISTATUS ); |
---|
| 1396 | for (iadj = 0; iadj < NADJ; iadj++) { |
---|
| 1397 | #ifdef FULL_ALGEBRA |
---|
| 1398 | for (i=0; i<NVAR; i++) |
---|
| 1399 | dFdT[iadj][i] = MATMUL(TRANSPOSE(dJdT),Y[iadj][i]); |
---|
| 1400 | #else |
---|
| 1401 | JacTR_SP_Vec(dJdT,&Y[iadj][0],&dFdT[iadj][0]); |
---|
| 1402 | #endif |
---|
| 1403 | WSCAL(NVAR,(-ONE),&dFdT[iadj][0],1); |
---|
| 1404 | } /* End for loop */ |
---|
| 1405 | } /* End if */ |
---|
| 1406 | |
---|
| 1407 | /*~~~> Ydot = -J^T*Y */ |
---|
| 1408 | #ifdef FULL_ALGEBRA |
---|
| 1409 | int i; |
---|
| 1410 | for(i=0; i<NVAR; i++) { |
---|
| 1411 | for(j=0; j<NVAR; j++) |
---|
| 1412 | Jac0[i][j] = -Jac0[i][j]; |
---|
| 1413 | } |
---|
| 1414 | #else |
---|
| 1415 | WSCAL(LU_NONZERO,(-ONE),Jac0,1); |
---|
| 1416 | #endif |
---|
| 1417 | |
---|
| 1418 | for ( iadj = 0; iadj < NADJ; iadj++ ) { |
---|
| 1419 | #ifdef FULL_ALGEBRA |
---|
| 1420 | int i; |
---|
| 1421 | for(i=0; i<NVAR; i++) |
---|
| 1422 | Fcn0[iadj][i] = MATMUL(TRANSPOSE(Jac0),Y[iadj][i]); |
---|
| 1423 | #else |
---|
| 1424 | JacTR_SP_Vec(Jac0,&Y[iadj][0],&Fcn0[iadj][0]); |
---|
| 1425 | #endif |
---|
| 1426 | } |
---|
| 1427 | |
---|
| 1428 | /*~~~> Repeat step calculation until current step accepted */ |
---|
| 1429 | do { /* UntilAccepted */ |
---|
| 1430 | |
---|
| 1431 | Singular = ros_PrepareMatrix(H,Direction,ros_Gamma[0], Jac0,Ghimj,Pivot, |
---|
| 1432 | ISTATUS); |
---|
| 1433 | if (Singular) /* More than 5 consecutive failed decompositions */ |
---|
| 1434 | return ros_ErrorMsg(-8,T,H); |
---|
| 1435 | |
---|
| 1436 | /*~~~> Compute the stages */ |
---|
| 1437 | for ( istage = 0; istage < ros_S; istage++ ) { /* Stage loop */ |
---|
| 1438 | |
---|
| 1439 | /* Current istage offset. Current istage vector |
---|
| 1440 | is K[ioffset+1:ioffset+NVAR] */ |
---|
| 1441 | ioffset = NVAR*(istage-1); |
---|
| 1442 | |
---|
| 1443 | /* For the 1st istage the function has been computed previously */ |
---|
| 1444 | if ( istage == 0 ) { |
---|
| 1445 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1446 | WCOPY(NVAR,&Fcn0[iadj][0],1,&Fcn[iadj][0],1); |
---|
| 1447 | |
---|
| 1448 | /* istage>0 and a new function evaluation is needed at |
---|
| 1449 | the current istage */ |
---|
| 1450 | } |
---|
| 1451 | else if ( ros_NewF[istage] ) { |
---|
| 1452 | WCOPY(NVAR*NADJ,&Y[0][0],1,&Ynew[0][0],1); |
---|
| 1453 | for (j = 0; j < istage-1; j++) { |
---|
| 1454 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1455 | WAXPY( NVAR,ros_A[(istage-1)*(istage-2)/2+j], |
---|
| 1456 | &K[iadj][NVAR*(j-1)+1],1,&Ynew[iadj][0],1); |
---|
| 1457 | } /* End for loop */ |
---|
| 1458 | Tau = T + ros_Alpha[istage]*Direction*H; |
---|
| 1459 | ros_cadj_Y( Tau, Y0 ); |
---|
| 1460 | JacTemplate(Tau, Y0, Jac); |
---|
| 1461 | ISTATUS[Njac]++; |
---|
| 1462 | |
---|
| 1463 | #ifdef FULL_ALGEBRA |
---|
| 1464 | for(i=0; i<NVAR; i++) { |
---|
| 1465 | for(j=0; j<NVAR; j++) |
---|
| 1466 | Jac[i][j] = -Jac[i][j]; |
---|
| 1467 | } |
---|
| 1468 | #else |
---|
| 1469 | WSCAL(LU_NONZERO,(-ONE),Jac,1); |
---|
| 1470 | #endif |
---|
| 1471 | |
---|
| 1472 | for ( iadj = 0; iadj < NADJ; iadj++ ) { |
---|
| 1473 | #ifdef FULL_ALGEBRA |
---|
| 1474 | for(i=0; i<NVAR; i++) |
---|
| 1475 | Fcn[iadj][i] = MATMUL(TRANSPOSE(Jac),Ynew[iadj][i]); |
---|
| 1476 | #else |
---|
| 1477 | JacTR_SP_Vec(Jac,&Ynew[iadj][0],&Fcn[iadj][0]); |
---|
| 1478 | #endif |
---|
| 1479 | } /* End for loop */ |
---|
| 1480 | } /* if istage == 1 elseif ros_NewF(istage) */ |
---|
| 1481 | |
---|
| 1482 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1483 | WCOPY(NVAR,&Fcn[iadj][0],1,&K[iadj][ioffset+1],1); |
---|
| 1484 | for ( j = 0; j < istage-1; j++ ) { |
---|
| 1485 | HC = ros_C[(istage-1)*(istage-2)/2+j]/(Direction*H); |
---|
| 1486 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1487 | WAXPY(NVAR,HC,&K[iadj][NVAR*(j-1)+1],1,&K[iadj][ioffset+1],1); |
---|
| 1488 | } /* End for loop */ |
---|
| 1489 | if ((!Autonomous) && (ros_Gamma[istage] != ZERO)) { |
---|
| 1490 | HG = Direction*H*ros_Gamma[istage]; |
---|
| 1491 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1492 | WAXPY(NVAR,HG,&dFdT[iadj][0],1,&K[iadj][ioffset+1],1); |
---|
| 1493 | } /* End if */ |
---|
| 1494 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1495 | ros_Solve('T', Ghimj, Pivot, &K[iadj][ioffset+1], ISTATUS); |
---|
| 1496 | |
---|
| 1497 | } /* End of Stage loop */ |
---|
| 1498 | |
---|
| 1499 | /*~~~> Compute the new solution */ |
---|
| 1500 | for ( iadj = 0; iadj < NADJ; iadj++ ) { |
---|
| 1501 | WCOPY(NVAR,&Y[iadj][0],1,&Ynew[iadj][0],1); |
---|
| 1502 | for ( j=0; j<ros_S; j++ ) |
---|
| 1503 | WAXPY(NVAR,ros_M[j],&K[iadj][NVAR*(j-1)+1],1,&Ynew[iadj][0],1); |
---|
| 1504 | } /* End for loop */ |
---|
| 1505 | |
---|
| 1506 | /*~~~> Compute the error estimation */ |
---|
| 1507 | WSCAL(NVAR*NADJ,ZERO,&Yerr[0][0],1); |
---|
| 1508 | for ( j=0; j<ros_S; j++ ) { |
---|
| 1509 | for ( iadj = 0; iadj < NADJ; iadj++ ) |
---|
| 1510 | WAXPY(NVAR,ros_E[j],&K[iadj][NVAR*(j-1)+1],1,&Yerr[iadj][0],1); |
---|
| 1511 | } /* End for loop */ |
---|
| 1512 | |
---|
| 1513 | /*~~~> Max error among all adjoint components */ |
---|
| 1514 | iadj = 1; |
---|
| 1515 | Err = ros_ErrorNorm ( &Y[iadj][0], &Ynew[iadj][0], &Yerr[iadj][0], |
---|
| 1516 | &AbsTol_adj[iadj][0], &RelTol_adj[iadj][0], |
---|
| 1517 | VectorTol ); |
---|
| 1518 | |
---|
| 1519 | /*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */ |
---|
| 1520 | Fac = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,(ONE/ros_ELO)))); |
---|
| 1521 | Hnew = H*Fac; |
---|
| 1522 | |
---|
| 1523 | /*~~~> Check the error magnitude and adjust step size */ |
---|
| 1524 | /* ISTATUS[Nstp] = ISTATUS[Nstp] + 1 */ |
---|
| 1525 | if ( (Err <= ONE) || (H <= Hmin) ) { /*~~~> Accept step */ |
---|
| 1526 | ISTATUS[Nacc] = ISTATUS[Nacc] + 1; |
---|
| 1527 | WCOPY(NVAR*NADJ,&Ynew[0][0],1,&Y[0][0],1); |
---|
| 1528 | T = T + Direction*H; |
---|
| 1529 | Hnew = MAX(Hmin,MIN(Hnew,Hmax)); |
---|
| 1530 | if (RejectLastH) /* No step size increase after a rejected step */ |
---|
| 1531 | Hnew = MIN(Hnew,H); |
---|
| 1532 | RSTATUS[Nhexit] = H; |
---|
| 1533 | RSTATUS[Nhnew] = Hnew; |
---|
| 1534 | RSTATUS[Ntexit] = T; |
---|
| 1535 | RejectLastH = FALSE; |
---|
| 1536 | RejectMoreH = FALSE; |
---|
| 1537 | H = Hnew; |
---|
| 1538 | break; /* UntilAccepted - EXIT THE LOOP: WHILE STEP NOT ACCEPTED */ |
---|
| 1539 | } |
---|
| 1540 | else { /*~~~> Reject step */ |
---|
| 1541 | if (RejectMoreH) |
---|
| 1542 | Hnew = H*FacRej; |
---|
| 1543 | RejectMoreH = RejectLastH; |
---|
| 1544 | RejectLastH = TRUE; |
---|
| 1545 | H = Hnew; |
---|
| 1546 | if (ISTATUS[Nacc] >= 1) |
---|
| 1547 | ISTATUS[Nrej]++; |
---|
| 1548 | } /* Err <= 1 */ |
---|
| 1549 | |
---|
| 1550 | } while(1); /* End of UntilAccepted do loop */ |
---|
| 1551 | |
---|
| 1552 | } /* End of TimeLoop */ |
---|
| 1553 | |
---|
| 1554 | /*~~~> Succesful exit */ |
---|
| 1555 | return 1; /*~~~> The integration was successful */ |
---|
| 1556 | } /* End of ros_CadjInt */ |
---|
| 1557 | |
---|
| 1558 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1559 | int ros_SimpleCadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, |
---|
| 1560 | KPP_REAL Tend, KPP_REAL T, int ISTATUS[], |
---|
| 1561 | int Autonomous, KPP_REAL Roundoff ) { |
---|
| 1562 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1563 | Template for the implementation of a generic RosenbrockADJ method |
---|
| 1564 | defined by ros_S (no of stages) |
---|
| 1565 | and its coefficients ros_{A,C,M,E,Alpha,Gamma} |
---|
| 1566 | ~~~> NADJ, Y[NADJ][NVAR] - Input: the initial condition at Tstart; |
---|
| 1567 | Output: the solution at T |
---|
| 1568 | ~~~> Tstart, Tend - Input: integration interval |
---|
| 1569 | ~~~> T - Output: time at which the solution is returned (T=Tend if success) |
---|
| 1570 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1571 | |
---|
| 1572 | /*~~~~ Local variables */ |
---|
| 1573 | KPP_REAL Y0[NVAR]; |
---|
| 1574 | KPP_REAL Ynew[NADJ][NVAR], Fcn0[NADJ][NVAR], Fcn[NADJ][NVAR]; |
---|
| 1575 | KPP_REAL K[NADJ][NVAR*ros_S], dFdT[NADJ][NVAR]; |
---|
| 1576 | #ifdef FULL_ALGEBRA |
---|
| 1577 | KPP_REAL Jac0[NVAR][NVAR], Ghimj[NVAR][NVAR], Jac[NVAR][NVAR], |
---|
| 1578 | dJdT[NVAR][NVAR]; |
---|
| 1579 | #else |
---|
| 1580 | KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO], Jac[LU_NONZERO], |
---|
| 1581 | dJdT[LU_NONZERO]; |
---|
| 1582 | #endif |
---|
| 1583 | KPP_REAL H, HC, HG, Tau; |
---|
| 1584 | KPP_REAL ghinv; |
---|
| 1585 | int Pivot[NVAR], Direction, ioffset, i, j, istage, iadj; |
---|
| 1586 | int istack; |
---|
| 1587 | |
---|
| 1588 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1589 | |
---|
| 1590 | /*~~~> INITIAL PREPARATIONS */ |
---|
| 1591 | if (Tend >= Tstart) |
---|
| 1592 | Direction = -1; |
---|
| 1593 | else |
---|
| 1594 | Direction = 1; |
---|
| 1595 | |
---|
| 1596 | /*~~~> Time loop begins below */ |
---|
| 1597 | for( istack = stack_ptr; istack >= 1; istack-- ) { /* TimeLoop */ |
---|
| 1598 | T = chk_T[istack]; |
---|
| 1599 | H = chk_H[istack-1]; |
---|
| 1600 | for(i=0; i<NVAR; i++) |
---|
| 1601 | Y0[i] = chk_Y[istack][i]; |
---|
| 1602 | |
---|
| 1603 | /*~~~> Compute the Jacobian at current time */ |
---|
| 1604 | JacTemplate(T, Y0, Jac0); |
---|
| 1605 | ISTATUS[Njac] = ISTATUS[Njac] + 1; |
---|
| 1606 | |
---|
| 1607 | /*~~~> Compute the function derivative with respect to T */ |
---|
| 1608 | if (!Autonomous) { |
---|
| 1609 | ros_JacTimeDerivative ( T, Roundoff, Y0, Jac0, dJdT, ISTATUS ); |
---|
| 1610 | for ( iadj = 0; iadj < NADJ; iadj++ ) { |
---|
| 1611 | #ifdef FULL_ALGEBRA |
---|
| 1612 | for( i=0; i<NVAR; i++) |
---|
| 1613 | dFdT[iadj][i] = MATMUL(TRANSPOSE(dJdT),Y[iadj][i]); |
---|
| 1614 | #else |
---|
| 1615 | JacTR_SP_Vec(dJdT,&Y[iadj][0],&dFdT[iadj][0]); |
---|
| 1616 | #endif |
---|
| 1617 | WSCAL(NVAR,(-ONE),&dFdT[iadj][0],1); |
---|
| 1618 | } |
---|
| 1619 | } |
---|
| 1620 | |
---|
| 1621 | /*~~~> Ydot = -J^T*Y */ |
---|
| 1622 | #ifdef FULL_ALGEBRA |
---|
| 1623 | for(i=0; i<NVAR; i++) { |
---|
| 1624 | for(j=0; j<NVAR; j++) |
---|
| 1625 | Jac0[i][j] = -Jac0[i][j]; |
---|
| 1626 | } |
---|
| 1627 | #else |
---|
| 1628 | WSCAL(LU_NONZERO,(-ONE),Jac0,1); |
---|
| 1629 | #endif |
---|
| 1630 | |
---|
| 1631 | for(iadj=0; iadj<NADJ; iadj++) { |
---|
| 1632 | #ifdef FULL_ALGEBRA |
---|
| 1633 | for(i=0; i<NVAR; i++) |
---|
| 1634 | Fcn0[iadj][i] = MATMUL(TRANSPOSE(Jac0),Y[iadj][i]); |
---|
| 1635 | #else |
---|
| 1636 | JacTR_SP_Vec(Jac0,&Y[iadj][0],&Fcn0[iadj][0]); |
---|
| 1637 | #endif |
---|
| 1638 | } |
---|
| 1639 | |
---|
| 1640 | /*~~~> Construct Ghimj = 1/(H*ham) - Jac0 */ |
---|
| 1641 | ghinv = ONE/(Direction*H*ros_Gamma[0]); |
---|
| 1642 | #ifdef FULL_ALGEBRA |
---|
| 1643 | for(i=0; i<NVAR; i++) { |
---|
| 1644 | for(j=0; j<NVAR; j++) |
---|
| 1645 | Ghimj[i][j] = -Jac0[i][j]; |
---|
| 1646 | } |
---|
| 1647 | for(i=0; i<NVAR; i++) |
---|
| 1648 | Ghimj[i][i] = Ghimj[i][i]+ghinv; |
---|
| 1649 | #else |
---|
| 1650 | WCOPY(LU_NONZERO,Jac0,1,Ghimj,1); |
---|
| 1651 | WSCAL(LU_NONZERO,(-ONE),Ghimj,1); |
---|
| 1652 | for(i=0; i<NVAR; i++) |
---|
| 1653 | Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv; |
---|
| 1654 | #endif |
---|
| 1655 | |
---|
| 1656 | /*~~~> Compute LU decomposition */ |
---|
| 1657 | ros_Decomp( Ghimj, Pivot, &j, ISTATUS ); |
---|
| 1658 | if (j != 0) { |
---|
| 1659 | ros_ErrorMsg(-8,T,H); |
---|
| 1660 | printf( "The matrix is singular !"); |
---|
| 1661 | exit(0); |
---|
| 1662 | } |
---|
| 1663 | |
---|
| 1664 | /*~~~> Compute the stages */ |
---|
| 1665 | for(istage=0; istage<ros_S; istage++) { /* Stage */ |
---|
| 1666 | /* Current istage offset. Current istage vector |
---|
| 1667 | is K(ioffset+1:ioffset+NVAR) */ |
---|
| 1668 | ioffset = NVAR*istage; |
---|
| 1669 | |
---|
| 1670 | /* For the 1st istage the function has been computed previously */ |
---|
| 1671 | if ( istage == 0 ) { |
---|
| 1672 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1673 | WCOPY(NVAR,&Fcn0[iadj][0],1,&Fcn[iadj][0],1); |
---|
| 1674 | } |
---|
| 1675 | /* istage>=1 and a new function evaluation is needed |
---|
| 1676 | at the current istage */ |
---|
| 1677 | else if ( ros_NewF[istage] ) { |
---|
| 1678 | WCOPY(NVAR*NADJ,&Y[0][0],1,&Ynew[0][0],1); |
---|
| 1679 | for(j=0; j<istage; j++) { |
---|
| 1680 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1681 | WAXPY(NVAR,ros_A[istage*(istage-1)/2+j], &K[iadj][NVAR*j],1, |
---|
| 1682 | &Ynew[iadj][0],1); |
---|
| 1683 | } |
---|
| 1684 | |
---|
| 1685 | Tau = T + ros_Alpha[istage]*Direction*H; |
---|
| 1686 | for(i=0; i<NVAR; i++) |
---|
| 1687 | ros_Hermite3( chk_T[istack-1], chk_T[istack], Tau, |
---|
| 1688 | &chk_Y[istack-1][i], &chk_Y[istack][i], |
---|
| 1689 | &chk_dY[istack-1][i], &chk_dY[istack][i], Y0 ); |
---|
| 1690 | JacTemplate(Tau, Y0, Jac); |
---|
| 1691 | ISTATUS[Njac]++; |
---|
| 1692 | |
---|
| 1693 | #ifdef FULL_ALGEBRA |
---|
| 1694 | for(i=0; i<NVAR; i++) { |
---|
| 1695 | for(j=0; j<NVAR; j++) |
---|
| 1696 | Jac[i][j] = -Jac[i][j]; |
---|
| 1697 | } |
---|
| 1698 | #else |
---|
| 1699 | WSCAL(LU_NONZERO,(-ONE),Jac,1); |
---|
| 1700 | #endif |
---|
| 1701 | |
---|
| 1702 | for(iadj=0; iadj<NADJ; iadj++) { |
---|
| 1703 | #ifdef FULL_ALGEBRA |
---|
| 1704 | for(i=0; i<NVAR; i++) |
---|
| 1705 | Fcn[iadj][i] = MATMUL(TRANSPOSE(Jac),Ynew[iadj][i]); |
---|
| 1706 | #else |
---|
| 1707 | JacTR_SP_Vec(Jac,&Ynew[iadj][0],&Fcn[iadj][0]); |
---|
| 1708 | #endif |
---|
| 1709 | } |
---|
| 1710 | } /* if istage == 1 elseif ros_NewF(istage) */ |
---|
| 1711 | |
---|
| 1712 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1713 | WCOPY(NVAR,&Fcn[iadj][0],1,&K[iadj][ioffset],1); |
---|
| 1714 | for(j=0; j<istage-1; j++) { |
---|
| 1715 | HC = ros_C[istage*(istage-1)/2+j]/(Direction*H); |
---|
| 1716 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1717 | WAXPY(NVAR,HC,&K[iadj][NVAR*j],1,&K[iadj][ioffset],1); |
---|
| 1718 | } |
---|
| 1719 | if((!Autonomous) && (ros_Gamma[istage] != ZERO)) { |
---|
| 1720 | HG = Direction*H*ros_Gamma[istage]; |
---|
| 1721 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1722 | WAXPY(NVAR,HG,&dFdT[iadj][0],1,&K[iadj][ioffset],1); |
---|
| 1723 | } |
---|
| 1724 | for(iadj=0; iadj<NADJ; iadj++) |
---|
| 1725 | ros_Solve('T', Ghimj, Pivot, &K[iadj][ioffset], ISTATUS); |
---|
| 1726 | } /* End of Stage loop */ |
---|
| 1727 | |
---|
| 1728 | /*~~~> Compute the new solution */ |
---|
| 1729 | for(iadj=0; iadj<NADJ; iadj++) { |
---|
| 1730 | for(j=0; j<ros_S; j++) |
---|
| 1731 | WAXPY(NVAR,ros_M[j],&K[iadj][NVAR*j],1,&Y[iadj][0],1); |
---|
| 1732 | } |
---|
| 1733 | } /* End of TimeLoop */ |
---|
| 1734 | |
---|
| 1735 | /*~~~> Succesful exit */ |
---|
| 1736 | return 1; /*~~~> The integration was successful */ |
---|
| 1737 | |
---|
| 1738 | } /* End of ros_SimpleCadjInt */ |
---|
| 1739 | |
---|
| 1740 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1741 | KPP_REAL ros_ErrorNorm ( KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], |
---|
| 1742 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int VectorTol ) { |
---|
| 1743 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1744 | ~~~> Computes the "scaled norm" of the error vector Yerr |
---|
| 1745 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1746 | |
---|
| 1747 | /* Local variables */ |
---|
| 1748 | KPP_REAL Err, Scale, Ymax; |
---|
| 1749 | int i; |
---|
| 1750 | |
---|
| 1751 | Err = ZERO; |
---|
| 1752 | for(i=0; i<NVAR; i++) { |
---|
| 1753 | Ymax = MAX(ABS(Y[i]),ABS(Ynew[i])); |
---|
| 1754 | if (VectorTol) |
---|
| 1755 | Scale = AbsTol[i]+RelTol[i]*Ymax; |
---|
| 1756 | else |
---|
| 1757 | Scale = AbsTol[0]+RelTol[0]*Ymax; |
---|
| 1758 | |
---|
| 1759 | Err = Err+pow((Yerr[i]/Scale),2); |
---|
| 1760 | } |
---|
| 1761 | Err = SQRT(Err/NVAR); |
---|
| 1762 | |
---|
| 1763 | return MAX(Err,(KPP_REAL)1.0e-10); |
---|
| 1764 | } /* End of ros_ErrorNorm */ |
---|
| 1765 | |
---|
| 1766 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1767 | void ros_FunTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], |
---|
| 1768 | KPP_REAL Fcn0[], KPP_REAL dFdT[], int ISTATUS[]) { |
---|
| 1769 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1770 | ~~~> The time partial derivative of the function by finite differences |
---|
| 1771 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1772 | |
---|
| 1773 | /*~~~> Local variables */ |
---|
| 1774 | KPP_REAL Delta; |
---|
| 1775 | |
---|
| 1776 | Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)); |
---|
| 1777 | FunTemplate(T+Delta,Y,dFdT); |
---|
| 1778 | ISTATUS[Nfun]++; |
---|
| 1779 | WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1); |
---|
| 1780 | WSCAL(NVAR,(ONE/Delta),dFdT,1); |
---|
| 1781 | |
---|
| 1782 | } /* End of ros_FunTimeDerivative */ |
---|
| 1783 | |
---|
| 1784 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1785 | void ros_JacTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], |
---|
| 1786 | KPP_REAL Jac0[], KPP_REAL dJdT[], int ISTATUS[]) { |
---|
| 1787 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1788 | ~~~> The time partial derivative of the Jacobian by finite differences |
---|
| 1789 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1790 | |
---|
| 1791 | /*~~~> Local variables */ |
---|
| 1792 | KPP_REAL Delta; |
---|
| 1793 | |
---|
| 1794 | Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)); |
---|
| 1795 | JacTemplate(T+Delta,Y,dJdT); |
---|
| 1796 | ISTATUS[Njac]++; |
---|
| 1797 | #ifdef FULL_ALGEBRA |
---|
| 1798 | WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1); |
---|
| 1799 | WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1); |
---|
| 1800 | #else |
---|
| 1801 | WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1); |
---|
| 1802 | WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1); |
---|
| 1803 | #endif |
---|
| 1804 | } /* End of ros_JacTimeDerivative */ |
---|
| 1805 | |
---|
| 1806 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1807 | int ros_PrepareMatrix ( KPP_REAL H, int Direction, KPP_REAL gam, |
---|
| 1808 | KPP_REAL Jac0[], KPP_REAL Ghimj[], int Pivot[], |
---|
| 1809 | int ISTATUS[] ) { |
---|
| 1810 | /* --- --- --- --- --- --- --- --- --- --- --- --- --- |
---|
| 1811 | Prepares the LHS matrix for stage calculations |
---|
| 1812 | 1. Construct Ghimj = 1/(H*gam) - Jac0 |
---|
| 1813 | "(Gamma H) Inverse Minus Jacobian" |
---|
| 1814 | 2. Repeat LU decomposition of Ghimj until successful. |
---|
| 1815 | -half the step size if LU decomposition fails and retry |
---|
| 1816 | -exit after 5 consecutive fails |
---|
| 1817 | --- --- --- --- --- --- --- --- --- --- --- --- --- */ |
---|
| 1818 | |
---|
| 1819 | /*~~~> Local variables */ |
---|
| 1820 | int i, ising, Nconsecutive; |
---|
| 1821 | int Singular; /* Boolean value */ |
---|
| 1822 | KPP_REAL ghinv; |
---|
| 1823 | |
---|
| 1824 | Nconsecutive = 0; |
---|
| 1825 | Singular = TRUE; |
---|
| 1826 | |
---|
| 1827 | while (Singular) { |
---|
| 1828 | |
---|
| 1829 | /*~~~> Construct Ghimj = 1/(H*gam) - Jac0 */ |
---|
| 1830 | #ifdef FULL_ALGEBRA |
---|
| 1831 | WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1); |
---|
| 1832 | WSCAL(NVAR*NVAR,(-ONE),Ghimj,1); |
---|
| 1833 | ghinv = ONE/(Direction*H*gam); |
---|
| 1834 | for(i=0; i<NVAR; i++) |
---|
| 1835 | Ghimj[i][i] = Ghimj[i][i]+ghinv; |
---|
| 1836 | #else |
---|
| 1837 | WCOPY(LU_NONZERO,Jac0,1,Ghimj,1); |
---|
| 1838 | WSCAL(LU_NONZERO,(-ONE),Ghimj,1); |
---|
| 1839 | ghinv = ONE/(Direction*H*gam); |
---|
| 1840 | for(i=0; i<NVAR; i++) |
---|
| 1841 | Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv; |
---|
| 1842 | #endif |
---|
| 1843 | |
---|
| 1844 | /*~~~> Compute LU decomposition */ |
---|
| 1845 | ros_Decomp( Ghimj, Pivot, &ising, ISTATUS ); |
---|
| 1846 | if (ising == 0) |
---|
| 1847 | /*~~~> If successful done */ |
---|
| 1848 | Singular = FALSE; |
---|
| 1849 | |
---|
| 1850 | else { /* ising != 0 */ |
---|
| 1851 | /*~~~> If unsuccessful half the step size; |
---|
| 1852 | if 5 consecutive fails then return */ |
---|
| 1853 | ISTATUS[Nsng]++; |
---|
| 1854 | Nconsecutive++; |
---|
| 1855 | Singular = TRUE; |
---|
| 1856 | printf( "Warning: LU Decomposition returned ising = %d", ising ); |
---|
| 1857 | if (Nconsecutive <= 5) /*Less than 5 consecutive failed decompositions*/ |
---|
| 1858 | H = H*HALF; |
---|
| 1859 | else /* More than 5 consecutive failed decompositions */ |
---|
| 1860 | return Singular; |
---|
| 1861 | |
---|
| 1862 | } /* End of ising */ |
---|
| 1863 | } /* while Singular */ |
---|
| 1864 | |
---|
| 1865 | return Singular; |
---|
| 1866 | } /* End of ros_PrepareMatrix */ |
---|
| 1867 | |
---|
| 1868 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1869 | void ros_Decomp( KPP_REAL A[], int Pivot[], int* ising, int ISTATUS[] ) { |
---|
| 1870 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1871 | Template for the LU decomposition |
---|
| 1872 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1873 | |
---|
| 1874 | #ifdef FULL_ALGEBRA |
---|
| 1875 | DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising ); |
---|
| 1876 | #else |
---|
| 1877 | *ising = KppDecomp ( A ); |
---|
| 1878 | Pivot[0] = 1; |
---|
| 1879 | #endif |
---|
| 1880 | ISTATUS[Ndec]++; |
---|
| 1881 | |
---|
| 1882 | } /* End of ros_Decomp */ |
---|
| 1883 | |
---|
| 1884 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1885 | void ros_Solve( char How, KPP_REAL A[], int Pivot[], KPP_REAL b[], |
---|
| 1886 | int ISTATUS[] ) { |
---|
| 1887 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1888 | Template for the forward/backward substitution |
---|
| 1889 | (using pre-computed LU decomposition) |
---|
| 1890 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1891 | switch (How) { |
---|
| 1892 | case 'N': |
---|
| 1893 | #ifdef FULL_ALGEBRA |
---|
| 1894 | DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 ); |
---|
| 1895 | #else |
---|
| 1896 | KppSolve( A, b ); |
---|
| 1897 | #endif |
---|
| 1898 | break; |
---|
| 1899 | case 'T': |
---|
| 1900 | #ifdef FULL_ALGEBRA |
---|
| 1901 | DGETRS( 'T', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 ); |
---|
| 1902 | #else |
---|
| 1903 | KppSolveTR( A, b, b ); |
---|
| 1904 | #endif |
---|
| 1905 | break; |
---|
| 1906 | default: |
---|
| 1907 | printf( "Error: unknown argument in ros_Solve: How=%d", How ); |
---|
| 1908 | exit(0); |
---|
| 1909 | } |
---|
| 1910 | ISTATUS[Nsol]++; |
---|
| 1911 | |
---|
| 1912 | } /* End of ros_Solve */ |
---|
| 1913 | |
---|
| 1914 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1915 | void ros_cadj_Y( KPP_REAL T, KPP_REAL Y[] ) { |
---|
| 1916 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1917 | Finds the solution Y at T by interpolating the stored forward trajectory |
---|
| 1918 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1919 | |
---|
| 1920 | /*~~~> Local variables */ |
---|
| 1921 | int i,j; |
---|
| 1922 | |
---|
| 1923 | if( (T < chk_T[0]) || (T> chk_T[stack_ptr]) ) { |
---|
| 1924 | printf( "Cannot locate solution at T = %f", T ); |
---|
| 1925 | printf( "Stored trajectory is between Tstart = %f", chk_T[0] ); |
---|
| 1926 | printf( " and Tend = %f", chk_T[stack_ptr] ); |
---|
| 1927 | exit(0); |
---|
| 1928 | } |
---|
| 1929 | for(i=0; i<stack_ptr-1; i++) { |
---|
| 1930 | if( (T >= chk_T[i]) &&(T <= chk_T[i+1]) ) |
---|
| 1931 | exit(0); |
---|
| 1932 | } |
---|
| 1933 | |
---|
| 1934 | for(j=0; j<NVAR; j++ ) |
---|
| 1935 | ros_Hermite3( chk_T[i], chk_T[i+1], T, &chk_Y[i][j], &chk_Y[i+1][j], |
---|
| 1936 | &chk_dY[i][j], &chk_dY[i+1][j], Y ); |
---|
| 1937 | |
---|
| 1938 | } /* End of ros_cadj_Y */ |
---|
| 1939 | |
---|
| 1940 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1941 | void ros_Hermite3( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], |
---|
| 1942 | KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Y[]) { |
---|
| 1943 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1944 | Template for Hermite interpolation of order 5 on the interval [a,b] |
---|
| 1945 | P = c(1) + c(2)*(x-a) + ... + c(4)*(x-a)^3 |
---|
| 1946 | P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb] |
---|
| 1947 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1948 | |
---|
| 1949 | /*~~~> Local variables */ |
---|
| 1950 | KPP_REAL Tau, amb[3], C[4][NVAR]; |
---|
| 1951 | int i, j; |
---|
| 1952 | |
---|
| 1953 | amb[0] = ((KPP_REAL)1.0)/(a-b); |
---|
| 1954 | for(i=1; i<3; i++) |
---|
| 1955 | amb[i] = amb[i-1]*amb[0]; |
---|
| 1956 | |
---|
| 1957 | /* c(1) = ya; */ |
---|
| 1958 | WCOPY(NVAR,Ya,1,&C[0][0],1); |
---|
| 1959 | /* c(2) = ja; */ |
---|
| 1960 | WCOPY(NVAR,Ja,1,&C[1][0],1); |
---|
| 1961 | /* c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb; */ |
---|
| 1962 | WCOPY(NVAR,Ya,1,&C[2][0],1); |
---|
| 1963 | WSCAL(NVAR,-3.0*amb[1],&C[2][0],1); |
---|
| 1964 | WAXPY(NVAR,3.0*amb[1],Yb,1,&C[2][0],1); |
---|
| 1965 | WAXPY(NVAR,2.0*amb[0],Ja,1,&C[2][0],1); |
---|
| 1966 | WAXPY(NVAR,amb[0],Jb,1,&C[2][0],1); |
---|
| 1967 | /* c(4) = 1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb */ |
---|
| 1968 | WCOPY(NVAR,Ya,1,&C[3][0],1); |
---|
| 1969 | WSCAL(NVAR,-2.0*amb[2],&C[3][0],1); |
---|
| 1970 | WAXPY(NVAR,2.0*amb[2],Yb,1,&C[3][0],1); |
---|
| 1971 | WAXPY(NVAR,amb[1],Ja,1,&C[3][0],1); |
---|
| 1972 | WAXPY(NVAR,amb[1],Jb,1,&C[3][0],1); |
---|
| 1973 | |
---|
| 1974 | Tau = T - a; |
---|
| 1975 | WCOPY(NVAR,&C[3][0],1,Y,1); |
---|
| 1976 | WSCAL(NVAR,pow(Tau,3),Y,1); |
---|
| 1977 | for(j=2; j>=0; j--) |
---|
| 1978 | WAXPY(NVAR,pow(Tau,(j-1)),&C[j][0],1,Y,1); |
---|
| 1979 | |
---|
| 1980 | } /* End of ros_Hermite3 */ |
---|
| 1981 | |
---|
| 1982 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1983 | void ros_Hermite5( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], |
---|
| 1984 | KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Ha[], |
---|
| 1985 | KPP_REAL Hb[], KPP_REAL Y[] ) { |
---|
| 1986 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1987 | Template for Hermite interpolation of order 5 on the interval [a,b] |
---|
| 1988 | P = c(1) + c(2)*(x-a) + ... + c(6)*(x-a)^5 |
---|
| 1989 | P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb] |
---|
| 1990 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1991 | |
---|
| 1992 | /*~~~> Local variables */ |
---|
| 1993 | KPP_REAL Tau, amb[5], C[6][NVAR]; |
---|
| 1994 | int i, j; |
---|
| 1995 | |
---|
| 1996 | amb[0] = ((KPP_REAL)1.0)/(a-b); |
---|
| 1997 | for(i=1; i<5; i++) |
---|
| 1998 | amb[i] = amb[i-1]*amb[0]; |
---|
| 1999 | |
---|
| 2000 | /* c(1) = ya; */ |
---|
| 2001 | WCOPY(NVAR,Ya,1,&C[0][0],1); |
---|
| 2002 | /* c(2) = ja; */ |
---|
| 2003 | WCOPY(NVAR,Ja,1,&C[1][0],1); |
---|
| 2004 | /* c(3) = ha/2; */ |
---|
| 2005 | WCOPY(NVAR,Ha,1,&C[2][0],1); |
---|
| 2006 | WSCAL(NVAR,HALF,&C[2][0],1); |
---|
| 2007 | |
---|
| 2008 | /* c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb |
---|
| 2009 | + 1.5*amb(1)*ha - 0.5*amb(1)*hb ; */ |
---|
| 2010 | WCOPY(NVAR,Ya,1,&C[3][0],1); |
---|
| 2011 | WSCAL(NVAR,10.0*amb[2],&C[3][0],1); |
---|
| 2012 | WAXPY(NVAR,-10.0*amb[2],Yb,1,&C[3][0],1); |
---|
| 2013 | WAXPY(NVAR,-6.0*amb[1],Ja,1,&C[3][0],1); |
---|
| 2014 | WAXPY(NVAR,-4.0*amb[1],Jb,1,&C[3][0],1); |
---|
| 2015 | WAXPY(NVAR, 1.5*amb[0],Ha,1,&C[3][0],1); |
---|
| 2016 | WAXPY(NVAR,-0.5*amb[0],Hb,1,&C[3][0],1); |
---|
| 2017 | |
---|
| 2018 | /* c(5) = 15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb |
---|
| 2019 | + 1.5*amb(2)*ha - 1*amb(2)*hb ; */ |
---|
| 2020 | WCOPY(NVAR,Ya,1,&C[4][0],1); |
---|
| 2021 | WSCAL(NVAR, 15.0*amb[3],&C[4][0],1); |
---|
| 2022 | WAXPY(NVAR,-15.0*amb[3],Yb,1,&C[4][0],1); |
---|
| 2023 | WAXPY(NVAR,-8.0*amb[2],Ja,1,&C[4][0],1); |
---|
| 2024 | WAXPY(NVAR,-7.0*amb[2],Jb,1,&C[4][0],1); |
---|
| 2025 | WAXPY(NVAR,1.5*amb[1],Ha,1,&C[4][0],1); |
---|
| 2026 | WAXPY(NVAR,-amb[1],Hb,1,&C[4][0],1); |
---|
| 2027 | |
---|
| 2028 | /* c(6) = 6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb |
---|
| 2029 | + 0.5*amb(3)*ha -0.5*amb(3)*hb ; */ |
---|
| 2030 | WCOPY(NVAR,Ya,1,&C[5][0],1); |
---|
| 2031 | WSCAL(NVAR, 6.0*amb[4],&C[5][0],1); |
---|
| 2032 | WAXPY(NVAR,-6.0*amb[4],Yb,1,&C[5][0],1); |
---|
| 2033 | WAXPY(NVAR,-3.0*amb[3],Ja,1,&C[5][0],1); |
---|
| 2034 | WAXPY(NVAR,-3.0*amb[3],Jb,1,&C[5][0],1); |
---|
| 2035 | WAXPY(NVAR, 0.5*amb[2],Ha,1,&C[5][0],1); |
---|
| 2036 | WAXPY(NVAR,-0.5*amb[2],Hb,1,&C[5][0],1); |
---|
| 2037 | |
---|
| 2038 | Tau = T - a; |
---|
| 2039 | WCOPY(NVAR,&C[5][0],1,Y,1); |
---|
| 2040 | for(j=4; j>=0; j--) { |
---|
| 2041 | WSCAL(NVAR,Tau,Y,1); |
---|
| 2042 | WAXPY(NVAR,ONE,&C[j][0],1,Y,1); |
---|
| 2043 | } |
---|
| 2044 | |
---|
| 2045 | } /* End of ros_Hermite5 */ |
---|
| 2046 | |
---|
| 2047 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2048 | void Ros2() { |
---|
| 2049 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2050 | --- AN L-STABLE METHOD, 2 stages, order 2 |
---|
| 2051 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2052 | |
---|
| 2053 | KPP_REAL g; |
---|
| 2054 | |
---|
| 2055 | g = (KPP_REAL)1.0 + ((KPP_REAL)1.0)/SQRT((KPP_REAL)2.0); |
---|
| 2056 | |
---|
| 2057 | rosMethod = RS2; |
---|
| 2058 | /*~~~> Name of the method */ |
---|
| 2059 | strcpy(ros_Name, "ROS-2"); |
---|
| 2060 | /*~~~> Number of stages */ |
---|
| 2061 | ros_S = 2; |
---|
| 2062 | |
---|
| 2063 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 2064 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 2065 | A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc. |
---|
| 2066 | The general mapping formula is: |
---|
| 2067 | A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ] |
---|
| 2068 | C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */ |
---|
| 2069 | |
---|
| 2070 | ros_A[0] = ((KPP_REAL)1.0)/g; |
---|
| 2071 | ros_C[0] = ((KPP_REAL)-2.0)/g; |
---|
| 2072 | |
---|
| 2073 | /*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE) |
---|
| 2074 | or does it re-use the function evaluation from stage i-1 |
---|
| 2075 | (ros_NewF[i]=FALSE) */ |
---|
| 2076 | ros_NewF[0] = TRUE; |
---|
| 2077 | ros_NewF[1] = TRUE; |
---|
| 2078 | |
---|
| 2079 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 2080 | ros_M[0]= ((KPP_REAL)3.0)/((KPP_REAL)2.0*g); |
---|
| 2081 | ros_M[1]= ((KPP_REAL)1.0)/((KPP_REAL)2.0*g); |
---|
| 2082 | |
---|
| 2083 | /* E_i = Coefficients for error estimator */ |
---|
| 2084 | ros_E[0] = ((KPP_REAL)1.0)/((KPP_REAL)2.0*g); |
---|
| 2085 | ros_E[1] = ((KPP_REAL)1.0)/((KPP_REAL)2.0*g); |
---|
| 2086 | |
---|
| 2087 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 2088 | main and the embedded scheme orders plus one */ |
---|
| 2089 | ros_ELO = (KPP_REAL)2.0; |
---|
| 2090 | |
---|
| 2091 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 2092 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 2093 | ros_Alpha[1] = (KPP_REAL)1.0; |
---|
| 2094 | |
---|
| 2095 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 2096 | ros_Gamma[0] = g; |
---|
| 2097 | ros_Gamma[1] =-g; |
---|
| 2098 | |
---|
| 2099 | } /* End of Ros2 */ |
---|
| 2100 | |
---|
| 2101 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2102 | void Ros3() { |
---|
| 2103 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2104 | --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations |
---|
| 2105 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2106 | |
---|
| 2107 | rosMethod = RS3; |
---|
| 2108 | /*~~~> Name of the method */ |
---|
| 2109 | strcpy(ros_Name, "ROS-3"); |
---|
| 2110 | /*~~~> Number of stages */ |
---|
| 2111 | ros_S = 3; |
---|
| 2112 | |
---|
| 2113 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 2114 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 2115 | A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc. |
---|
| 2116 | The general mapping formula is: |
---|
| 2117 | A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ] |
---|
| 2118 | C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */ |
---|
| 2119 | |
---|
| 2120 | ros_A[0]= (KPP_REAL)1.0; |
---|
| 2121 | ros_A[1]= (KPP_REAL)1.0; |
---|
| 2122 | ros_A[2]= (KPP_REAL)0.0; |
---|
| 2123 | ros_C[0] = (KPP_REAL)-0.10156171083877702091975600115545e01; |
---|
| 2124 | ros_C[1] = (KPP_REAL) 0.40759956452537699824805835358067e01; |
---|
| 2125 | ros_C[2] = (KPP_REAL) 0.92076794298330791242156818474003e01; |
---|
| 2126 | |
---|
| 2127 | /*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE) |
---|
| 2128 | or does it re-use the function evaluation from stage i-1 |
---|
| 2129 | (ros_NewF[i]=FALSE) */ |
---|
| 2130 | ros_NewF[0] = TRUE; |
---|
| 2131 | ros_NewF[1] = TRUE; |
---|
| 2132 | ros_NewF[2] = FALSE; |
---|
| 2133 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 2134 | ros_M[0] = (KPP_REAL) 0.1e01; |
---|
| 2135 | ros_M[1] = (KPP_REAL) 0.61697947043828245592553615689730e01; |
---|
| 2136 | ros_M[2] = (KPP_REAL)-0.42772256543218573326238373806514; |
---|
| 2137 | /* E_i = Coefficients for error estimator */ |
---|
| 2138 | ros_E[0] = (KPP_REAL) 0.5; |
---|
| 2139 | ros_E[1] = (KPP_REAL)-0.29079558716805469821718236208017e01; |
---|
| 2140 | ros_E[2] = (KPP_REAL) 0.22354069897811569627360909276199; |
---|
| 2141 | |
---|
| 2142 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 2143 | main and the embedded scheme orders plus 1 */ |
---|
| 2144 | ros_ELO = (KPP_REAL)3.0; |
---|
| 2145 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 2146 | ros_Alpha[0]= (KPP_REAL)0.0; |
---|
| 2147 | ros_Alpha[1]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 2148 | ros_Alpha[2]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 2149 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 2150 | ros_Gamma[0]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 2151 | ros_Gamma[1]= (KPP_REAL)0.24291996454816804366592249683314; |
---|
| 2152 | ros_Gamma[2]= (KPP_REAL)0.21851380027664058511513169485832e01; |
---|
| 2153 | |
---|
| 2154 | } /* End of Ros3 */ |
---|
| 2155 | |
---|
| 2156 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2157 | void Ros4() { |
---|
| 2158 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2159 | L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES |
---|
| 2160 | L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 |
---|
| 2161 | |
---|
| 2162 | E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
| 2163 | EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
| 2164 | SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, |
---|
| 2165 | SPRINGER-VERLAG (1990) |
---|
| 2166 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2167 | |
---|
| 2168 | rosMethod = RS4; |
---|
| 2169 | /*~~~> Name of the method */ |
---|
| 2170 | strcpy(ros_Name, "ROS-4"); |
---|
| 2171 | /*~~~> Number of stages */ |
---|
| 2172 | ros_S = 4; |
---|
| 2173 | |
---|
| 2174 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 2175 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 2176 | A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc. |
---|
| 2177 | The general mapping formula is: |
---|
| 2178 | A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ] |
---|
| 2179 | C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */ |
---|
| 2180 | |
---|
| 2181 | ros_A[0] = (KPP_REAL)0.2000000000000000e01; |
---|
| 2182 | ros_A[1] = (KPP_REAL)0.1867943637803922e01; |
---|
| 2183 | ros_A[2] = (KPP_REAL)0.2344449711399156; |
---|
| 2184 | ros_A[3] = ros_A[1]; |
---|
| 2185 | ros_A[4] = ros_A[2]; |
---|
| 2186 | ros_A[5] = (KPP_REAL)0.0; |
---|
| 2187 | |
---|
| 2188 | ros_C[0] = (KPP_REAL)-0.7137615036412310e01; |
---|
| 2189 | ros_C[1] = (KPP_REAL) 0.2580708087951457e01; |
---|
| 2190 | ros_C[2] = (KPP_REAL) 0.6515950076447975; |
---|
| 2191 | ros_C[3] = (KPP_REAL)-0.2137148994382534e01; |
---|
| 2192 | ros_C[4] = (KPP_REAL)-0.3214669691237626; |
---|
| 2193 | ros_C[5] = (KPP_REAL)-0.6949742501781779; |
---|
| 2194 | |
---|
| 2195 | /*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE) |
---|
| 2196 | or does it re-use the function evaluation from stage i-1 |
---|
| 2197 | (ros_NewF[i]=FALSE) */ |
---|
| 2198 | ros_NewF[0] = TRUE; |
---|
| 2199 | ros_NewF[1] = TRUE; |
---|
| 2200 | ros_NewF[2] = TRUE; |
---|
| 2201 | ros_NewF[3] = FALSE; |
---|
| 2202 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 2203 | ros_M[0] = (KPP_REAL)0.2255570073418735e01; |
---|
| 2204 | ros_M[1] = (KPP_REAL)0.2870493262186792; |
---|
| 2205 | ros_M[2] = (KPP_REAL)0.4353179431840180; |
---|
| 2206 | ros_M[3] = (KPP_REAL)0.1093502252409163e01; |
---|
| 2207 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 2208 | ros_E[0] = (KPP_REAL)-0.2815431932141155; |
---|
| 2209 | ros_E[1] = (KPP_REAL)-0.7276199124938920e-01; |
---|
| 2210 | ros_E[2] = (KPP_REAL)-0.1082196201495311; |
---|
| 2211 | ros_E[3] = (KPP_REAL)-0.1093502252409163e01; |
---|
| 2212 | |
---|
| 2213 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 2214 | main and the embedded scheme orders plus 1 */ |
---|
| 2215 | ros_ELO = (KPP_REAL)4.0; |
---|
| 2216 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 2217 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 2218 | ros_Alpha[1] = (KPP_REAL)0.1145640000000000e01; |
---|
| 2219 | ros_Alpha[2] = (KPP_REAL)0.6552168638155900; |
---|
| 2220 | ros_Alpha[3] = ros_Alpha[2]; |
---|
| 2221 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 2222 | ros_Gamma[0] = (KPP_REAL) 0.5728200000000000; |
---|
| 2223 | ros_Gamma[1] = (KPP_REAL)-0.1769193891319233e01; |
---|
| 2224 | ros_Gamma[2] = (KPP_REAL) 0.7592633437920482; |
---|
| 2225 | ros_Gamma[3] = (KPP_REAL)-0.1049021087100450; |
---|
| 2226 | |
---|
| 2227 | } /* End of Ros4 */ |
---|
| 2228 | |
---|
| 2229 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2230 | void Rodas3() { |
---|
| 2231 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2232 | --- A STIFFLY-STABLE METHOD, 4 stages, order 3 |
---|
| 2233 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2234 | |
---|
| 2235 | rosMethod = RD3; |
---|
| 2236 | /*~~~> Name of the method */ |
---|
| 2237 | strcpy(ros_Name, "RODAS-3"); |
---|
| 2238 | /*~~~> Number of stages */ |
---|
| 2239 | ros_S = 4; |
---|
| 2240 | |
---|
| 2241 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 2242 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 2243 | A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc. |
---|
| 2244 | The general mapping formula is: |
---|
| 2245 | A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ] |
---|
| 2246 | C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */ |
---|
| 2247 | |
---|
| 2248 | ros_A[0] = (KPP_REAL)0.0; |
---|
| 2249 | ros_A[1] = (KPP_REAL)2.0; |
---|
| 2250 | ros_A[2] = (KPP_REAL)0.0; |
---|
| 2251 | ros_A[3] = (KPP_REAL)2.0; |
---|
| 2252 | ros_A[4] = (KPP_REAL)0.0; |
---|
| 2253 | ros_A[5] = (KPP_REAL)1.0; |
---|
| 2254 | |
---|
| 2255 | ros_C[0] = (KPP_REAL) 4.0; |
---|
| 2256 | ros_C[1] = (KPP_REAL) 1.0; |
---|
| 2257 | ros_C[2] = (KPP_REAL)-1.0; |
---|
| 2258 | ros_C[3] = (KPP_REAL) 1.0; |
---|
| 2259 | ros_C[4] = (KPP_REAL)-1.0; |
---|
| 2260 | ros_C[5] = -(((KPP_REAL)8.0)/((KPP_REAL)3.0)); |
---|
| 2261 | |
---|
| 2262 | /*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE) |
---|
| 2263 | or does it re-use the function evaluation from stage i-1 |
---|
| 2264 | (ros_NewF[i]=FALSE) */ |
---|
| 2265 | ros_NewF[0] = TRUE; |
---|
| 2266 | ros_NewF[1] = FALSE; |
---|
| 2267 | ros_NewF[2] = TRUE; |
---|
| 2268 | ros_NewF[3] = TRUE; |
---|
| 2269 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 2270 | ros_M[0] = (KPP_REAL)2.0; |
---|
| 2271 | ros_M[1] = (KPP_REAL)0.0; |
---|
| 2272 | ros_M[2] = (KPP_REAL)1.0; |
---|
| 2273 | ros_M[3] = (KPP_REAL)1.0; |
---|
| 2274 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 2275 | ros_E[0] = (KPP_REAL)0.0; |
---|
| 2276 | ros_E[1] = (KPP_REAL)0.0; |
---|
| 2277 | ros_E[2] = (KPP_REAL)0.0; |
---|
| 2278 | ros_E[3] = (KPP_REAL)1.0; |
---|
| 2279 | |
---|
| 2280 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 2281 | main and the embedded scheme orders plus 1 */ |
---|
| 2282 | ros_ELO = (KPP_REAL)3.0; |
---|
| 2283 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 2284 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 2285 | ros_Alpha[1] = (KPP_REAL)0.0; |
---|
| 2286 | ros_Alpha[2] = (KPP_REAL)1.0; |
---|
| 2287 | ros_Alpha[3] = (KPP_REAL)1.0; |
---|
| 2288 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 2289 | ros_Gamma[0] = (KPP_REAL)0.5; |
---|
| 2290 | ros_Gamma[1] = (KPP_REAL)1.5; |
---|
| 2291 | ros_Gamma[2] = (KPP_REAL)0.0; |
---|
| 2292 | ros_Gamma[3] = (KPP_REAL)0.0; |
---|
| 2293 | |
---|
| 2294 | } /* End of Rodas3 */ |
---|
| 2295 | |
---|
| 2296 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2297 | void Rodas4() { |
---|
| 2298 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2299 | STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES |
---|
| 2300 | |
---|
| 2301 | E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
| 2302 | EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
| 2303 | SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, |
---|
| 2304 | SPRINGER-VERLAG (1996) |
---|
| 2305 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2306 | |
---|
| 2307 | rosMethod = RD4; |
---|
| 2308 | /*~~~> Name of the method */ |
---|
| 2309 | strcpy(ros_Name, "RODAS-4"); |
---|
| 2310 | /*~~~> Number of stages */ |
---|
| 2311 | ros_S = 6; |
---|
| 2312 | |
---|
| 2313 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 2314 | ros_Alpha[0] = (KPP_REAL)0.000; |
---|
| 2315 | ros_Alpha[1] = (KPP_REAL)0.386; |
---|
| 2316 | ros_Alpha[2] = (KPP_REAL)0.210; |
---|
| 2317 | ros_Alpha[3] = (KPP_REAL)0.630; |
---|
| 2318 | ros_Alpha[4] = (KPP_REAL)1.000; |
---|
| 2319 | ros_Alpha[5] = (KPP_REAL)1.000; |
---|
| 2320 | |
---|
| 2321 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 2322 | ros_Gamma[0] = (KPP_REAL) 0.2500000000000000; |
---|
| 2323 | ros_Gamma[1] = (KPP_REAL)-0.1043000000000000; |
---|
| 2324 | ros_Gamma[2] = (KPP_REAL) 0.1035000000000000; |
---|
| 2325 | ros_Gamma[3] = (KPP_REAL)-0.3620000000000023e-01; |
---|
| 2326 | ros_Gamma[4] = (KPP_REAL) 0.0; |
---|
| 2327 | ros_Gamma[5] = (KPP_REAL) 0.0; |
---|
| 2328 | |
---|
| 2329 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 2330 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 2331 | A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc. |
---|
| 2332 | The general mapping formula is: A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ] |
---|
| 2333 | C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */ |
---|
| 2334 | |
---|
| 2335 | ros_A[0] = (KPP_REAL) 0.1544000000000000e01; |
---|
| 2336 | ros_A[1] = (KPP_REAL) 0.9466785280815826; |
---|
| 2337 | ros_A[2] = (KPP_REAL) 0.2557011698983284; |
---|
| 2338 | ros_A[3] = (KPP_REAL) 0.3314825187068521e01; |
---|
| 2339 | ros_A[4] = (KPP_REAL) 0.2896124015972201e01; |
---|
| 2340 | ros_A[5] = (KPP_REAL) 0.9986419139977817; |
---|
| 2341 | ros_A[6] = (KPP_REAL) 0.1221224509226641e01; |
---|
| 2342 | ros_A[7] = (KPP_REAL) 0.6019134481288629e01; |
---|
| 2343 | ros_A[8] = (KPP_REAL) 0.1253708332932087e02; |
---|
| 2344 | ros_A[9] = (KPP_REAL)-0.6878860361058950; |
---|
| 2345 | ros_A[10] = ros_A[6]; |
---|
| 2346 | ros_A[11] = ros_A[7]; |
---|
| 2347 | ros_A[12] = ros_A[8]; |
---|
| 2348 | ros_A[13] = ros_A[9]; |
---|
| 2349 | ros_A[14] = (KPP_REAL)1.0; |
---|
| 2350 | |
---|
| 2351 | ros_C[0] = (KPP_REAL)-0.5668800000000000e01; |
---|
| 2352 | ros_C[1] = (KPP_REAL)-0.2430093356833875e01; |
---|
| 2353 | ros_C[2] = (KPP_REAL)-0.2063599157091915; |
---|
| 2354 | ros_C[3] = (KPP_REAL)-0.1073529058151375; |
---|
| 2355 | ros_C[4] = (KPP_REAL)-0.9594562251023355e01; |
---|
| 2356 | ros_C[5] = (KPP_REAL)-0.2047028614809616e02; |
---|
| 2357 | ros_C[6] = (KPP_REAL) 0.7496443313967647e01; |
---|
| 2358 | ros_C[7] = (KPP_REAL)-0.1024680431464352e02; |
---|
| 2359 | ros_C[8] = (KPP_REAL)-0.3399990352819905e02; |
---|
| 2360 | ros_C[9] = (KPP_REAL) 0.1170890893206160e02; |
---|
| 2361 | ros_C[10] = (KPP_REAL) 0.8083246795921522e01; |
---|
| 2362 | ros_C[11] = (KPP_REAL)-0.7981132988064893e01; |
---|
| 2363 | ros_C[12] = (KPP_REAL)-0.3152159432874371e02; |
---|
| 2364 | ros_C[13] = (KPP_REAL) 0.1631930543123136e02; |
---|
| 2365 | ros_C[14] = (KPP_REAL)-0.6058818238834054e01; |
---|
| 2366 | |
---|
| 2367 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 2368 | ros_M[0] = ros_A[6]; |
---|
| 2369 | ros_M[1] = ros_A[7]; |
---|
| 2370 | ros_M[2] = ros_A[8]; |
---|
| 2371 | ros_M[3] = ros_A[9]; |
---|
| 2372 | ros_M[4] = (KPP_REAL)1.0; |
---|
| 2373 | ros_M[5] = (KPP_REAL)1.0; |
---|
| 2374 | |
---|
| 2375 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 2376 | ros_E[0] = (KPP_REAL)0.0; |
---|
| 2377 | ros_E[1] = (KPP_REAL)0.0; |
---|
| 2378 | ros_E[2] = (KPP_REAL)0.0; |
---|
| 2379 | ros_E[3] = (KPP_REAL)0.0; |
---|
| 2380 | ros_E[4] = (KPP_REAL)0.0; |
---|
| 2381 | ros_E[5] = (KPP_REAL)1.0; |
---|
| 2382 | |
---|
| 2383 | /*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE) |
---|
| 2384 | or does it re-use the function evaluation from stage i-1 |
---|
| 2385 | (ros_NewF[i]=FALSE) */ |
---|
| 2386 | ros_NewF[0] = TRUE; |
---|
| 2387 | ros_NewF[1] = TRUE; |
---|
| 2388 | ros_NewF[2] = TRUE; |
---|
| 2389 | ros_NewF[3] = TRUE; |
---|
| 2390 | ros_NewF[4] = TRUE; |
---|
| 2391 | ros_NewF[5] = TRUE; |
---|
| 2392 | |
---|
| 2393 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 2394 | main and the embedded scheme orders plus 1 */ |
---|
| 2395 | ros_ELO = (KPP_REAL)4.0; |
---|
| 2396 | |
---|
| 2397 | } /* End of Rodas4 */ |
---|
| 2398 | |
---|
| 2399 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2400 | void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ) { |
---|
| 2401 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2402 | Template for the ODE function call. |
---|
| 2403 | Updates the rate coefficients (and possibly the fixed species) at each call |
---|
| 2404 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2405 | |
---|
| 2406 | /*~~~> Local variables */ |
---|
| 2407 | KPP_REAL Told; |
---|
| 2408 | |
---|
| 2409 | Told = TIME; |
---|
| 2410 | TIME = T; |
---|
| 2411 | Update_SUN(); |
---|
| 2412 | Update_RCONST(); |
---|
| 2413 | Fun( Y, FIX, RCONST, Ydot ); |
---|
| 2414 | TIME = Told; |
---|
| 2415 | |
---|
| 2416 | } /* End of FunTemplate */ |
---|
| 2417 | |
---|
| 2418 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2419 | void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ) { |
---|
| 2420 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2421 | Template for the ODE Jacobian call. |
---|
| 2422 | Updates the rate coefficients (and possibly the fixed species) at each call |
---|
| 2423 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2424 | |
---|
| 2425 | /*~~~> Local variables */ |
---|
| 2426 | KPP_REAL Told; |
---|
| 2427 | #ifdef FULL_ALGEBRA |
---|
| 2428 | KPP_REAL JV[LU_NONZERO]; |
---|
| 2429 | int i, j; |
---|
| 2430 | #endif |
---|
| 2431 | |
---|
| 2432 | Told = TIME; |
---|
| 2433 | TIME = T; |
---|
| 2434 | Update_SUN(); |
---|
| 2435 | Update_RCONST(); |
---|
| 2436 | #ifdef FULL_ALGEBRA |
---|
| 2437 | Jac_SP(Y, FIX, RCONST, JV); |
---|
| 2438 | for(j=0; j<NVAR; j++) { |
---|
| 2439 | for(i=0; i<NVAR; i++) |
---|
| 2440 | Jcb[i][j] = (KPP_REAL)0.0; |
---|
| 2441 | } |
---|
| 2442 | for(i=0; i<LU_NONZERO; i++) |
---|
| 2443 | Jcb[LU_ICOL[i]][LU_IROW[i]] = JV[i]; |
---|
| 2444 | #else |
---|
| 2445 | Jac_SP( Y, FIX, RCONST, Jcb ); |
---|
| 2446 | #endif |
---|
| 2447 | TIME = Told; |
---|
| 2448 | } /* End of JacTemplate */ |
---|
| 2449 | |
---|
| 2450 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2451 | void HessTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Hes[] ) { |
---|
| 2452 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2453 | Template for the ODE Hessian call. |
---|
| 2454 | Updates the rate coefficients (and possibly the fixed species) at each call |
---|
| 2455 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2456 | |
---|
| 2457 | /*~~~> Local variables */ |
---|
| 2458 | KPP_REAL Told; |
---|
| 2459 | |
---|
| 2460 | Told = TIME; |
---|
| 2461 | TIME = T; |
---|
| 2462 | Update_SUN(); |
---|
| 2463 | Update_RCONST(); |
---|
| 2464 | Hessian( Y, FIX, RCONST, Hes ); |
---|
| 2465 | TIME = Told; |
---|
| 2466 | |
---|
| 2467 | } /* End of HessTemplate */ |
---|
| 2468 | |
---|
| 2469 | /* End of INTEGRATE function */ |
---|
| 2470 | /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|