[2696] | 1 | /****************************************************************************** |
---|
| 2 | |
---|
| 3 | KPP - The Kinetic PreProcessor |
---|
| 4 | Builds simulation code for chemical kinetic systems |
---|
| 5 | |
---|
[4843] | 6 | Copyright (C) -2021 996 Valeriu Damian and Adrian Sandu |
---|
| 7 | Copyright (C) -2021 005 Adrian Sandu |
---|
[2696] | 8 | |
---|
| 9 | KPP is free software; you can redistribute it and/or modify it under the |
---|
| 10 | terms of the GNU General Public License as published by the Free Software |
---|
| 11 | Foundation (http://www.gnu.org/copyleft/gpl.html); either version 2 of the |
---|
| 12 | License, or (at your option) any later version. |
---|
| 13 | |
---|
| 14 | KPP is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 15 | WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
---|
| 16 | FOR A PARTICULAR PURPOSE. See the GNU General Public License for more |
---|
| 17 | details. |
---|
| 18 | |
---|
| 19 | You should have received a copy of the GNU General Public License along |
---|
| 20 | with this program; if not, consult http://www.gnu.org/copyleft/gpl.html or |
---|
| 21 | write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, |
---|
| 22 | Boston, MA 02111-1307, USA. |
---|
| 23 | |
---|
| 24 | Adrian Sandu |
---|
| 25 | Computer Science Department |
---|
| 26 | Virginia Polytechnic Institute and State University |
---|
| 27 | Blacksburg, VA 24060 |
---|
| 28 | E-mail: sandu@cs.vt.edu |
---|
| 29 | |
---|
| 30 | ******************************************************************************/ |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | #include <string.h> |
---|
| 34 | #include <math.h> |
---|
| 35 | #include "gdata.h" |
---|
| 36 | #include "code.h" |
---|
| 37 | #include "scan.h" |
---|
| 38 | |
---|
| 39 | #define MAX_MONITOR 8 |
---|
| 40 | enum strutypes { PLAIN, LU }; |
---|
| 41 | |
---|
| 42 | int **structB; |
---|
| 43 | int **structJ; |
---|
| 44 | int **LUstructJ; |
---|
| 45 | |
---|
| 46 | ICODE InlineCode[ INLINE_OPT ]; |
---|
| 47 | |
---|
| 48 | int NSPEC, NVAR, NVARACT, NFIX, NREACT; |
---|
| 49 | int NVARST, NFIXST; |
---|
| 50 | /* int PI; */ |
---|
| 51 | int C_DEFAULT, C; |
---|
| 52 | int DC; |
---|
| 53 | int ARP, JVRP, NJVRP, CROW_JVRP, IROW_JVRP, ICOL_JVRP; |
---|
| 54 | int V, F, VAR, FIX; |
---|
| 55 | int RCONST, RCT; |
---|
| 56 | int Vdot, P_VAR, D_VAR; |
---|
| 57 | int KR, A, BV, BR, IV; |
---|
| 58 | int JV, UV, JUV, JTUV, JVS; |
---|
| 59 | int JR, UR, JUR, JRS; |
---|
| 60 | int U1, U2, HU, HTU; |
---|
| 61 | int X, XX, NTMPB; |
---|
| 62 | int D2A, NTMPD2A, NHESS, HESS, IHESS_I, IHESS_J, IHESS_K; |
---|
| 63 | int DDMTYPE; |
---|
| 64 | int STOICM, NSTOICM, IROW_STOICM, ICOL_STOICM, CCOL_STOICM, CNEQN; |
---|
| 65 | int IROW, ICOL, CROW, DIAG; |
---|
| 66 | int LU_IROW, LU_ICOL, LU_CROW, LU_DIAG, CNVAR; |
---|
| 67 | int LOOKAT, NLOOKAT, MONITOR, NMONITOR; |
---|
| 68 | int NMASS, SMASS; |
---|
| 69 | int SPC_NAMES, EQN_NAMES; |
---|
| 70 | int EQN_TAGS; |
---|
| 71 | int NONZERO, LU_NONZERO; |
---|
| 72 | int TIME, SUN, TEMP; |
---|
| 73 | int RTOLS, TSTART, TEND, DT; |
---|
| 74 | int ATOL, RTOL, STEPMIN, STEPMAX, CFACTOR; |
---|
| 75 | int V_USER, CL; |
---|
| 76 | int NMLCV, NMLCF, SCT, PROPENSITY, VOLUME, IRCT; |
---|
| 77 | |
---|
| 78 | int Jac_NZ, LU_Jac_NZ, nzr; |
---|
| 79 | |
---|
| 80 | NODE *sum, *prod; |
---|
| 81 | int real; |
---|
| 82 | int nlookat; |
---|
| 83 | int nmoni; |
---|
| 84 | int ntrans; |
---|
| 85 | int nmass; |
---|
| 86 | char * CommonName; |
---|
| 87 | |
---|
| 88 | int Hess_NZ, *iHess_i, *iHess_j, *iHess_k; |
---|
| 89 | int nnz_stoicm; |
---|
| 90 | |
---|
| 91 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 92 | char * ascii(int x) |
---|
| 93 | { |
---|
| 94 | static char s[40]; |
---|
| 95 | |
---|
| 96 | sprintf(s, "%d", x); |
---|
| 97 | return s; |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 101 | char * ascid(double x) |
---|
| 102 | { |
---|
| 103 | static char s[40]; |
---|
| 104 | |
---|
| 105 | sprintf(s, "%12.6e", x); |
---|
| 106 | /* if (useDouble && ( (useLang==F77_LANG)||(useLang==F90_LANG) ) ) { */ |
---|
| 107 | if (useDouble && (useLang==F77_LANG)) |
---|
| 108 | s[strlen(s)-4] = 'd'; |
---|
| 109 | if (useDouble && (useLang==F90_LANG)) |
---|
| 110 | sprintf(s, "%s_dp",s); |
---|
| 111 | return s; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 115 | NODE * RConst( int n ) |
---|
| 116 | { |
---|
| 117 | switch( kr[n].type ) { |
---|
| 118 | case NUMBER: return Const( kr[n].val.f ); |
---|
| 119 | case PHOTO: |
---|
| 120 | case EXPRESION: return Elm( RCT, n ); |
---|
| 121 | } |
---|
| 122 | return 0; |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 128 | void InitGen() |
---|
| 129 | { |
---|
| 130 | int i,j; |
---|
| 131 | |
---|
| 132 | NSPEC = DefConst( "NSPEC", INT, "Number of chemical species" ); |
---|
| 133 | NVAR = DefConst( "NVAR", INT, "Number of Variable species" ); |
---|
| 134 | NVARACT = DefConst( "NVARACT", INT, "Number of Active species" ); |
---|
| 135 | NFIX = DefConst( "NFIX", INT, "Number of Fixed species" ); |
---|
| 136 | NREACT = DefConst( "NREACT", INT, "Number of reactions" ); |
---|
| 137 | NVARST = DefConst( "NVARST", INT, "Starting of variables in conc. vect." ); |
---|
| 138 | NFIXST = DefConst( "NFIXST", INT, "Starting of fixed in conc. vect." ); |
---|
| 139 | NONZERO = DefConst( "NONZERO", INT, "Number of nonzero entries in Jacobian" ); |
---|
| 140 | LU_NONZERO = DefConst( "LU_NONZERO", INT, "Number of nonzero entries in LU factoriz. of Jacobian" ); |
---|
| 141 | CNVAR = DefConst( "CNVAR", INT, "(NVAR+1) Number of elements in compressed row format" ); |
---|
| 142 | CNEQN = DefConst( "CNEQN", INT, "(NREACT+1) Number stoicm elements in compressed col format" ); |
---|
| 143 | |
---|
| 144 | /* PI = DefConst( "PI", real, "Value of pi" ); */ |
---|
| 145 | |
---|
| 146 | VAR = DefvElm( "VAR", real, -NVAR, "Concentrations of variable species (global)" ); |
---|
| 147 | FIX = DefvElm( "FIX", real, -NFIX, "Concentrations of fixed species (global)" ); |
---|
| 148 | |
---|
| 149 | V = DefvElm( "V", real, -NVAR, "Concentrations of variable species (local)" ); |
---|
| 150 | F = DefvElm( "F", real, -NFIX, "Concentrations of fixed species (local)" ); |
---|
| 151 | |
---|
| 152 | V_USER = DefvElm( "V_USER", real, -NVAR, "Concentration of variable species in USER's order" ); |
---|
| 153 | |
---|
| 154 | RCONST = DefvElm( "RCONST", real, -NREACT, "Rate constants (global)" ); |
---|
| 155 | RCT = DefvElm( "RCT", real, -NREACT, "Rate constants (local)" ); |
---|
| 156 | |
---|
| 157 | Vdot = DefvElm( "Vdot", real, -NVAR, "Time derivative of variable species concentrations" ); |
---|
| 158 | P_VAR = DefvElm( "P_VAR", real, -NVAR, "Production term" ); |
---|
| 159 | D_VAR = DefvElm( "D_VAR", real, -NVAR, "Destruction term" ); |
---|
| 160 | |
---|
| 161 | |
---|
| 162 | JVS = DefvElm( "JVS", real, -LU_NONZERO, "sparse Jacobian of variables" ); |
---|
| 163 | |
---|
| 164 | JV = DefmElm( "JV", real, -NVAR, -NVAR, "full Jacobian of variables" ); |
---|
| 165 | |
---|
| 166 | UV = DefvElm( "UV", real, -NVAR, "User vector for variables" ); |
---|
| 167 | JUV = DefvElm( "JUV", real, -NVAR, "Jacobian times user vector" ); |
---|
| 168 | JTUV = DefvElm( "JTUV",real, -NVAR, "Jacobian transposed times user vector" ); |
---|
| 169 | |
---|
| 170 | X = DefvElm( "X", real, -NVAR, "Vector for variables" ); |
---|
| 171 | XX = DefvElm( "XX", real, -NVAR, "Vector for output variables" ); |
---|
| 172 | |
---|
| 173 | TIME = DefElm( "TIME", real, "Current integration time"); |
---|
| 174 | SUN = DefElm( "SUN", real, "Sunlight intensity between [0,1]"); |
---|
| 175 | TEMP = DefElm( "TEMP", real, "Temperature"); |
---|
| 176 | |
---|
| 177 | RTOLS = DefElm( "RTOLS", real, "(scalar) Relative tolerance"); |
---|
| 178 | TSTART = DefElm( "TSTART", real, "Integration start time"); |
---|
| 179 | TEND = DefElm( "TEND", real, "Integration end time"); |
---|
| 180 | DT = DefElm( "DT", real, "Integration step"); |
---|
| 181 | |
---|
| 182 | A = DefvElm( "A", real, -NREACT, "Rate for each equation" ); |
---|
| 183 | |
---|
| 184 | ARP = DefvElm( "ARP", real, -NREACT, "Reactant product in each equation" ); |
---|
| 185 | NJVRP = DefConst( "NJVRP", INT, "Length of sparse Jacobian JVRP" ); |
---|
| 186 | JVRP = DefvElm( "JVRP", real, -NJVRP, "d ARP(1:NREACT)/d VAR (1:NVAR)" ); |
---|
| 187 | CROW_JVRP= DefvElm( "CROW_JVRP", INT, -CNEQN, "Beginning of rows in JVRP" ); |
---|
| 188 | ICOL_JVRP= DefvElm( "ICOL_JVRP", INT, -NJVRP, "Column indices in JVRP" ); |
---|
| 189 | IROW_JVRP= DefvElm( "IROW_JVRP", INT, -NJVRP, "Row indices in JVRP" ); |
---|
| 190 | |
---|
| 191 | NTMPB = DefConst( "NTMPB", INT, "Length of Temporary Array B" ); |
---|
| 192 | BV = DefvElm( "B", real, -NTMPB, "Temporary array" ); |
---|
| 193 | |
---|
| 194 | NSTOICM = DefConst("NSTOICM", INT, "Length of Sparse Stoichiometric Matrix" ); |
---|
| 195 | STOICM = DefvElm( "STOICM", real, -NSTOICM, "Stoichiometric Matrix in compressed column format" ); |
---|
| 196 | IROW_STOICM = DefvElm( "IROW_STOICM", INT, -NSTOICM, "Row indices in STOICM" ); |
---|
| 197 | ICOL_STOICM = DefvElm( "ICOL_STOICM", INT, -NSTOICM, "Column indices in STOICM" ); |
---|
| 198 | CCOL_STOICM = DefvElm( "CCOL_STOICM", INT, -CNEQN, "Beginning of columns in STOICM" ); |
---|
| 199 | |
---|
| 200 | DDMTYPE = DefElm( "DDMTYPE", INT, "DDM sensitivity w.r.t.: 0=init.val., 1=params" ); |
---|
| 201 | |
---|
| 202 | NTMPD2A= DefConst( "NTMPD2A", INT, "Length of Temporary Array D2A" ); |
---|
| 203 | D2A = DefvElm( "D2A", real, -NTMPD2A, "Second derivatives of equation rates" ); |
---|
| 204 | NHESS = DefConst( "NHESS", INT, "Length of Sparse Hessian" ); |
---|
| 205 | HESS = DefvElm( "HESS", real, -NHESS, "Hessian of Var (i.e. the 3-tensor d Jac / d Var)" ); |
---|
| 206 | IHESS_I = DefvElm( "IHESS_I", INT, -NHESS, "Index i of Hessian element d^2 f_i/dv_j.dv_k" ); |
---|
| 207 | IHESS_J = DefvElm( "IHESS_J", INT, -NHESS, "Index j of Hessian element d^2 f_i/dv_j.dv_k" ); |
---|
| 208 | IHESS_K = DefvElm( "IHESS_K", INT, -NHESS, "Index k of Hessian element d^2 f_i/dv_j.dv_k" ); |
---|
| 209 | U1 = DefvElm( "U1", real, -NVAR, "User vector" ); |
---|
| 210 | U2 = DefvElm( "U2", real, -NVAR, "User vector" ); |
---|
| 211 | HU = DefvElm( "HU", real, -NVAR, "Hessian times user vectors: (Hess x U2) * U1 = [d (Jac*U1)/d Var] * U2" ); |
---|
| 212 | HTU = DefvElm( "HTU", real, -NVAR, "Transposed Hessian times user vectors: (Hess x U2)^T * U1 = [d (Jac^T*U1)/d Var] * U2 " ); |
---|
| 213 | |
---|
| 214 | KR = DefeElm( "KR", 0 ); |
---|
| 215 | |
---|
| 216 | IROW = DefvElm( "IROW", INT, -NONZERO, "Row indexes of the Jacobian of variables" ); |
---|
| 217 | ICOL = DefvElm( "ICOL", INT, -NONZERO, "Column indexes of the Jacobian of variables" ); |
---|
| 218 | CROW = DefvElm( "CROW", INT, -CNVAR, "Compressed row indexes of the Jacobian of variables" ); |
---|
| 219 | DIAG = DefvElm( "DIAG", INT, -CNVAR, "Diagonal indexes of the Jacobian of variables" ); |
---|
| 220 | LU_IROW = DefvElm( "LU_IROW", INT, -LU_NONZERO, "Row indexes of the LU Jacobian of variables" ); |
---|
| 221 | LU_ICOL = DefvElm( "LU_ICOL", INT, -LU_NONZERO, "Column indexes of the LU Jacobian of variables" ); |
---|
| 222 | LU_CROW = DefvElm( "LU_CROW", INT, -CNVAR, "Compressed row indexes of the LU Jacobian of variables" ); |
---|
| 223 | LU_DIAG = DefvElm( "LU_DIAG", INT, -CNVAR, "Diagonal indexes of the LU Jacobian of variables" ); |
---|
| 224 | |
---|
| 225 | IV = DefeElm( "IV", 0 ); |
---|
| 226 | |
---|
| 227 | C_DEFAULT = DefvElm( "C_DEFAULT", real, -NSPEC, "Default concentration for all species" ); |
---|
| 228 | C = DefvElm( "C", real, -NSPEC, "Concentration of all species" ); |
---|
| 229 | CL = DefvElm( "CL", real, -NSPEC, "Concentration of all species (local)" ); |
---|
| 230 | DC = DefvElm( "DC", real, -NSPEC, "Fluxes of all species" ); |
---|
| 231 | ATOL = DefvElm( "ATOL", real, -NVAR, "Absolute tolerance" ); |
---|
| 232 | RTOL = DefvElm( "RTOL", real, -NVAR, "Relative tolerance" ); |
---|
| 233 | |
---|
| 234 | STEPMIN = DefElm( "STEPMIN", real, "Lower bound for integration step"); |
---|
| 235 | STEPMAX = DefElm( "STEPMAX", real, "Upper bound for integration step"); |
---|
| 236 | |
---|
| 237 | NLOOKAT = DefConst( "NLOOKAT", INT, "Number of species to look at" ); |
---|
| 238 | LOOKAT = DefvElm( "LOOKAT", INT, -NLOOKAT, "Indexes of species to look at" ); |
---|
| 239 | |
---|
| 240 | NMONITOR = DefConst( "NMONITOR", INT, "Number of species to monitor" ); |
---|
| 241 | MONITOR = DefvElm( "MONITOR", INT, -NMONITOR, "Indexes of species to monitor" ); |
---|
| 242 | |
---|
| 243 | NMASS = DefConst( "NMASS", INT, "Number of atoms to check mass balance" ); |
---|
| 244 | SMASS = DefvElm( "SMASS", STRING, -NMASS, "Names of atoms for mass balance" ); |
---|
| 245 | |
---|
| 246 | EQN_TAGS = DefvElm( "EQN_TAGS", STRING, -NREACT, "Equation tags" ); |
---|
| 247 | EQN_NAMES = DefvElm( "EQN_NAMES", DOUBLESTRING, -NREACT, "Equation names" ); |
---|
| 248 | SPC_NAMES = DefvElm( "SPC_NAMES", STRING, -NSPEC, "Names of chemical species" ); |
---|
| 249 | |
---|
| 250 | CFACTOR = DefElm( "CFACTOR", real, "Conversion factor for concentration units"); |
---|
| 251 | |
---|
| 252 | /* Elements of Stochastic simulation*/ |
---|
| 253 | NMLCV = DefvElm( "NmlcV", INT, -NVAR, "No. molecules of variable species" ); |
---|
| 254 | NMLCF = DefvElm( "NmlcF", INT, -NFIX, "No. molecules of fixed species" ); |
---|
| 255 | SCT = DefvElm( "SCT", real, -NREACT, "Stochastic rate constants" ); |
---|
| 256 | PROPENSITY = DefvElm( "Prop", real, -NREACT, "Propensity vector" ); |
---|
| 257 | VOLUME = DefElm( "Volume", real, "Volume of the reaction container" ); |
---|
| 258 | IRCT = DefElm( "IRCT", INT, "Index of chemical reaction" ); |
---|
| 259 | |
---|
| 260 | for ( i=0; i<EqnNr; i++ ) |
---|
| 261 | for ( j=0; j<SpcNr; j++ ) |
---|
| 262 | structB[i][j] = ( Stoich_Left[j][i] != 0 ) ? 1 : 0; |
---|
| 263 | |
---|
| 264 | /* Constant values are useful to declare vectors of this size */ |
---|
| 265 | if (useDeclareValues) { |
---|
| 266 | varTable[ NSPEC ] -> value = max(SpcNr,1); |
---|
| 267 | varTable[ NVAR ] -> value = max(VarNr,1); |
---|
| 268 | varTable[ NVARACT ] -> value = max(VarActiveNr,1); |
---|
| 269 | varTable[ NFIX ] -> value = max(FixNr,1); |
---|
| 270 | varTable[ NREACT ] -> value = max(EqnNr,1); |
---|
| 271 | varTable[ NVARST ] -> value = Index(0); |
---|
| 272 | varTable[ NFIXST ] -> value = Index(VarNr); |
---|
| 273 | } |
---|
| 274 | } |
---|
| 275 | |
---|
| 276 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 277 | int NonZero( int stru, int start, int end, |
---|
| 278 | int *row, int *col, int *crow, int *diag ) |
---|
| 279 | { |
---|
| 280 | int nElm; |
---|
| 281 | int i,j; |
---|
| 282 | |
---|
| 283 | nElm = 0; |
---|
| 284 | for (i = 0; i < end-start; i++) { |
---|
| 285 | crow[i] = Index(nElm); |
---|
| 286 | for (j = 0; j < end-start; j++) { |
---|
| 287 | if( (i == j) || ( (stru) ? LUstructJ[i+start][j+start] |
---|
| 288 | : structJ[i+start][j+start] ) ) { |
---|
| 289 | row[nElm] = Index(i); |
---|
| 290 | col[nElm] = Index(j); |
---|
| 291 | nElm++; |
---|
| 292 | } |
---|
| 293 | if( i == j ) { |
---|
| 294 | diag[i] = Index(nElm-1); |
---|
| 295 | } |
---|
| 296 | } |
---|
| 297 | } |
---|
| 298 | crow[i] = Index(nElm); |
---|
| 299 | diag[i] = Index(nElm); |
---|
| 300 | return nElm; |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 305 | void GenerateGData() |
---|
| 306 | { |
---|
| 307 | int i,j,k; |
---|
| 308 | int *crow; |
---|
| 309 | int *diag; |
---|
| 310 | int nElm; |
---|
| 311 | int *lookat; |
---|
| 312 | int *moni; |
---|
| 313 | char *snames[MAX_SPECIES]; |
---|
| 314 | int *trans; |
---|
| 315 | char *strans[MAX_SPECIES]; |
---|
| 316 | char *smass[MAX_ATOMS]; |
---|
| 317 | char *EQN_NAMES[MAX_EQN]; |
---|
| 318 | char *EQN_TAGS[MAX_EQN]; |
---|
| 319 | char *bufeqn, *p; |
---|
| 320 | int dim; |
---|
| 321 | |
---|
| 322 | if ( (useLang != C_LANG)&&(useLang != MATLAB_LANG) ) return; |
---|
| 323 | |
---|
| 324 | UseFile( driverFile ); |
---|
| 325 | |
---|
| 326 | NewLines(1); |
---|
| 327 | |
---|
| 328 | GlobalDeclare( C ); |
---|
| 329 | C_Inline("%s * %s = & %s[%d];", C_types[real], |
---|
| 330 | varTable[VAR]->name, varTable[C]->name, 0 ); |
---|
| 331 | C_Inline("%s * %s = & %s[%d];", C_types[real], |
---|
| 332 | varTable[FIX]->name, varTable[C]->name, VarNr ); |
---|
| 333 | |
---|
| 334 | |
---|
| 335 | GlobalDeclare( RCONST ); |
---|
| 336 | GlobalDeclare( TIME ); |
---|
| 337 | GlobalDeclare( SUN ); |
---|
| 338 | GlobalDeclare( TEMP ); |
---|
| 339 | GlobalDeclare( RTOLS ); |
---|
| 340 | GlobalDeclare( TSTART ); |
---|
| 341 | GlobalDeclare( TEND ); |
---|
| 342 | GlobalDeclare( DT ); |
---|
| 343 | GlobalDeclare( ATOL ); |
---|
| 344 | GlobalDeclare( RTOL ); |
---|
| 345 | GlobalDeclare( STEPMIN ); |
---|
| 346 | GlobalDeclare( STEPMAX ); |
---|
| 347 | GlobalDeclare( CFACTOR ); |
---|
| 348 | if (useStochastic) |
---|
| 349 | GlobalDeclare( VOLUME ); |
---|
| 350 | |
---|
| 351 | MATLAB_Inline(" %s_Parameters;",rootFileName); |
---|
| 352 | MATLAB_Inline(" %s_Global_defs;",rootFileName); |
---|
| 353 | MATLAB_Inline(" %s_Sparse;",rootFileName); |
---|
| 354 | MATLAB_Inline(" %s_Monitor;",rootFileName); |
---|
| 355 | if (useJacSparse ) |
---|
| 356 | MATLAB_Inline(" %s_JacobianSP;",rootFileName); |
---|
| 357 | if (useHessian ) |
---|
| 358 | MATLAB_Inline(" %s_HessianSP;",rootFileName); |
---|
| 359 | if (useStoicmat ) |
---|
| 360 | MATLAB_Inline(" %s_StoichiomSP;",rootFileName); |
---|
| 361 | |
---|
| 362 | NewLines(1); |
---|
| 363 | |
---|
| 364 | } |
---|
| 365 | |
---|
| 366 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 367 | void GenerateMonitorData() |
---|
| 368 | { |
---|
| 369 | int i,j,k; |
---|
| 370 | int *crow; |
---|
| 371 | int *diag; |
---|
| 372 | int nElm; |
---|
| 373 | int *lookat; |
---|
| 374 | int *moni; |
---|
| 375 | char *snames[MAX_SPECIES]; |
---|
| 376 | int *trans; |
---|
| 377 | char *strans[MAX_SPECIES]; |
---|
| 378 | char *smass[MAX_ATOMS]; |
---|
| 379 | char *seqn[MAX_EQN]; |
---|
| 380 | char *bufeqn, *p; |
---|
| 381 | int dim; |
---|
| 382 | |
---|
| 383 | |
---|
| 384 | /* Allocate local data structures */ |
---|
| 385 | dim = SpcNr+2; |
---|
| 386 | crow = AllocIntegerVector( dim, "crow in GenerateMonitorData"); |
---|
| 387 | diag = AllocIntegerVector( dim, "diag in GenerateMonitorData"); |
---|
| 388 | lookat = AllocIntegerVector( dim, "lookat in GenerateMonitorData"); |
---|
| 389 | moni = AllocIntegerVector( dim, "moni in GenerateMonitorData"); |
---|
| 390 | trans = AllocIntegerVector( dim, "trans in GenerateMonitorData"); |
---|
| 391 | |
---|
| 392 | UseFile( monitorFile ); |
---|
| 393 | |
---|
| 394 | F77_Inline("%6sBLOCK DATA MONITOR_DATA\n", " " ); |
---|
| 395 | F77_Inline("%6sINCLUDE '%s_Parameters.h'", " ",rootFileName); |
---|
| 396 | F77_Inline("%6sINCLUDE '%s_Global.h'", " ",rootFileName); |
---|
| 397 | F77_Inline("%6sINTEGER i", " " ); |
---|
| 398 | |
---|
| 399 | /* InitDeclare( CFACTOR, 0, (void*)&cfactor ); */ |
---|
| 400 | |
---|
| 401 | NewLines(1); |
---|
| 402 | |
---|
| 403 | for (i = 0; i < SpcNr; i++) { |
---|
| 404 | snames[i] = SpeciesTable[Code[i]].name; |
---|
| 405 | } |
---|
| 406 | InitDeclare( SPC_NAMES, SpcNr, (void*)snames ); |
---|
| 407 | |
---|
| 408 | nlookat = 0; |
---|
| 409 | for (i = 0; i < SpcNr; i++) |
---|
| 410 | if ( SpeciesTable[Code[i]].lookat ) { |
---|
| 411 | lookat[nlookat] = Index(i); |
---|
| 412 | nlookat++; |
---|
| 413 | } |
---|
| 414 | |
---|
| 415 | if (useDeclareValues) |
---|
| 416 | varTable[ NLOOKAT ] -> value = max(nlookat,1); |
---|
| 417 | InitDeclare( LOOKAT, nlookat, (void*)lookat ); |
---|
| 418 | |
---|
| 419 | nmoni = 0; |
---|
| 420 | for (i = 0; i < SpcNr; i++) |
---|
| 421 | if ( SpeciesTable[Code[i]].moni ) { |
---|
| 422 | moni[nmoni] = Index(i); |
---|
| 423 | nmoni++; |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | if( nmoni > MAX_MONITOR ) { |
---|
| 427 | Warning( "%d species to monitorize. Too many, keeping %d.", |
---|
| 428 | nmoni, MAX_MONITOR ); |
---|
| 429 | nmoni = MAX_MONITOR; |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | if (useDeclareValues) |
---|
| 433 | varTable[ NMONITOR ] -> value = max(nmoni,1); |
---|
| 434 | InitDeclare( MONITOR, nmoni, (void*)moni ); |
---|
| 435 | |
---|
| 436 | ntrans = 0; |
---|
| 437 | for (i = 0; i < SpcNr; i++) |
---|
| 438 | if ( SpeciesTable[Code[i]].trans ) { |
---|
| 439 | trans[ntrans] = Index(i); |
---|
| 440 | strans[ntrans] = SpeciesTable[Code[i]].name; |
---|
| 441 | ntrans++; |
---|
| 442 | } |
---|
| 443 | |
---|
| 444 | nmass = 0; |
---|
| 445 | for (i = 0; i < AtomNr; i++) |
---|
| 446 | if ( AtomTable[i].masscheck ) { |
---|
| 447 | smass[nmass] = AtomTable[i].name; |
---|
| 448 | nmass++; |
---|
| 449 | } |
---|
| 450 | if (useDeclareValues) |
---|
| 451 | varTable[ NMASS ] -> value = max(nmass,1); |
---|
| 452 | InitDeclare( SMASS, nmass, (void*)smass ); |
---|
| 453 | |
---|
| 454 | if ( (bufeqn = (char*)malloc(MAX_EQNLEN*EqnNr+2))==NULL ) { |
---|
| 455 | FatalError(-30,"GenerateMonitorData: Cannot allocate bufeqn (%d chars)", |
---|
| 456 | MAX_EQNLEN*EqnNr); |
---|
| 457 | } |
---|
| 458 | |
---|
| 459 | p = bufeqn; |
---|
| 460 | for (i = 0; i < EqnNr; i++) { |
---|
| 461 | EqnString(i, p); |
---|
| 462 | seqn[i] = p; |
---|
| 463 | p += MAX_EQNLEN; |
---|
| 464 | } |
---|
| 465 | InitDeclare( EQN_NAMES, EqnNr, (void*)seqn ); |
---|
| 466 | |
---|
| 467 | free( bufeqn ); |
---|
| 468 | |
---|
| 469 | if (useEqntags==1) { |
---|
| 470 | for (i = 0; i < EqnNr; i++) { |
---|
| 471 | seqn[i] = kr[i].label; |
---|
| 472 | } |
---|
| 473 | InitDeclare( EQN_TAGS, EqnNr, (void*)seqn ); |
---|
| 474 | } |
---|
| 475 | |
---|
| 476 | NewLines(1); |
---|
| 477 | WriteComment("INLINED global variables"); |
---|
| 478 | |
---|
| 479 | switch( useLang ) { |
---|
| 480 | case C_LANG: bprintf( InlineCode[ C_DATA ].code ); |
---|
| 481 | break; |
---|
| 482 | case F77_LANG: bprintf( InlineCode[ F77_DATA ].code ); |
---|
| 483 | break; |
---|
| 484 | case F90_LANG: bprintf( InlineCode[ F90_DATA ].code ); |
---|
| 485 | break; |
---|
| 486 | case MATLAB_LANG: bprintf( InlineCode[ MATLAB_DATA ].code ); |
---|
| 487 | break; |
---|
| 488 | } |
---|
| 489 | FlushBuf(); |
---|
| 490 | |
---|
| 491 | NewLines(1); |
---|
| 492 | WriteComment("End INLINED global variables"); |
---|
| 493 | NewLines(1); |
---|
| 494 | |
---|
| 495 | F77_Inline( "%6sEND\n\n", " " ); |
---|
| 496 | |
---|
| 497 | /* Free local data structures */ |
---|
| 498 | free(crow); free(diag); free(lookat); free(moni); free(trans); |
---|
| 499 | |
---|
| 500 | } |
---|
| 501 | |
---|
| 502 | |
---|
| 503 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 504 | void GenerateJacobianSparseData() |
---|
| 505 | { |
---|
| 506 | int* irow; |
---|
| 507 | int* icol; |
---|
| 508 | int* crow; |
---|
| 509 | int* diag; |
---|
| 510 | int nElm; |
---|
| 511 | int dim; |
---|
| 512 | |
---|
| 513 | if( !useJacSparse ) return; |
---|
| 514 | |
---|
| 515 | /* Allocate local arrays */ |
---|
| 516 | dim=MAX_SPECIES; |
---|
| 517 | irow = AllocIntegerVector( dim*dim, "irow in GenerateJacobianSparseData" ); |
---|
| 518 | icol = AllocIntegerVector( dim*dim, "icol in GenerateJacobianSparseData" ); |
---|
| 519 | crow = AllocIntegerVector( dim, "crow in GenerateJacobianSparseData" ); |
---|
| 520 | diag = AllocIntegerVector( dim, "diag in GenerateJacobianSparseData" ); |
---|
| 521 | |
---|
| 522 | UseFile( sparse_jacFile ); |
---|
| 523 | |
---|
| 524 | NewLines(1); |
---|
| 525 | WriteComment("Sparse Jacobian Data"); |
---|
| 526 | NewLines(1); |
---|
| 527 | |
---|
| 528 | F77_Inline("%6sBLOCK DATA JACOBIAN_SPARSE_DATA\n", " " ); |
---|
| 529 | F77_Inline("%6sINCLUDE '%s_Sparse.h'", " ",rootFileName); |
---|
| 530 | F77_Inline("%6sINTEGER i"," "); |
---|
| 531 | /* F90_Inline(" USE %s_Sparse", rootFileName); */ |
---|
| 532 | |
---|
| 533 | |
---|
| 534 | Jac_NZ = NonZero( PLAIN, 0, VarNr, irow, icol, crow, diag ); |
---|
| 535 | LU_Jac_NZ = NonZero( LU, 0, VarNr, irow, icol, crow, diag ); |
---|
| 536 | if (useDeclareValues) { |
---|
| 537 | varTable[NONZERO] -> value = Jac_NZ; |
---|
| 538 | varTable[LU_NONZERO] -> value = LU_Jac_NZ; |
---|
| 539 | } |
---|
| 540 | |
---|
| 541 | switch (useJacobian) { |
---|
| 542 | case JAC_ROW: |
---|
| 543 | Jac_NZ = NonZero( PLAIN, 0, VarNr, irow, icol, crow, diag ); |
---|
| 544 | InitDeclare( IROW, Jac_NZ, (void*)irow ); |
---|
| 545 | InitDeclare( ICOL, Jac_NZ, (void*)icol ); |
---|
| 546 | InitDeclare( CROW, VarNr+1, (void*)crow ); |
---|
| 547 | InitDeclare( DIAG, VarNr+1, (void*)diag ); |
---|
| 548 | break; |
---|
| 549 | case JAC_LU_ROW: |
---|
| 550 | LU_Jac_NZ = NonZero( LU, 0, VarNr, irow, icol, crow, diag ); |
---|
| 551 | InitDeclare( LU_IROW, LU_Jac_NZ, (void*)irow ); |
---|
| 552 | InitDeclare( LU_ICOL, LU_Jac_NZ, (void*)icol ); |
---|
| 553 | InitDeclare( LU_CROW, VarNr+1, (void*)crow ); |
---|
| 554 | InitDeclare( LU_DIAG, VarNr+1, (void*)diag ); |
---|
| 555 | } |
---|
| 556 | NewLines(1); |
---|
| 557 | F77_Inline( "%6sEND\n\n", " " ); |
---|
| 558 | |
---|
| 559 | /* Free local arrays */ |
---|
| 560 | free(irow); free(icol); free(crow); free(diag); |
---|
| 561 | } |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | |
---|
| 565 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 566 | void GenerateJacobianSparseHeader() |
---|
| 567 | { |
---|
| 568 | UseFile( sparse_dataFile ); |
---|
| 569 | |
---|
| 570 | CommonName = "SDATA"; |
---|
| 571 | |
---|
| 572 | NewLines(1); |
---|
| 573 | WriteComment(" ----------> Sparse Jacobian Data"); |
---|
| 574 | NewLines(1); |
---|
| 575 | |
---|
| 576 | switch (useJacobian) { |
---|
| 577 | case JAC_ROW: |
---|
| 578 | ExternDeclare( IROW ); |
---|
| 579 | ExternDeclare( ICOL ); |
---|
| 580 | ExternDeclare( CROW ); |
---|
| 581 | ExternDeclare( DIAG ); |
---|
| 582 | break; |
---|
| 583 | case JAC_LU_ROW: |
---|
| 584 | ExternDeclare( LU_IROW ); |
---|
| 585 | ExternDeclare( LU_ICOL ); |
---|
| 586 | ExternDeclare( LU_CROW ); |
---|
| 587 | ExternDeclare( LU_DIAG ); |
---|
| 588 | } |
---|
| 589 | |
---|
| 590 | NewLines(1); |
---|
| 591 | } |
---|
| 592 | |
---|
| 593 | |
---|
| 594 | |
---|
| 595 | |
---|
| 596 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 597 | void GenerateFun() |
---|
| 598 | { |
---|
| 599 | int i, j, k; |
---|
| 600 | int used; |
---|
| 601 | int l, m; |
---|
| 602 | int F_VAR, FSPLIT_VAR; |
---|
| 603 | |
---|
| 604 | if( VarNr == 0 ) return; |
---|
| 605 | |
---|
| 606 | if (useLang != MATLAB_LANG) /* Matlab generates an additional file per function */ |
---|
| 607 | UseFile( functionFile ); |
---|
| 608 | |
---|
| 609 | F_VAR = DefFnc( "Fun", 4, "time derivatives of variables - Agregate form"); |
---|
| 610 | FSPLIT_VAR = DefFnc( "Fun_SPLIT", 5, "time derivatives of variables - Split form"); |
---|
| 611 | |
---|
| 612 | if( useAggregate ) |
---|
| 613 | FunctionBegin( F_VAR, V, F, RCT, Vdot ); |
---|
| 614 | else |
---|
| 615 | FunctionBegin( FSPLIT_VAR, V, F, RCT, P_VAR, D_VAR ); |
---|
| 616 | |
---|
| 617 | if ( (useLang==MATLAB_LANG)&&(!useAggregate) ) |
---|
| 618 | printf("\nWarning: in the function definition move P_VAR to output vars\n"); |
---|
| 619 | |
---|
| 620 | |
---|
| 621 | if ( useLang!=F90_LANG ) { /* A is a module variable in F90 */ |
---|
| 622 | NewLines(1); |
---|
| 623 | WriteComment("Local variables"); |
---|
| 624 | Declare( A ); |
---|
| 625 | } |
---|
| 626 | NewLines(1); |
---|
| 627 | WriteComment("Computation of equation rates"); |
---|
| 628 | |
---|
| 629 | for(j=0; j<EqnNr; j++) { |
---|
| 630 | used = 0; |
---|
| 631 | if( useAggregate ) { |
---|
| 632 | for (i = 0; i < VarNr; i++) |
---|
| 633 | if ( Stoich[i][j] != 0 ) { |
---|
| 634 | used = 1; |
---|
| 635 | break; |
---|
| 636 | } |
---|
| 637 | } else { |
---|
| 638 | for (i = 0; i < VarNr; i++) |
---|
| 639 | if ( Stoich_Right[i][j] != 0 ) { |
---|
| 640 | used = 1; |
---|
| 641 | break; |
---|
| 642 | } |
---|
| 643 | } |
---|
| 644 | |
---|
| 645 | if ( used ) { |
---|
| 646 | prod = RConst( j ); |
---|
| 647 | for (i = 0; i < VarNr; i++) |
---|
| 648 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 649 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 650 | for ( ; i < SpcNr; i++) |
---|
| 651 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 652 | prod = Mul( prod, Elm( F, i - VarNr ) ); |
---|
| 653 | Assign( Elm( A, j ), prod ); |
---|
| 654 | } |
---|
| 655 | } |
---|
| 656 | |
---|
| 657 | if( useAggregate ) { |
---|
| 658 | |
---|
| 659 | NewLines(1); |
---|
| 660 | WriteComment("Aggregate function"); |
---|
| 661 | |
---|
| 662 | for (i = 0; i < VarNr; i++) { |
---|
| 663 | sum = Const(0); |
---|
| 664 | for (j = 0; j < EqnNr; j++) |
---|
| 665 | sum = Add( sum, Mul( Const( Stoich[i][j] ), Elm( A, j ) ) ); |
---|
| 666 | Assign( Elm( Vdot, i ), sum ); |
---|
| 667 | } |
---|
| 668 | |
---|
| 669 | } else { |
---|
| 670 | |
---|
| 671 | NewLines(1); |
---|
| 672 | WriteComment("Production function"); |
---|
| 673 | |
---|
| 674 | for (i = 0; i < VarNr; i++) { |
---|
| 675 | sum = Const(0); |
---|
| 676 | for (j = 0; j < EqnNr; j++) |
---|
| 677 | sum = Add( sum, Mul( Const( Stoich_Right[i][j] ), Elm( A, j ) ) ); |
---|
| 678 | Assign( Elm( P_VAR, i ), sum ); |
---|
| 679 | } |
---|
| 680 | |
---|
| 681 | NewLines(1); |
---|
| 682 | WriteComment("Destruction function"); |
---|
| 683 | |
---|
| 684 | for (i = 0; i < VarNr; i++) { |
---|
| 685 | sum = Const(0); |
---|
| 686 | for(j=0; j<EqnNr; j++) { |
---|
| 687 | if ( Stoich_Left[i][j] == 0 ) continue; |
---|
| 688 | prod = Mul( RConst( j ), Const( Stoich_Left[i][j] ) ); |
---|
| 689 | for (l = 0; l < VarNr; l++) { |
---|
| 690 | m=(int)Stoich_Left[l][j] - (l==i); |
---|
| 691 | for (k = 1; k <= m; k++ ) |
---|
| 692 | prod = Mul( prod, Elm( V, l ) ); |
---|
| 693 | } |
---|
| 694 | for ( ; l < SpcNr; l++) |
---|
| 695 | for (k = 1; k <= (int)Stoich_Left[l][j]; k++ ) |
---|
| 696 | prod = Mul( prod, Elm( F, l - VarNr ) ); |
---|
| 697 | sum = Add( sum, prod ); |
---|
| 698 | } |
---|
| 699 | Assign( Elm( D_VAR, i ), sum ); |
---|
| 700 | } |
---|
| 701 | } |
---|
| 702 | |
---|
| 703 | if( useAggregate ) |
---|
| 704 | MATLAB_Inline("\n Vdot = Vdot(:);\n"); |
---|
| 705 | else |
---|
| 706 | MATLAB_Inline("\n P_VAR = P_VAR(:);\n D_VAR = D_VAR(:);\n"); |
---|
| 707 | |
---|
| 708 | if( useAggregate ) |
---|
| 709 | FunctionEnd( F_VAR ); |
---|
| 710 | else |
---|
| 711 | FunctionEnd( FSPLIT_VAR ); |
---|
| 712 | |
---|
| 713 | FreeVariable( F_VAR ); |
---|
| 714 | FreeVariable( FSPLIT_VAR ); |
---|
| 715 | } |
---|
| 716 | |
---|
| 717 | |
---|
| 718 | |
---|
| 719 | |
---|
| 720 | |
---|
| 721 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 722 | void GenerateStochastic() |
---|
| 723 | { |
---|
| 724 | int i, j, k, l, m, n, jnr; |
---|
| 725 | int used; |
---|
| 726 | int F_VAR; |
---|
| 727 | |
---|
| 728 | if( VarNr == 0 ) return; |
---|
| 729 | |
---|
| 730 | if (useLang != MATLAB_LANG) /* Matlab generates an additional file per function */ |
---|
| 731 | UseFile( stochasticFile ); |
---|
| 732 | |
---|
| 733 | /* ~~~~~~~> 1. PROPENSITY FUNCTION ~~~~~~~~~~~~ */ |
---|
| 734 | F_VAR = DefFnc( "Propensity", 4, "Propensity function"); |
---|
| 735 | FunctionBegin( F_VAR, NMLCV, NMLCF, SCT, PROPENSITY ); |
---|
| 736 | |
---|
| 737 | if ( (useLang==MATLAB_LANG)&&(!useAggregate) ) |
---|
| 738 | printf("\nWarning: in the function definition move P_VAR to output vars\n"); |
---|
| 739 | |
---|
| 740 | NewLines(1); |
---|
| 741 | |
---|
| 742 | for(j=0; j<EqnNr; j++) { |
---|
| 743 | used = 0; |
---|
| 744 | for (i = 0; i < VarNr; i++) |
---|
| 745 | if ( Stoich[i][j] != 0 ) { |
---|
| 746 | used = 1; |
---|
| 747 | break; |
---|
| 748 | } |
---|
| 749 | if ( used ) { |
---|
| 750 | prod = Elm( SCT, j ); |
---|
| 751 | for (i = 0; i < VarNr; i++) |
---|
| 752 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 753 | if (k==1) |
---|
| 754 | prod = Mul( prod, Elm( NMLCV, i ) ); |
---|
| 755 | else |
---|
| 756 | prod = Mul( prod, Add( Elm( NMLCV, i ), Const(-k+1) ) ); |
---|
| 757 | |
---|
| 758 | for ( ; i < SpcNr; i++) |
---|
| 759 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 760 | if (k==1) |
---|
| 761 | prod = Mul( prod, Elm( NMLCF, i - VarNr ) ); |
---|
| 762 | else |
---|
| 763 | prod = Mul( prod, Add( Elm( NMLCF, i - VarNr ), Const(-k+1) ) ); |
---|
| 764 | Assign( Elm( PROPENSITY, j ), prod ); |
---|
| 765 | } /* if used */ |
---|
| 766 | } /* for j */ |
---|
| 767 | |
---|
| 768 | MATLAB_Inline("\n Prop = Prop(:);\n"); |
---|
| 769 | FunctionEnd( F_VAR ); |
---|
| 770 | FreeVariable( F_VAR ); |
---|
| 771 | |
---|
| 772 | /* ~~~~~~~> 2. RATE CONVERSION ~~~~~~~~~~~~ */ |
---|
| 773 | F_VAR = DefFnc( "StochasticRates", 3, "Convert deterministic rates to stochastic"); |
---|
| 774 | FunctionBegin( F_VAR, RCT, VOLUME, SCT ); |
---|
| 775 | WriteComment("No. of molecules = Concentration x Volume"); |
---|
| 776 | WriteComment("For a reaction with k reactants:"); |
---|
| 777 | WriteComment(" RCT [ (molec/Volume)^(1-k) * sec^(-1) ]"); |
---|
| 778 | WriteComment(" SCT [ (molec)^(1-k) * sec^(-1) ] = RCT*Volume^(k-1)"); |
---|
| 779 | WriteComment("For p molecules of the same type: SCT = SCT/(p!)"); |
---|
| 780 | |
---|
| 781 | NewLines(1); |
---|
| 782 | |
---|
| 783 | for(j=0; j<EqnNr; j++) { |
---|
| 784 | prod = Elm( RCT, j ); |
---|
| 785 | m = 0; |
---|
| 786 | for (i = 0; i < SpcNr; i++) |
---|
| 787 | m += (int)Stoich_Left[i][j]; |
---|
| 788 | for ( i=2 ; i <= m; i++) |
---|
| 789 | prod = Mul( prod, Elm(VOLUME, 1) ); |
---|
| 790 | n = 1; |
---|
| 791 | for (i = 0; i < SpcNr; i++) |
---|
| 792 | for (k = 2; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 793 | n *= k; |
---|
| 794 | prod = Div( prod, Const( n ) ); |
---|
| 795 | Assign( Elm( SCT, j ), prod ); |
---|
| 796 | } /* for j */ |
---|
| 797 | |
---|
| 798 | MATLAB_Inline("\n SCT = SCT(:);\n"); |
---|
| 799 | FunctionEnd( F_VAR ); |
---|
| 800 | FreeVariable( F_VAR ); |
---|
| 801 | |
---|
| 802 | /* ~~~~~~~> 3. THE CHANGE IN NUMBER OF MOLECULES */ |
---|
| 803 | if (useLang == MATLAB_LANG) { |
---|
| 804 | F_VAR = DefFnc( "MoleculeChange", 3, "Change in the number of molecules"); |
---|
| 805 | FunctionBegin( F_VAR, IRCT, NMLCV, NMLCV ); |
---|
| 806 | } else { |
---|
| 807 | F_VAR = DefFnc( "MoleculeChange", 2, "Change in the number of molecules"); |
---|
| 808 | FunctionBegin( F_VAR, IRCT, NMLCV ); |
---|
| 809 | } |
---|
| 810 | |
---|
| 811 | NewLines(1); |
---|
| 812 | |
---|
| 813 | F90_Inline("\n SELECT CASE (IRCT)\n"); |
---|
| 814 | C_Inline ("\n switch (IRCT) { \n"); |
---|
| 815 | MATLAB_Inline("\n switch (IRCT) \n"); |
---|
| 816 | for(j=0; j<EqnNr; j++) { |
---|
| 817 | jnr = j+1; |
---|
| 818 | if (j==0) |
---|
| 819 | F77_Inline("\n%6sIF (IRCT.EQ.%d) THEN"," ",jnr); |
---|
| 820 | else |
---|
| 821 | F77_Inline("\n%6sELSEIF (IRCT.EQ.%d) THEN "," ",jnr); |
---|
| 822 | F90_Inline("\n CASE (%d) ",jnr); |
---|
| 823 | C_Inline("\n case %d: ",jnr); |
---|
| 824 | MATLAB_Inline("\n case %d, ",jnr); |
---|
| 825 | for (i = 0; i < VarNr; i++) { |
---|
| 826 | if ( Stoich_Left[i][j] || Stoich_Right[i][j] ) |
---|
| 827 | Assign( Elm( NMLCV, i ), Add(Elm( NMLCV, i ), |
---|
| 828 | Const(Stoich_Right[i][j]-Stoich_Left[i][j])) ); |
---|
| 829 | } /* for i */ |
---|
| 830 | C_Inline(" break;",j); |
---|
| 831 | } /* for j */ |
---|
| 832 | F77_Inline("\n%6sEND IF ! n\n"," "); |
---|
| 833 | F90_Inline("\n END SELECT\n"); |
---|
| 834 | C_Inline("\n } /* switch (IRCT) */ \n"); |
---|
| 835 | MATLAB_Inline("\n end\n"); |
---|
| 836 | |
---|
| 837 | FunctionEnd( F_VAR ); |
---|
| 838 | FreeVariable( F_VAR ); |
---|
| 839 | |
---|
| 840 | } |
---|
| 841 | |
---|
| 842 | |
---|
| 843 | |
---|
| 844 | |
---|
| 845 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 846 | void GenerateReactantProd() |
---|
| 847 | { |
---|
| 848 | int i, j, k; |
---|
| 849 | int used; |
---|
| 850 | int l, m; |
---|
| 851 | int F_STOIC; |
---|
| 852 | |
---|
| 853 | if( VarNr == 0 ) return; |
---|
| 854 | |
---|
| 855 | UseFile( stoichiomFile ); |
---|
| 856 | |
---|
| 857 | F_STOIC = DefFnc( "ReactantProd",3, "Reactant Products in each equation"); |
---|
| 858 | |
---|
| 859 | FunctionBegin( F_STOIC, V, F, ARP ); |
---|
| 860 | |
---|
| 861 | NewLines(1); |
---|
| 862 | WriteComment("Reactant Products in each equation are useful in the"); |
---|
| 863 | WriteComment(" stoichiometric formulation of mass action law"); |
---|
| 864 | |
---|
| 865 | for(j=0; j<EqnNr; j++) { |
---|
| 866 | |
---|
| 867 | prod = Const( 1 ); |
---|
| 868 | for (i = 0; i < VarNr; i++) |
---|
| 869 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 870 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 871 | for ( ; i < SpcNr; i++) |
---|
| 872 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 873 | prod = Mul( prod, Elm( F, i - VarNr ) ); |
---|
| 874 | Assign( Elm( ARP, j ), prod ); |
---|
| 875 | |
---|
| 876 | } /* for j EqnNr */ |
---|
| 877 | |
---|
| 878 | FunctionEnd( F_STOIC ); |
---|
| 879 | FreeVariable( F_STOIC ); |
---|
| 880 | } |
---|
| 881 | |
---|
| 882 | |
---|
| 883 | |
---|
| 884 | |
---|
| 885 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 886 | void GenerateJacReactantProd() |
---|
| 887 | { |
---|
| 888 | int i, j, k, l, m, JVRP_NZ, newrow; |
---|
| 889 | int used; |
---|
| 890 | int F_STOIC; |
---|
| 891 | int crow_JVRP[MAX_EQN], icol_JVRP[MAX_EQN*MAX_SPECIES]; |
---|
| 892 | int irow_JVRP[MAX_EQN*MAX_SPECIES]; |
---|
| 893 | |
---|
| 894 | if( VarNr == 0 ) return; |
---|
| 895 | |
---|
| 896 | if (useDeclareValues) { |
---|
| 897 | JVRP_NZ = -1; |
---|
| 898 | for ( i=0; i<EqnNr; i++ ) |
---|
| 899 | for ( j=0; j<VarNr; j++ ) |
---|
| 900 | if ( Stoich_Left[j][i] != 0 ) JVRP_NZ++; |
---|
| 901 | varTable[ NJVRP ] -> value = JVRP_NZ + 1; |
---|
| 902 | } |
---|
| 903 | |
---|
| 904 | UseFile( stoichiomFile ); |
---|
| 905 | |
---|
| 906 | F_STOIC = DefFnc( "JacReactantProd",3, "Jacobian of Reactant Products vector"); |
---|
| 907 | |
---|
| 908 | FunctionBegin( F_STOIC, V, F, JVRP ); |
---|
| 909 | |
---|
| 910 | |
---|
| 911 | NewLines(1); |
---|
| 912 | WriteComment("Reactant Products in each equation are useful in the"); |
---|
| 913 | WriteComment(" stoichiometric formulation of mass action law"); |
---|
| 914 | WriteComment("Below we compute the Jacobian of the Reactant Products vector"); |
---|
| 915 | WriteComment(" w.r.t. variable species: d ARP(1:NREACT) / d Var(1:NVAR)"); |
---|
| 916 | |
---|
| 917 | NewLines(1); |
---|
| 918 | |
---|
| 919 | JVRP_NZ = -1; |
---|
| 920 | for ( i=0; i<EqnNr; i++ ) { |
---|
| 921 | newrow = 0; |
---|
| 922 | crow_JVRP[i] = JVRP_NZ+1; |
---|
| 923 | for ( j=0; j<VarNr; j++ ) { |
---|
| 924 | if ( Stoich_Left[j][i] != 0 ) { |
---|
| 925 | JVRP_NZ++; |
---|
| 926 | icol_JVRP[JVRP_NZ] = j; |
---|
| 927 | irow_JVRP[JVRP_NZ] = i; |
---|
| 928 | if ( newrow == 0 ) { /* a new row begins here */ |
---|
| 929 | crow_JVRP[i] = JVRP_NZ; |
---|
| 930 | newrow = 1; |
---|
| 931 | } |
---|
| 932 | prod = Const( Stoich_Left[j][i] ) ; |
---|
| 933 | for (l = 0; l < VarNr; l++) { |
---|
| 934 | m = (int)Stoich_Left[l][i] - (l==j); |
---|
| 935 | for (k = 1; k <= m; k++ ) |
---|
| 936 | prod = Mul( prod, Elm( V, l ) ); |
---|
| 937 | } |
---|
| 938 | for ( ; l < SpcNr; l++) |
---|
| 939 | for (k = 1; k <= (int)Stoich_Left[l][i]; k++ ) |
---|
| 940 | prod = Mul( prod, Elm( F, l - VarNr ) ); |
---|
| 941 | /* Comment the B */ |
---|
| 942 | WriteComment("JVRP(%d) = dARP(%d)/dV(%d)",Index(JVRP_NZ),Index(i),Index(j)); |
---|
| 943 | Assign( Elm( JVRP, JVRP_NZ ), prod ); |
---|
| 944 | } |
---|
| 945 | } |
---|
| 946 | } |
---|
| 947 | crow_JVRP[EqnNr] = JVRP_NZ + 1; |
---|
| 948 | |
---|
| 949 | FunctionEnd( F_STOIC ); |
---|
| 950 | FreeVariable( F_STOIC ); |
---|
| 951 | |
---|
| 952 | |
---|
| 953 | UseFile( sparse_stoicmFile ); |
---|
| 954 | NewLines(1); |
---|
| 955 | WriteComment("Row-compressed sparse data for the Jacobian of reaction products JVRP"); |
---|
| 956 | F77_Inline("%6sBLOCK DATA JVRP_SPARSE_DATA\n", " " ); |
---|
| 957 | F77_Inline("%6sINCLUDE '%s_Sparse.h'", " ", rootFileName); |
---|
| 958 | F77_Inline("%6sINTEGER i", " "); |
---|
| 959 | /* F90_Inline(" USE %s_Sparse", rootFileName); */ |
---|
| 960 | if( (useLang==F77_LANG)||(useLang==F90_LANG) ) { |
---|
| 961 | for (k=0; k<JVRP_NZ+1; k++) { |
---|
| 962 | irow_JVRP[k]++; |
---|
| 963 | icol_JVRP[k]++; |
---|
| 964 | } |
---|
| 965 | for (k=0; k<EqnNr+1; k++) |
---|
| 966 | crow_JVRP[k]++; |
---|
| 967 | } |
---|
| 968 | InitDeclare( CROW_JVRP, EqnNr+1, (void*)crow_JVRP ); |
---|
| 969 | InitDeclare( ICOL_JVRP, JVRP_NZ + 1, (void*)icol_JVRP ); |
---|
| 970 | InitDeclare( IROW_JVRP, JVRP_NZ + 1, (void*)irow_JVRP ); |
---|
| 971 | NewLines(1); |
---|
| 972 | F77_Inline( "%6sEND\n\n", " " ); |
---|
| 973 | NewLines(1); |
---|
| 974 | |
---|
| 975 | |
---|
| 976 | UseFile( param_headerFile ); |
---|
| 977 | CommonName = "GDATA"; |
---|
| 978 | NewLines(1); |
---|
| 979 | DeclareConstant( NJVRP, ascii( JVRP_NZ + 1 ) ); |
---|
| 980 | |
---|
| 981 | } |
---|
| 982 | |
---|
| 983 | |
---|
| 984 | |
---|
| 985 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 986 | void GenerateJac() |
---|
| 987 | { |
---|
| 988 | int i,j,k,l,m; |
---|
| 989 | int nElm, nonzeros_B; |
---|
| 990 | int Jac_SP, Jac; |
---|
| 991 | |
---|
| 992 | if( VarNr == 0 ) return; |
---|
| 993 | if (useJacobian == JAC_OFF) return; |
---|
| 994 | |
---|
| 995 | if (useLang != MATLAB_LANG) /* Matlab generates an additional file per function */ |
---|
| 996 | UseFile( jacobianFile ); |
---|
| 997 | |
---|
| 998 | Jac_SP = DefFnc( "Jac_SP", 4, |
---|
| 999 | "the Jacobian of Variables in sparse matrix representation"); |
---|
| 1000 | Jac = DefFnc( "Jac", 4, "the Jacobian of Variables"); |
---|
| 1001 | |
---|
| 1002 | if( useJacSparse ) |
---|
| 1003 | FunctionBegin( Jac_SP, V, F, RCT, JVS ); |
---|
| 1004 | else |
---|
| 1005 | FunctionBegin( Jac, V, F, RCT, JV ); |
---|
| 1006 | |
---|
| 1007 | if (useLang == MATLAB_LANG) { |
---|
| 1008 | switch (useJacobian) { |
---|
| 1009 | case JAC_ROW: |
---|
| 1010 | ExternDeclare(IROW); ExternDeclare(ICOL); |
---|
| 1011 | break; |
---|
| 1012 | case JAC_LU_ROW: |
---|
| 1013 | ExternDeclare(LU_IROW); ExternDeclare(LU_ICOL); |
---|
| 1014 | break; |
---|
| 1015 | } |
---|
| 1016 | } |
---|
| 1017 | |
---|
| 1018 | /* Each nonzero entry of B now counts its rank */ |
---|
| 1019 | nonzeros_B = 0; |
---|
| 1020 | for ( i=0; i<EqnNr; i++ ) |
---|
| 1021 | for ( j=0; j<SpcNr; j++ ) |
---|
| 1022 | if ( structB[i][j] != 0 ) { |
---|
| 1023 | nonzeros_B++; |
---|
| 1024 | structB[i][j] = nonzeros_B; |
---|
| 1025 | } |
---|
| 1026 | |
---|
| 1027 | if ( (useLang==C_LANG)||(useLang==F77_LANG)||(useLang==F90_LANG) ) { |
---|
| 1028 | NewLines(1); |
---|
| 1029 | WriteComment("Local variables"); |
---|
| 1030 | /* DeclareConstant( NTMPB, ascii( nonzeros_B ) ); */ |
---|
| 1031 | varTable[ NTMPB ] -> value = nonzeros_B; |
---|
| 1032 | Declare( BV ); |
---|
| 1033 | } |
---|
| 1034 | |
---|
| 1035 | NewLines(1); |
---|
| 1036 | |
---|
| 1037 | for ( i=0; i<EqnNr; i++ ) { |
---|
| 1038 | for ( j=0; j<VarNr; j++ ) { |
---|
| 1039 | if ( Stoich_Left[j][i] != 0 ) { |
---|
| 1040 | prod = Mul( RConst( i ), Const( Stoich_Left[j][i] ) ); |
---|
| 1041 | for (l = 0; l < VarNr; l++) { |
---|
| 1042 | m = (int)Stoich_Left[l][i] - (l==j); |
---|
| 1043 | for (k = 1; k <= m; k++ ) |
---|
| 1044 | prod = Mul( prod, Elm( V, l ) ); |
---|
| 1045 | } |
---|
| 1046 | for ( ; l < SpcNr; l++) |
---|
| 1047 | for (k = 1; k <= (int)Stoich_Left[l][i]; k++ ) |
---|
| 1048 | prod = Mul( prod, Elm( F, l - VarNr ) ); |
---|
| 1049 | /* Comment the B */ |
---|
| 1050 | WriteComment("B(%d) = dA(%d)/dV(%d)",Index(structB[i][j]-1),Index(i),Index(j)); |
---|
| 1051 | Assign( Elm( BV, structB[i][j]-1 ), prod ); |
---|
| 1052 | } |
---|
| 1053 | } |
---|
| 1054 | } |
---|
| 1055 | |
---|
| 1056 | nElm = 0; |
---|
| 1057 | NewLines(1); |
---|
| 1058 | WriteComment("Construct the Jacobian terms from B's"); |
---|
| 1059 | |
---|
| 1060 | if ( useJacSparse ) { |
---|
| 1061 | for (i = 0; i < VarNr; i++) { |
---|
| 1062 | for (j = 0; j < VarNr; j++) { |
---|
| 1063 | if( LUstructJ[i][j] ) { |
---|
| 1064 | sum = Const(0); |
---|
| 1065 | for (k = 0; k < EqnNr; k++) { |
---|
| 1066 | if( Stoich[i][k]*structB[k][j] != 0 ) |
---|
| 1067 | sum = Add( sum, Mul( Const( Stoich[i][k] ), Elm( BV, structB[k][j]-1 ) ) ); |
---|
| 1068 | } |
---|
| 1069 | /* Comment the B */ |
---|
| 1070 | WriteComment("JVS(%d) = Jac_FULL(%d,%d)", |
---|
| 1071 | Index(nElm),Index(i),Index(j)); |
---|
| 1072 | Assign( Elm( JVS, nElm ), sum ); |
---|
| 1073 | nElm++; |
---|
| 1074 | } else { |
---|
| 1075 | if( i == j ) { |
---|
| 1076 | Assign( Elm( JVS, nElm ), Const(0) ); |
---|
| 1077 | nElm++; |
---|
| 1078 | } |
---|
| 1079 | } |
---|
| 1080 | } |
---|
| 1081 | } |
---|
| 1082 | } else { /* full Jacobian */ |
---|
| 1083 | for (i = 0; i < VarNr; i++) { |
---|
| 1084 | for (j = 0; j < VarNr; j++) { |
---|
| 1085 | if( structJ[i][j] ) { |
---|
| 1086 | sum = Const(0); |
---|
| 1087 | for (k = 0; k < EqnNr; k++) { |
---|
| 1088 | if( Stoich[i][k]*structB[k][j] != 0 ) |
---|
| 1089 | sum = Add( sum, Mul( Const( Stoich[i][k] ), Elm( BV, structB[k][j]-1 ) ) ); |
---|
| 1090 | } |
---|
| 1091 | Assign( Elm( JV, i, j ), sum ); |
---|
| 1092 | } |
---|
| 1093 | } |
---|
| 1094 | } |
---|
| 1095 | } |
---|
| 1096 | |
---|
| 1097 | if (useLang == MATLAB_LANG) { |
---|
| 1098 | switch (useJacobian) { |
---|
| 1099 | case JAC_ROW: |
---|
| 1100 | MATLAB_Inline("\n JVS = sparse(IROW,ICOL,JVS,%d,%d);\n",VarNr,VarNr); |
---|
| 1101 | break; |
---|
| 1102 | case JAC_LU_ROW: |
---|
| 1103 | MATLAB_Inline("\n JVS = sparse(LU_IROW,LU_ICOL,JVS,%d,%d);\n",VarNr,VarNr); |
---|
| 1104 | break; |
---|
| 1105 | } |
---|
| 1106 | } |
---|
| 1107 | |
---|
| 1108 | if( useJacSparse ) |
---|
| 1109 | FunctionEnd( Jac_SP ); |
---|
| 1110 | else |
---|
| 1111 | FunctionEnd( Jac ); |
---|
| 1112 | |
---|
| 1113 | FreeVariable( Jac_SP ); |
---|
| 1114 | FreeVariable( Jac ); |
---|
| 1115 | |
---|
| 1116 | } |
---|
| 1117 | |
---|
| 1118 | |
---|
| 1119 | |
---|
| 1120 | |
---|
| 1121 | |
---|
| 1122 | |
---|
| 1123 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1124 | void GenerateHessian() |
---|
| 1125 | /* Unlike Hess, this function deffers the sparse Data structure generation */ |
---|
| 1126 | { |
---|
| 1127 | int i, j, k; |
---|
| 1128 | int used; |
---|
| 1129 | int l, m, i1, i2, nElm; |
---|
| 1130 | int F_Hess, F_Hess_VEC, F_HessTR_VEC; |
---|
| 1131 | int *coeff_j, *coeff_i1, *coeff_i2; |
---|
| 1132 | int Djv_isElm; |
---|
| 1133 | |
---|
| 1134 | if ( VarNr == 0 ) return; |
---|
| 1135 | |
---|
| 1136 | if (useLang != MATLAB_LANG) /* Matlab generates an additional file per function */ |
---|
| 1137 | UseFile( hessianFile ); |
---|
| 1138 | |
---|
| 1139 | /* Calculate the number of nonzero terms of the form d^2 A(j)/ ( d v(i1) d v(i2) )*/ |
---|
| 1140 | nElm = 0; |
---|
| 1141 | for(j=0; j<EqnNr; j++) |
---|
| 1142 | for (i1 = 0; i1 < VarNr; i1++) |
---|
| 1143 | for (i2 = i1; i2 < VarNr; i2++) |
---|
| 1144 | if (i1==i2) { |
---|
| 1145 | if (Stoich_Left[i1][j]>=2) |
---|
| 1146 | nElm++; |
---|
| 1147 | } else { /* i1 != i2 */ |
---|
| 1148 | if ( (Stoich_Left[i1][j]>=1)&&(Stoich_Left[i2][j]>=1) ) |
---|
| 1149 | nElm++; |
---|
| 1150 | } |
---|
| 1151 | |
---|
| 1152 | /* Allocate temporary index arrays */ |
---|
| 1153 | coeff_j = AllocIntegerVector(nElm, "coeff_j in GenerateHess"); |
---|
| 1154 | coeff_i1 = AllocIntegerVector(nElm, "coeff_i1 in GenerateHess"); |
---|
| 1155 | coeff_i2 = AllocIntegerVector(nElm, "coeff_i2 in GenerateHess"); |
---|
| 1156 | |
---|
| 1157 | /* Fill in temporary index arrays */ |
---|
| 1158 | nElm = 0; |
---|
| 1159 | for(j=0; j<EqnNr; j++) |
---|
| 1160 | for (i1 = 0; i1 < VarNr; i1++) |
---|
| 1161 | for (i2 = i1; i2 < VarNr; i2++) |
---|
| 1162 | if (i1==i2) { |
---|
| 1163 | if (Stoich_Left[i1][j]>=2) { |
---|
| 1164 | coeff_j[nElm] = j; coeff_i1[nElm] = i1; coeff_i2[nElm] = i2; |
---|
| 1165 | nElm++; |
---|
| 1166 | } |
---|
| 1167 | } else { /* i1 != i2 */ |
---|
| 1168 | if ( (Stoich_Left[i1][j]>=1)&&(Stoich_Left[i2][j]>=1) ) { |
---|
| 1169 | coeff_j[nElm] = j; coeff_i1[nElm] = i1; coeff_i2[nElm] = i2; |
---|
| 1170 | nElm++; |
---|
| 1171 | } |
---|
| 1172 | } |
---|
| 1173 | /* Number of nonzero terms of the form d^2 f(i)/ ( d v(i1) d v(i2) ) */ |
---|
| 1174 | Hess_NZ = 0; |
---|
| 1175 | for (i = 0; i < VarNr; i++) |
---|
| 1176 | for (i1 = 0; i1 < VarNr; i1++) |
---|
| 1177 | for (i2 = i1; i2 < VarNr; i2++) { |
---|
| 1178 | Djv_isElm = 0; |
---|
| 1179 | for (j = 0; j < EqnNr; j++) |
---|
| 1180 | if ( Stoich[i][j] != 0 ) |
---|
| 1181 | for (k = 0; k < nElm; k++) |
---|
| 1182 | if ( (coeff_j[k]==j) && (coeff_i1[k]==i1) |
---|
| 1183 | && (coeff_i2[k]==i2) ) { |
---|
| 1184 | Djv_isElm = 1; |
---|
| 1185 | } |
---|
| 1186 | if (Djv_isElm == 1) Hess_NZ++ ; |
---|
| 1187 | } /* for i, i1, i2 */ |
---|
| 1188 | if (useDeclareValues) |
---|
| 1189 | varTable[ NHESS ] -> value = max( Hess_NZ, 1 ); |
---|
| 1190 | |
---|
| 1191 | /* Allocate temporary index arrays */ |
---|
| 1192 | iHess_i = AllocIntegerVector(Hess_NZ, "iHess_i in GenerateHess"); |
---|
| 1193 | iHess_j = AllocIntegerVector(Hess_NZ, "iHess_j in GenerateHess"); |
---|
| 1194 | iHess_k = AllocIntegerVector(Hess_NZ, "iHess_k in GenerateHess"); |
---|
| 1195 | |
---|
| 1196 | F_Hess = DefFnc( "Hessian", 4, "function for Hessian (Jac derivative w.r.t. variables)"); |
---|
| 1197 | FunctionBegin( F_Hess, V, F, RCT, HESS ); |
---|
| 1198 | |
---|
| 1199 | WriteComment("--------------------------------------------------------"); |
---|
| 1200 | WriteComment("Note: HESS is represented in coordinate sparse format: "); |
---|
| 1201 | WriteComment(" HESS(m) = d^2 f_i / dv_j dv_k = d Jac_{i,j} / dv_k"); |
---|
| 1202 | WriteComment(" where i = IHESS_I(m), j = IHESS_J(m), k = IHESS_K(m)."); |
---|
| 1203 | WriteComment("--------------------------------------------------------"); |
---|
| 1204 | WriteComment("Note: d^2 f_i / dv_j dv_k = d^2 f_i / dv_k dv_j, "); |
---|
| 1205 | WriteComment(" therefore only the terms d^2 f_i / dv_j dv_k"); |
---|
| 1206 | WriteComment(" with j <= k are computed and stored in HESS."); |
---|
| 1207 | WriteComment("--------------------------------------------------------"); |
---|
| 1208 | |
---|
| 1209 | if ( (useLang==C_LANG)||(useLang==F77_LANG)||(useLang==F90_LANG) ) { |
---|
| 1210 | NewLines(1); |
---|
| 1211 | WriteComment("Local variables"); |
---|
| 1212 | /* DeclareConstant( NTMPD2A, ascii( max( nElm, 1 ) ) ); */ |
---|
| 1213 | varTable[ NTMPD2A ] -> value = max( nElm, 1 ); |
---|
| 1214 | Declare( D2A ); |
---|
| 1215 | } |
---|
| 1216 | |
---|
| 1217 | NewLines(1); |
---|
| 1218 | WriteComment("Computation of the second derivatives of equation rates"); |
---|
| 1219 | |
---|
| 1220 | /* Generate d^2 A(j)/ ( d v(i1) d v(i2) )*/ |
---|
| 1221 | nElm = 0; |
---|
| 1222 | for(j=0; j<EqnNr; j++) |
---|
| 1223 | for (i1 = 0; i1 < VarNr; i1++) |
---|
| 1224 | for (i2 = i1; i2 < VarNr; i2++) { |
---|
| 1225 | |
---|
| 1226 | if (i1==i2) { |
---|
| 1227 | |
---|
| 1228 | if (Stoich_Left[i1][j]>=2) { |
---|
| 1229 | prod = RConst( j ); |
---|
| 1230 | for (i = 0; i < i1; i++) |
---|
| 1231 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1232 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 1233 | prod = Mul( prod, Const( Stoich_Left[i1][j] ) ); |
---|
| 1234 | prod = Mul( prod, Const( Stoich_Left[i1][j]-1 ) ); |
---|
| 1235 | for (k = 1; k <= (int)Stoich_Left[i1][j]-2; k++ ) |
---|
| 1236 | prod = Mul( prod, Elm( V, i1 ) ); |
---|
| 1237 | for (i = i1+1; i < VarNr; i++) |
---|
| 1238 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1239 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 1240 | for ( ; i < SpcNr; i++) |
---|
| 1241 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1242 | prod = Mul( prod, Elm( F, i - VarNr ) ); |
---|
| 1243 | /* Comment the D2A */ |
---|
| 1244 | WriteComment("D2A(%d) = d^2 A(%d)/{dV(%d)dV(%d)}",Index(nElm),Index(j),Index(i1),Index(i2)); |
---|
| 1245 | Assign( Elm( D2A, nElm ), prod ); |
---|
| 1246 | nElm++; |
---|
| 1247 | } /* if (Stoich_Left[i1][j]>=2) */ |
---|
| 1248 | |
---|
| 1249 | } else { /* i1 != i2 */ |
---|
| 1250 | if ( (Stoich_Left[i1][j]>=1)&&(Stoich_Left[i2][j]>=1) ) { |
---|
| 1251 | prod = RConst( j ); |
---|
| 1252 | for (i = 0; i < i1; i++) |
---|
| 1253 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1254 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 1255 | prod = Mul( prod, Const( Stoich_Left[i1][j] ) ); |
---|
| 1256 | for (k = 1; k <= (int)Stoich_Left[i1][j]-1; k++ ) |
---|
| 1257 | prod = Mul( prod, Elm( V, i1 ) ); |
---|
| 1258 | for (i = i1+1; i < i2; i++) |
---|
| 1259 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1260 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 1261 | prod = Mul( prod, Const( Stoich_Left[i2][j] ) ); |
---|
| 1262 | for (k = 1; k <= (int)Stoich_Left[i2][j]-1; k++ ) |
---|
| 1263 | prod = Mul( prod, Elm( V, i2 ) ); |
---|
| 1264 | for (i = i2+1; i < VarNr; i++) |
---|
| 1265 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1266 | prod = Mul( prod, Elm( V, i ) ); |
---|
| 1267 | for ( ; i < SpcNr; i++) |
---|
| 1268 | for (k = 1; k <= (int)Stoich_Left[i][j]; k++ ) |
---|
| 1269 | prod = Mul( prod, Elm( F, i - VarNr ) ); |
---|
| 1270 | /* Comment the D2A */ |
---|
| 1271 | WriteComment("D2A(%d) = d^2 A(%d) / dV(%d)dV(%d)", |
---|
| 1272 | Index(nElm),Index(j),Index(i1),Index(i2)); |
---|
| 1273 | Assign( Elm( D2A, nElm ), prod ); |
---|
| 1274 | nElm++; |
---|
| 1275 | } /* if ( (Stoich_Left[i1][j]>=1)&&(Stoich_Left[i2][j]>=1) ) */ |
---|
| 1276 | } /* if i1==i2 */ |
---|
| 1277 | |
---|
| 1278 | } /* for j, i1, i2 */ |
---|
| 1279 | |
---|
| 1280 | NewLines(1); |
---|
| 1281 | WriteComment("Computation of the Jacobian derivative"); |
---|
| 1282 | |
---|
| 1283 | /* Generate d^2 f(i)/ ( d v(i1) d v(i2) ) */ |
---|
| 1284 | Hess_NZ = 0; |
---|
| 1285 | for (i = 0; i < VarNr; i++) |
---|
| 1286 | for (i1 = 0; i1 < VarNr; i1++) |
---|
| 1287 | for (i2 = i1; i2 < VarNr; i2++) { |
---|
| 1288 | sum = Const(0); |
---|
| 1289 | Djv_isElm = 0; |
---|
| 1290 | for (j = 0; j < EqnNr; j++) |
---|
| 1291 | if ( Stoich[i][j] != 0 ) |
---|
| 1292 | for (k = 0; k < nElm; k++) |
---|
| 1293 | if ( (coeff_j[k]==j) && (coeff_i1[k]==i1) |
---|
| 1294 | && (coeff_i2[k]==i2) ) { |
---|
| 1295 | sum = Add( sum, |
---|
| 1296 | Mul( Const( Stoich[i][j] ), Elm( D2A, k ) ) ); |
---|
| 1297 | Djv_isElm = 1; |
---|
| 1298 | } |
---|
| 1299 | if (Djv_isElm == 1) { |
---|
| 1300 | WriteComment("HESS(%d) = d^2 Vdot(%d)/{dV(%d)dV(%d)} = d^2 Vdot(%d)/{dV(%d)dV(%d)}", |
---|
| 1301 | Index(Hess_NZ),Index(i),Index(i1),Index(i2),Index(i),Index(i2),Index(i1)); |
---|
| 1302 | Assign( Elm( HESS, Hess_NZ ), sum ); |
---|
| 1303 | iHess_i[ Hess_NZ ] = i; |
---|
| 1304 | iHess_j[ Hess_NZ ] = i1; |
---|
| 1305 | iHess_k[ Hess_NZ ] = i2; |
---|
| 1306 | Hess_NZ++; |
---|
| 1307 | } |
---|
| 1308 | |
---|
| 1309 | } /* for i, i1, i2 */ |
---|
| 1310 | |
---|
| 1311 | |
---|
| 1312 | /* free temporary index arrays */ |
---|
| 1313 | free(coeff_j); free(coeff_i1); free(coeff_i2); |
---|
| 1314 | |
---|
| 1315 | MATLAB_Inline("\n HESS = HESS(:);"); |
---|
| 1316 | |
---|
| 1317 | FunctionEnd( F_Hess ); |
---|
| 1318 | |
---|
| 1319 | FreeVariable( F_Hess ); |
---|
| 1320 | |
---|
| 1321 | |
---|
| 1322 | F_HessTR_VEC = DefFnc( "HessTR_Vec", 4, "Hessian transposed times user vectors"); |
---|
| 1323 | FunctionBegin( F_HessTR_VEC, HESS, U1, U2, HTU ); |
---|
| 1324 | WriteComment("Compute the vector HTU =(Hess x U2)^T * U1 = d (Jac^T*U1)/d Var * U2 "); |
---|
| 1325 | |
---|
| 1326 | for (i=0; i<VarNr; i++) { |
---|
| 1327 | sum = Const(0); |
---|
| 1328 | for (k=0; k<Hess_NZ; k++) { |
---|
| 1329 | if (iHess_j[k]==i) |
---|
| 1330 | sum = Add( sum, |
---|
| 1331 | Mul( Elm( HESS, k ), |
---|
| 1332 | Mul( Elm( U1, iHess_i[k] ), Elm( U2, iHess_k[k] ) ) ) ); |
---|
| 1333 | else if (iHess_k[k]==i) |
---|
| 1334 | sum = Add( sum, |
---|
| 1335 | Mul( Elm( HESS, k ), |
---|
| 1336 | Mul( Elm( U1, iHess_i[k] ), Elm( U2, iHess_j[k] ) ) ) ); |
---|
| 1337 | } |
---|
| 1338 | Assign( Elm( HTU, i ), sum ); |
---|
| 1339 | } |
---|
| 1340 | |
---|
| 1341 | MATLAB_Inline("\n HTU = HTU(:);"); |
---|
| 1342 | |
---|
| 1343 | FunctionEnd( F_HessTR_VEC ); |
---|
| 1344 | FreeVariable( F_HessTR_VEC ); |
---|
| 1345 | |
---|
| 1346 | |
---|
| 1347 | F_Hess_VEC = DefFnc( "Hess_Vec", 4, "Hessian times user vectors"); |
---|
| 1348 | FunctionBegin( F_HessTR_VEC, HESS, U1, U2, HU ); |
---|
| 1349 | WriteComment("Compute the vector HU =(Hess x U2) * U1 = d (Jac*U1)/d Var * U2 "); |
---|
| 1350 | |
---|
| 1351 | for (i=0; i<VarNr; i++) { |
---|
| 1352 | sum = Const(0); |
---|
| 1353 | for (m=0; m<Hess_NZ; m++) { |
---|
| 1354 | if (iHess_i[m]==i) { |
---|
| 1355 | j = iHess_j[m]; |
---|
| 1356 | k = iHess_k[m]; |
---|
| 1357 | if (j==k) { |
---|
| 1358 | sum = Add( sum, |
---|
| 1359 | Mul( Elm( HESS, m ), |
---|
| 1360 | Mul( Elm( U1, k ), Elm( U2, k ) ) ) ); |
---|
| 1361 | } |
---|
| 1362 | else { |
---|
| 1363 | /* This is the optimized code. It can lead to problems when |
---|
| 1364 | splitting for continuation lines. Therefore we |
---|
| 1365 | use the non-optimized code below the comment. |
---|
| 1366 | sum = Add( sum, |
---|
| 1367 | Mul( Elm( HESS, m ), |
---|
| 1368 | Add( Mul( Elm( U1, j ), Elm( U2, k ) ), |
---|
| 1369 | Mul( Elm( U1, k ), Elm( U2, j ) ) ) ) ); */ |
---|
| 1370 | sum = Add( sum, |
---|
| 1371 | Mul( Elm( HESS, m ), |
---|
| 1372 | Mul( Elm( U1, j ), Elm( U2, k ) ) ) ); |
---|
| 1373 | sum = Add( sum, |
---|
| 1374 | Mul( Elm( HESS, m ), |
---|
| 1375 | Mul( Elm( U1, k ), Elm( U2, j ) ) ) ); |
---|
| 1376 | } |
---|
| 1377 | } |
---|
| 1378 | } |
---|
| 1379 | Assign( Elm( HU, i ), sum ); |
---|
| 1380 | } |
---|
| 1381 | |
---|
| 1382 | MATLAB_Inline("\n HU = HU(:);"); |
---|
| 1383 | |
---|
| 1384 | FunctionEnd( F_Hess_VEC ); |
---|
| 1385 | FreeVariable( F_Hess_VEC ); |
---|
| 1386 | } |
---|
| 1387 | |
---|
| 1388 | |
---|
| 1389 | |
---|
| 1390 | |
---|
| 1391 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1392 | void GenerateHessianSparseData() |
---|
| 1393 | { |
---|
| 1394 | int k; |
---|
| 1395 | |
---|
| 1396 | |
---|
| 1397 | UseFile( sparse_hessFile ); |
---|
| 1398 | |
---|
| 1399 | NewLines(1); |
---|
| 1400 | |
---|
| 1401 | |
---|
| 1402 | WriteComment("Hessian Sparse Data"); |
---|
| 1403 | WriteComment(""); |
---|
| 1404 | |
---|
| 1405 | F77_Inline("%6sBLOCK DATA HESSIAN_SPARSE_DATA\n", " " ); |
---|
| 1406 | F77_Inline("%6sINCLUDE '%s_Sparse.h'", " ", rootFileName); |
---|
| 1407 | F77_Inline("%6sINTEGER i", " "); |
---|
| 1408 | /* F90_Inline(" USE %s_Sparse", rootFileName); */ |
---|
| 1409 | |
---|
| 1410 | if( (useLang==F77_LANG)||(useLang==F90_LANG)||(useLang==MATLAB_LANG) ) { |
---|
| 1411 | for (k=0; k<Hess_NZ; k++) { |
---|
| 1412 | iHess_i[k]++; iHess_j[k]++; iHess_k[k]++; |
---|
| 1413 | } |
---|
| 1414 | } |
---|
| 1415 | |
---|
| 1416 | InitDeclare( IHESS_I, Hess_NZ, (void*)iHess_i ); |
---|
| 1417 | InitDeclare( IHESS_J, Hess_NZ, (void*)iHess_j ); |
---|
| 1418 | InitDeclare( IHESS_K, Hess_NZ, (void*)iHess_k ); |
---|
| 1419 | |
---|
| 1420 | if( (useLang==F77_LANG)||(useLang==F90_LANG) ) { |
---|
| 1421 | for (k=0; k<Hess_NZ; k++) { |
---|
| 1422 | iHess_i[k]--; iHess_j[k]--; iHess_k[k]--; |
---|
| 1423 | } |
---|
| 1424 | } |
---|
| 1425 | |
---|
| 1426 | NewLines(1); |
---|
| 1427 | F77_Inline( "%6sEND\n\n", " " ); |
---|
| 1428 | |
---|
| 1429 | } |
---|
| 1430 | |
---|
| 1431 | |
---|
| 1432 | |
---|
| 1433 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1434 | void GenerateHessianSparseHeader() |
---|
| 1435 | { |
---|
| 1436 | UseFile( sparse_dataFile ); |
---|
| 1437 | |
---|
| 1438 | CommonName = "HESSDATA"; |
---|
| 1439 | |
---|
| 1440 | NewLines(1); |
---|
| 1441 | WriteComment(" ----------> Sparse Hessian Data"); |
---|
| 1442 | NewLines(1); |
---|
| 1443 | |
---|
| 1444 | ExternDeclare( IHESS_I ); |
---|
| 1445 | ExternDeclare( IHESS_J ); |
---|
| 1446 | ExternDeclare( IHESS_K ); |
---|
| 1447 | |
---|
| 1448 | NewLines(1); |
---|
| 1449 | } |
---|
| 1450 | |
---|
| 1451 | |
---|
| 1452 | |
---|
| 1453 | |
---|
| 1454 | |
---|
| 1455 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1456 | void GenerateStoicmSparseData() |
---|
| 1457 | { |
---|
| 1458 | int i,j,k, nnz_stoicm; |
---|
| 1459 | /* |
---|
| 1460 | int irow_stoicm[MAX_SPECIES*MAX_EQN]; |
---|
| 1461 | int ccol_stoicm[MAX_EQN+2]; |
---|
| 1462 | int icol_stoicm[MAX_SPECIES*MAX_EQN]; |
---|
| 1463 | double stoicm[MAX_SPECIES*MAX_EQN]; |
---|
| 1464 | */ |
---|
| 1465 | |
---|
| 1466 | int *irow_stoicm; |
---|
| 1467 | int *ccol_stoicm; |
---|
| 1468 | int *icol_stoicm; |
---|
| 1469 | double *stoicm; |
---|
| 1470 | |
---|
| 1471 | /* Compute the sparsity structure and allocate data structure vectors */ |
---|
| 1472 | nnz_stoicm = 0; |
---|
| 1473 | for (j=0; j<EqnNr; j++) |
---|
| 1474 | for (i=0; i<VarNr; i++) |
---|
| 1475 | if ( Stoich[i][j] != 0.0 ) |
---|
| 1476 | nnz_stoicm++; |
---|
| 1477 | if ( (irow_stoicm=(int*)calloc(nnz_stoicm+2,sizeof(int)) ) == NULL ) |
---|
| 1478 | FatalError(-30,"GenerateStoicmSparseData: Cannot allocate irow_stoicm"); |
---|
| 1479 | if ( (ccol_stoicm=(int*)calloc(EqnNr+2,sizeof(int)) ) == NULL ) |
---|
| 1480 | FatalError(-30,"GenerateStoicmSparseData: Cannot allocate ccol_stoicm"); |
---|
| 1481 | if ( (icol_stoicm=(int*)calloc(nnz_stoicm+2,sizeof(int)) ) == NULL ) |
---|
| 1482 | FatalError(-30,"GenerateStoicmSparseData: Cannot allocate icol_stoicm"); |
---|
| 1483 | if ( (stoicm=(double*)calloc(nnz_stoicm+2,sizeof(double)) ) == NULL ) |
---|
| 1484 | FatalError(-30,"GenerateStoicmSparseData: Cannot allocate stoicm"); |
---|
| 1485 | |
---|
| 1486 | UseFile( sparse_stoicmFile ); |
---|
| 1487 | |
---|
| 1488 | nnz_stoicm = 0; |
---|
| 1489 | for (j=0; j<EqnNr; j++) { |
---|
| 1490 | ccol_stoicm[ j ] = nnz_stoicm; |
---|
| 1491 | for (i=0; i<VarNr; i++) { |
---|
| 1492 | if ( Stoich[i][j] != 0 ) { |
---|
| 1493 | irow_stoicm[ nnz_stoicm ] = i; |
---|
| 1494 | icol_stoicm[ nnz_stoicm ] = j; |
---|
| 1495 | stoicm[ nnz_stoicm ] = Stoich[i][j]; |
---|
| 1496 | nnz_stoicm++; |
---|
| 1497 | } |
---|
| 1498 | } |
---|
| 1499 | } |
---|
| 1500 | ccol_stoicm[ EqnNr ] = nnz_stoicm; |
---|
| 1501 | |
---|
| 1502 | if( (useLang==F77_LANG)||(useLang==F90_LANG) ) { |
---|
| 1503 | for (k=0; k<nnz_stoicm; k++) { |
---|
| 1504 | irow_stoicm[k]++; icol_stoicm[k]++; |
---|
| 1505 | } |
---|
| 1506 | for (k=0; k<=EqnNr; k++) { |
---|
| 1507 | ccol_stoicm[k]++; |
---|
| 1508 | } |
---|
| 1509 | } |
---|
| 1510 | |
---|
| 1511 | |
---|
| 1512 | NewLines(1); |
---|
| 1513 | WriteComment(" Stoichiometric Matrix in Compressed Column Sparse Format"); |
---|
| 1514 | NewLines(1); |
---|
| 1515 | F77_Inline("%6sBLOCK DATA STOICM_MATRIX\n", " " ); |
---|
| 1516 | F77_Inline("%6sINCLUDE '%s_Sparse.h'", " ", rootFileName); |
---|
| 1517 | F77_Inline("%6sINTEGER i", " "); |
---|
| 1518 | /* F90_Inline(" USE %s_Sparse", rootFileName); */ |
---|
| 1519 | InitDeclare( CCOL_STOICM, EqnNr+1, (void*)ccol_stoicm ); |
---|
| 1520 | InitDeclare( IROW_STOICM, nnz_stoicm, (void*)irow_stoicm ); |
---|
| 1521 | InitDeclare( ICOL_STOICM, nnz_stoicm, (void*)icol_stoicm ); |
---|
| 1522 | InitDeclare( STOICM, nnz_stoicm, (void*)stoicm ); |
---|
| 1523 | NewLines(1); |
---|
| 1524 | F77_Inline( "%6sEND\n\n", " " ); |
---|
| 1525 | |
---|
| 1526 | |
---|
| 1527 | UseFile( param_headerFile ); |
---|
| 1528 | CommonName = "GDATA"; |
---|
| 1529 | NewLines(1); |
---|
| 1530 | DeclareConstant( NSTOICM, ascii( max( nnz_stoicm, 1 ) ) ); |
---|
| 1531 | |
---|
| 1532 | /* Free data structure vectors */ |
---|
| 1533 | free(irow_stoicm); free(ccol_stoicm); free(icol_stoicm); free(stoicm); |
---|
| 1534 | } |
---|
| 1535 | |
---|
| 1536 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1537 | void GenerateStoicmSparseHeader() |
---|
| 1538 | { |
---|
| 1539 | UseFile( sparse_dataFile ); |
---|
| 1540 | |
---|
| 1541 | NewLines(1); |
---|
| 1542 | WriteComment(" ----------> Sparse Stoichiometric Matrix"); |
---|
| 1543 | NewLines(1); |
---|
| 1544 | CommonName = "STOICM_VALUES"; |
---|
| 1545 | ExternDeclare( STOICM ); |
---|
| 1546 | CommonName = "STOICM_DATA"; |
---|
| 1547 | ExternDeclare( IROW_STOICM ); |
---|
| 1548 | ExternDeclare( CCOL_STOICM ); |
---|
| 1549 | ExternDeclare( ICOL_STOICM ); |
---|
| 1550 | NewLines(1); |
---|
| 1551 | |
---|
| 1552 | NewLines(1); |
---|
| 1553 | WriteComment(" ----------> Sparse Data for Jacobian of Reactant Products"); |
---|
| 1554 | NewLines(1); |
---|
| 1555 | CommonName = "JVRP"; |
---|
| 1556 | ExternDeclare( ICOL_JVRP ); |
---|
| 1557 | ExternDeclare( IROW_JVRP ); |
---|
| 1558 | ExternDeclare( CROW_JVRP ); |
---|
| 1559 | NewLines(1); |
---|
| 1560 | |
---|
| 1561 | } |
---|
| 1562 | |
---|
| 1563 | |
---|
| 1564 | |
---|
| 1565 | |
---|
| 1566 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1567 | void GenerateJacVect() |
---|
| 1568 | { |
---|
| 1569 | int i, j, nElm; |
---|
| 1570 | int Jac_VEC; |
---|
| 1571 | int Jac_SP_VEC; |
---|
| 1572 | |
---|
| 1573 | if( useLang == MATLAB_LANG ) return; |
---|
| 1574 | |
---|
| 1575 | if( VarNr == 0 ) return; |
---|
| 1576 | |
---|
| 1577 | UseFile( jacobianFile ); |
---|
| 1578 | Jac_VEC = DefFnc( "Jac_Vec", 3, |
---|
| 1579 | "function for sparse multiplication: square Jacobian times vector"); |
---|
| 1580 | Jac_SP_VEC = DefFnc( "Jac_SP_Vec", 3, |
---|
| 1581 | "function for sparse multiplication: sparse Jacobian times vector"); |
---|
| 1582 | |
---|
| 1583 | if ( useJacSparse ) { |
---|
| 1584 | FunctionBegin( Jac_SP_VEC, JVS, UV, JUV ); |
---|
| 1585 | nElm = 0; |
---|
| 1586 | for( i = 0; i < VarNr; i++) { |
---|
| 1587 | sum = Const(0); |
---|
| 1588 | for( j = 0; j < VarNr; j++ ) |
---|
| 1589 | if( LUstructJ[i][j] ) { |
---|
| 1590 | if( structJ[i][j] != 0 ) |
---|
| 1591 | sum = Add( sum, Mul( Elm( JVS, nElm ), Elm( UV, j ) ) ); |
---|
| 1592 | nElm++; |
---|
| 1593 | } |
---|
| 1594 | Assign( Elm( JUV, i ), sum ); |
---|
| 1595 | } |
---|
| 1596 | FunctionEnd( Jac_SP_VEC ); |
---|
| 1597 | } |
---|
| 1598 | |
---|
| 1599 | else { |
---|
| 1600 | FunctionBegin( Jac_VEC, JV, UV, JUV ); |
---|
| 1601 | for( i = 0; i < VarNr; i++) { |
---|
| 1602 | sum = Const(0); |
---|
| 1603 | for( j = 0; j < VarNr; j++ ) |
---|
| 1604 | if( structJ[i][j] != 0 ) |
---|
| 1605 | sum = Add( sum, Mul( Elm( JV, i, j ), Elm( UV, j ) ) ); |
---|
| 1606 | Assign( Elm( JUV, i ), sum ); |
---|
| 1607 | } |
---|
| 1608 | FunctionEnd( Jac_VEC ); |
---|
| 1609 | } |
---|
| 1610 | |
---|
| 1611 | FreeVariable( Jac_VEC ); |
---|
| 1612 | FreeVariable( Jac_SP_VEC ); |
---|
| 1613 | } |
---|
| 1614 | |
---|
| 1615 | |
---|
| 1616 | |
---|
| 1617 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1618 | void GenerateJacTRVect() |
---|
| 1619 | { |
---|
| 1620 | int i, j, nElm; |
---|
| 1621 | int JacTR_VEC; |
---|
| 1622 | int JacTR_SP_VEC; |
---|
| 1623 | int **TmpStruct; |
---|
| 1624 | |
---|
| 1625 | if( useLang == MATLAB_LANG ) return; |
---|
| 1626 | |
---|
| 1627 | if ( VarNr == 0 ) return; |
---|
| 1628 | |
---|
| 1629 | UseFile( jacobianFile ); |
---|
| 1630 | |
---|
| 1631 | JacTR_VEC = DefFnc( "JacTR_Vec", 3, |
---|
| 1632 | "sparse multiplication: square Jacobian transposed times vector"); |
---|
| 1633 | JacTR_SP_VEC = DefFnc( "JacTR_SP_Vec", 3, |
---|
| 1634 | "sparse multiplication: sparse Jacobian transposed times vector"); |
---|
| 1635 | |
---|
| 1636 | if ( useJacSparse ) { |
---|
| 1637 | |
---|
| 1638 | /* The temporary array of structure */ |
---|
| 1639 | TmpStruct = AllocIntegerMatrix( VarNr, VarNr, "TmpStruct in GenerateJacTRVect" ); |
---|
| 1640 | |
---|
| 1641 | nElm = 0; |
---|
| 1642 | for( i = 0; i < VarNr; i++) |
---|
| 1643 | for( j = 0; j < VarNr; j++ ) |
---|
| 1644 | if( LUstructJ[i][j] ) { |
---|
| 1645 | TmpStruct[i][j] = nElm; |
---|
| 1646 | nElm++; |
---|
| 1647 | } |
---|
| 1648 | |
---|
| 1649 | FunctionBegin( JacTR_SP_VEC, JVS, UV, JTUV ); |
---|
| 1650 | nElm = 0; |
---|
| 1651 | for( i = 0; i < VarNr; i++) { |
---|
| 1652 | sum = Const(0); |
---|
| 1653 | for( j = 0; j < VarNr; j++ ) |
---|
| 1654 | if( structJ[j][i] != 0 ) |
---|
| 1655 | sum = Add( sum, Mul( Elm( JVS, TmpStruct[j][i] ), Elm( UV, j ) ) ); |
---|
| 1656 | Assign( Elm( JTUV, i ), sum ); |
---|
| 1657 | } |
---|
| 1658 | FunctionEnd( JacTR_SP_VEC ); |
---|
| 1659 | |
---|
| 1660 | /* Free the temporary array of structure */ |
---|
| 1661 | FreeIntegerMatrix( TmpStruct, VarNr, VarNr ); |
---|
| 1662 | |
---|
| 1663 | } /* useJacSparse*/ |
---|
| 1664 | |
---|
| 1665 | else { |
---|
| 1666 | FunctionBegin( JacTR_VEC, JV, UV, JTUV ); |
---|
| 1667 | for( i = 0; i < VarNr; i++) { |
---|
| 1668 | sum = Const(0); |
---|
| 1669 | for( j = 0; j < VarNr; j++ ) |
---|
| 1670 | if( structJ[j][i] != 0 ) |
---|
| 1671 | sum = Add( sum, Mul( Elm( JV, j, i ), Elm( UV, j ) ) ); |
---|
| 1672 | Assign( Elm( JTUV, i ), sum ); |
---|
| 1673 | } |
---|
| 1674 | FunctionEnd( JacTR_VEC ); |
---|
| 1675 | } |
---|
| 1676 | |
---|
| 1677 | FreeVariable( JacTR_VEC ); |
---|
| 1678 | FreeVariable( JacTR_SP_VEC ); |
---|
| 1679 | } |
---|
| 1680 | |
---|
| 1681 | |
---|
| 1682 | |
---|
| 1683 | |
---|
| 1684 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1685 | void GenerateSparseUtil() |
---|
| 1686 | { |
---|
| 1687 | int SUTIL; |
---|
| 1688 | |
---|
| 1689 | if ( useLang == MATLAB_LANG ) return; |
---|
| 1690 | |
---|
| 1691 | UseFile( linalgFile ); |
---|
| 1692 | |
---|
| 1693 | SUTIL = DefFnc( "SPARSE_UTIL", 0, "SPARSE utility functions"); |
---|
| 1694 | CommentFunctionBegin( SUTIL ); |
---|
| 1695 | |
---|
| 1696 | IncludeCode( "%s/util/sutil", Home ); |
---|
| 1697 | |
---|
| 1698 | CommentFunctionEnd( SUTIL ); |
---|
| 1699 | FreeVariable( SUTIL ); |
---|
| 1700 | } |
---|
| 1701 | |
---|
| 1702 | |
---|
| 1703 | |
---|
| 1704 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1705 | void GenerateBlas() |
---|
| 1706 | { |
---|
| 1707 | int BLAS; |
---|
| 1708 | |
---|
| 1709 | if ( useLang == MATLAB_LANG ) return; |
---|
| 1710 | |
---|
| 1711 | UseFile( linalgFile ); |
---|
| 1712 | |
---|
| 1713 | BLAS = DefFnc( "BLAS_UTIL", 0, "BLAS-LIKE utility functions"); |
---|
| 1714 | CommentFunctionBegin( BLAS ); |
---|
| 1715 | |
---|
| 1716 | IncludeCode( "%s/util/blas", Home ); |
---|
| 1717 | |
---|
| 1718 | CommentFunctionEnd( BLAS ); |
---|
| 1719 | FreeVariable( BLAS ); |
---|
| 1720 | } |
---|
| 1721 | |
---|
| 1722 | |
---|
| 1723 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1724 | void GenerateDFunDRcoeff() |
---|
| 1725 | { |
---|
| 1726 | |
---|
| 1727 | UseFile( stoichiomFile ); |
---|
| 1728 | |
---|
| 1729 | NewLines(1); |
---|
| 1730 | WriteComment("Begin Derivative w.r.t. Rate Coefficients"); |
---|
| 1731 | NewLines(1); |
---|
| 1732 | |
---|
| 1733 | IncludeCode( "%s/util/dFun_dRcoeff", Home ); |
---|
| 1734 | |
---|
| 1735 | NewLines(1); |
---|
| 1736 | WriteComment("End Derivative w.r.t. Rate Coefficients"); |
---|
| 1737 | NewLines(1); |
---|
| 1738 | |
---|
| 1739 | } |
---|
| 1740 | |
---|
| 1741 | |
---|
| 1742 | |
---|
| 1743 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1744 | void GenerateDJacDRcoeff() |
---|
| 1745 | { |
---|
| 1746 | |
---|
| 1747 | UseFile( stoichiomFile ); |
---|
| 1748 | |
---|
| 1749 | NewLines(1); |
---|
| 1750 | WriteComment("Begin Jacobian Derivative w.r.t. Rate Coefficients"); |
---|
| 1751 | NewLines(1); |
---|
| 1752 | |
---|
| 1753 | IncludeCode( "%s/util/dJac_dRcoeff", Home ); |
---|
| 1754 | |
---|
| 1755 | NewLines(1); |
---|
| 1756 | WriteComment("End Jacobian Derivative w.r.t. Rate Coefficients"); |
---|
| 1757 | NewLines(1); |
---|
| 1758 | |
---|
| 1759 | } |
---|
| 1760 | |
---|
| 1761 | |
---|
| 1762 | |
---|
| 1763 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1764 | void GenerateSolve() |
---|
| 1765 | { |
---|
| 1766 | int i, j; |
---|
| 1767 | int SOLVE; |
---|
| 1768 | int *irow; |
---|
| 1769 | int *icol; |
---|
| 1770 | int *crow; |
---|
| 1771 | int *diag; |
---|
| 1772 | int nElm; |
---|
| 1773 | int ibgn, iend; |
---|
| 1774 | int useLangOld; |
---|
| 1775 | int dim; |
---|
| 1776 | |
---|
| 1777 | if( useLang == MATLAB_LANG ) return; |
---|
| 1778 | |
---|
| 1779 | /* Allocate local arrays for dimension dim */ |
---|
| 1780 | dim = VarNr+2; |
---|
| 1781 | irow = AllocIntegerVector( dim*dim, "irow in GenerateSolve" ); |
---|
| 1782 | icol = AllocIntegerVector( dim*dim, "icol in GenerateSolve" ); |
---|
| 1783 | crow = AllocIntegerVector( dim, "crow in GenerateSolve" ); |
---|
| 1784 | diag = AllocIntegerVector( dim, "diag in GenerateSolve" ); |
---|
| 1785 | |
---|
| 1786 | useLangOld = useLang; |
---|
| 1787 | useLang = C_LANG; |
---|
| 1788 | nElm = NonZero( LU, 0, VarNr, irow, icol, crow, diag ); |
---|
| 1789 | useLang = useLangOld; |
---|
| 1790 | |
---|
| 1791 | UseFile( linalgFile ); |
---|
| 1792 | |
---|
| 1793 | SOLVE = DefFnc( "KppSolve", 2, "sparse back substitution"); |
---|
| 1794 | FunctionBegin( SOLVE, JVS, X ); |
---|
| 1795 | |
---|
| 1796 | for( i = 0; i < VarNr; i++) { |
---|
| 1797 | ibgn = crow[i]; |
---|
| 1798 | iend = diag[i]; |
---|
| 1799 | if( ibgn <= iend ) { |
---|
| 1800 | sum = Elm( X, i ); |
---|
| 1801 | if ( ibgn < iend ) { |
---|
| 1802 | for( j = ibgn; j < iend; j++ ) |
---|
| 1803 | sum = Sub( sum, Mul( Elm( JVS, j ), Elm( X, icol[j] ) ) ); |
---|
| 1804 | Assign( Elm( X, i ), sum ); |
---|
| 1805 | } |
---|
| 1806 | } |
---|
| 1807 | } |
---|
| 1808 | |
---|
| 1809 | for( i = VarNr-1; i >=0; i--) { |
---|
| 1810 | ibgn = diag[i] + 1; |
---|
| 1811 | iend = crow[i+1]; |
---|
| 1812 | sum = Elm( X, i ); |
---|
| 1813 | for( j = ibgn; j < iend; j++ ) |
---|
| 1814 | sum = Sub( sum, Mul( Elm( JVS, j ), Elm( X, icol[j] ) ) ); |
---|
| 1815 | sum = Div( sum, Elm( JVS, diag[i] ) ); |
---|
| 1816 | Assign( Elm( X, i ), sum ); |
---|
| 1817 | } |
---|
| 1818 | |
---|
| 1819 | FunctionEnd( SOLVE ); |
---|
| 1820 | FreeVariable( SOLVE ); |
---|
| 1821 | |
---|
| 1822 | /* Free Local Arrays */ |
---|
| 1823 | free(irow); |
---|
| 1824 | free(icol); |
---|
| 1825 | free(crow); |
---|
| 1826 | free(diag); |
---|
| 1827 | } |
---|
| 1828 | |
---|
| 1829 | |
---|
| 1830 | |
---|
| 1831 | |
---|
| 1832 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1833 | void GenerateTRSolve() |
---|
| 1834 | { |
---|
| 1835 | int i, j; |
---|
| 1836 | int SOLVETR; |
---|
| 1837 | int *irow; |
---|
| 1838 | int *icol; |
---|
| 1839 | int *crow; |
---|
| 1840 | int *diag; |
---|
| 1841 | int nElm; |
---|
| 1842 | int ibgn, iend; |
---|
| 1843 | int useLangOld; |
---|
| 1844 | int **pos; |
---|
| 1845 | int dim; |
---|
| 1846 | |
---|
| 1847 | if( useLang == MATLAB_LANG ) return; |
---|
| 1848 | |
---|
| 1849 | /* Allocate local arrays for dimension dim */ |
---|
| 1850 | dim = VarNr+2; |
---|
| 1851 | irow = AllocIntegerVector( dim*dim, "irow in GenerateTRSolve" ); |
---|
| 1852 | icol = AllocIntegerVector( dim*dim, "icol in GenerateTRSolve" ); |
---|
| 1853 | crow = AllocIntegerVector( dim, "crow in GenerateTRSolve" ); |
---|
| 1854 | diag = AllocIntegerVector( dim, "diag in GenerateTRSolve" ); |
---|
| 1855 | pos = AllocIntegerMatrix( dim+1, dim+1, "pos in GenerateTRSolve"); |
---|
| 1856 | |
---|
| 1857 | useLangOld = useLang; |
---|
| 1858 | useLang = C_LANG; |
---|
| 1859 | nElm = NonZero( LU, 0, VarNr, irow, icol, crow, diag ); |
---|
| 1860 | useLang = useLangOld; |
---|
| 1861 | |
---|
| 1862 | UseFile( linalgFile ); |
---|
| 1863 | |
---|
| 1864 | SOLVETR = DefFnc( "KppSolveTR", 3, "sparse, transposed back substitution"); |
---|
| 1865 | FunctionBegin( SOLVETR, JVS, X, XX ); |
---|
| 1866 | for( i = 0; i < VarNr; i++) { |
---|
| 1867 | for( j = 0; j < VarNr; j++) |
---|
| 1868 | pos[i][j]=-1; |
---|
| 1869 | } |
---|
| 1870 | for( i = 0; i < VarNr; i++) { |
---|
| 1871 | ibgn = crow[i]; |
---|
| 1872 | iend = diag[i]; |
---|
| 1873 | if( ibgn <= iend ) { |
---|
| 1874 | if ( ibgn < iend ) { |
---|
| 1875 | for( j = ibgn; j < iend; j++ ) |
---|
| 1876 | pos[icol[j]][i]=j; |
---|
| 1877 | } |
---|
| 1878 | } |
---|
| 1879 | } |
---|
| 1880 | |
---|
| 1881 | for( i = VarNr-1; i >=0; i--) { |
---|
| 1882 | ibgn = diag[i] + 1; |
---|
| 1883 | iend = crow[i+1]; |
---|
| 1884 | for( j = ibgn; j < iend; j++ ) |
---|
| 1885 | pos[icol[j]][i]=j; |
---|
| 1886 | pos[i][i]=diag[i]; |
---|
| 1887 | } |
---|
| 1888 | |
---|
| 1889 | for( i = 0; i<VarNr; i++) { |
---|
| 1890 | sum = Elm( X, i ); |
---|
| 1891 | for (j=0; j<i; j++ ){ |
---|
| 1892 | if(pos[i][j] >= 0) { |
---|
| 1893 | sum=Sub( sum, Mul ( Elm(JVS,pos[i][j] ), Elm( XX, j ) ) ); |
---|
| 1894 | } |
---|
| 1895 | } |
---|
| 1896 | sum=Div( sum, Elm(JVS, diag[i] ) ); |
---|
| 1897 | Assign( Elm( XX, i ), sum ); |
---|
| 1898 | } |
---|
| 1899 | for( i = VarNr-1; i >=0; i--) { |
---|
| 1900 | sum = Elm( XX, i ); |
---|
| 1901 | for (j=i+1; j<VarNr; j++) { |
---|
| 1902 | if(pos[i][j] >= 0) { |
---|
| 1903 | sum=Sub( sum, Mul ( Elm(JVS,pos[i][j] ), Elm( XX, j ) ) ); |
---|
| 1904 | } |
---|
| 1905 | } |
---|
| 1906 | Assign( Elm( XX, i ), sum ); |
---|
| 1907 | } |
---|
| 1908 | |
---|
| 1909 | FunctionEnd( SOLVETR ); |
---|
| 1910 | FreeVariable( SOLVETR ); |
---|
| 1911 | /* Free Local Arrays */ |
---|
| 1912 | free(irow); free(icol); free(crow); free(diag); |
---|
| 1913 | FreeIntegerMatrix(pos, dim+1, dim+1); |
---|
| 1914 | } |
---|
| 1915 | |
---|
| 1916 | |
---|
| 1917 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1918 | void GenerateRateLaws() |
---|
| 1919 | { |
---|
| 1920 | |
---|
| 1921 | UseFile( rateFile ); |
---|
| 1922 | |
---|
| 1923 | NewLines(1); |
---|
| 1924 | WriteComment("Begin Rate Law Functions from KPP_HOME/util/UserRateLaws"); |
---|
| 1925 | NewLines(1); |
---|
| 1926 | IncludeCode( "%s/util/UserRateLaws", Home ); |
---|
| 1927 | NewLines(1); |
---|
| 1928 | WriteComment("End Rate Law Functions from KPP_HOME/util/UserRateLaws"); |
---|
| 1929 | NewLines(1); |
---|
| 1930 | |
---|
| 1931 | NewLines(1); |
---|
| 1932 | WriteComment("Begin INLINED Rate Law Functions"); |
---|
| 1933 | NewLines(1); |
---|
| 1934 | |
---|
| 1935 | switch( useLang ) { |
---|
| 1936 | case C_LANG: bprintf( InlineCode[ C_RATES ].code ); |
---|
| 1937 | break; |
---|
| 1938 | case F77_LANG: bprintf( InlineCode[ F77_RATES ].code ); |
---|
| 1939 | break; |
---|
| 1940 | case F90_LANG: bprintf( InlineCode[ F90_RATES ].code ); |
---|
| 1941 | break; |
---|
| 1942 | case MATLAB_LANG: bprintf( InlineCode[ MATLAB_RATES ].code ); |
---|
| 1943 | break; |
---|
| 1944 | } |
---|
| 1945 | FlushBuf(); |
---|
| 1946 | |
---|
| 1947 | NewLines(1); |
---|
| 1948 | WriteComment("End INLINED Rate Law Functions"); |
---|
| 1949 | NewLines(1); |
---|
| 1950 | |
---|
| 1951 | |
---|
| 1952 | } |
---|
| 1953 | |
---|
| 1954 | |
---|
| 1955 | |
---|
| 1956 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1957 | void GenerateUpdateSun() |
---|
| 1958 | { |
---|
| 1959 | int UPDATE_SUN; |
---|
| 1960 | |
---|
| 1961 | UseFile( rateFile ); |
---|
| 1962 | |
---|
| 1963 | UPDATE_SUN = DefFnc( "Update_SUN", 0, "update SUN light using TIME"); |
---|
| 1964 | CommentFunctionBegin( UPDATE_SUN ); |
---|
| 1965 | |
---|
| 1966 | IncludeCode( "%s/util/UpdateSun", Home ); |
---|
| 1967 | |
---|
| 1968 | CommentFunctionEnd( UPDATE_SUN ); |
---|
| 1969 | FreeVariable( UPDATE_SUN ); |
---|
| 1970 | } |
---|
| 1971 | |
---|
| 1972 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1973 | void GenerateUpdateRconst() |
---|
| 1974 | { |
---|
| 1975 | int i; |
---|
| 1976 | int UPDATE_RCONST; |
---|
| 1977 | |
---|
| 1978 | UseFile( rateFile ); |
---|
| 1979 | |
---|
| 1980 | UPDATE_RCONST = DefFnc( "Update_RCONST", 0, "function to update rate constants"); |
---|
| 1981 | |
---|
| 1982 | FunctionBegin( UPDATE_RCONST ); |
---|
| 1983 | F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); |
---|
| 1984 | MATLAB_Inline("global SUN TEMP RCONST"); |
---|
| 1985 | |
---|
| 1986 | if ( (useLang==F77_LANG) ) |
---|
| 1987 | IncludeCode( "%s/util/UserRateLaws_FcnHeader", Home ); |
---|
| 1988 | |
---|
| 1989 | NewLines(1); |
---|
| 1990 | |
---|
| 1991 | NewLines(1); |
---|
| 1992 | WriteComment("Begin INLINED RCONST"); |
---|
| 1993 | NewLines(1); |
---|
| 1994 | |
---|
| 1995 | switch( useLang ) { |
---|
| 1996 | case C_LANG: bprintf( InlineCode[ C_RCONST ].code ); |
---|
| 1997 | break; |
---|
| 1998 | case F77_LANG: bprintf( InlineCode[ F77_RCONST ].code ); |
---|
| 1999 | break; |
---|
| 2000 | case F90_LANG: bprintf( InlineCode[ F90_RCONST ].code ); |
---|
| 2001 | break; |
---|
| 2002 | case MATLAB_LANG: bprintf( InlineCode[ MATLAB_RCONST ].code ); |
---|
| 2003 | break; |
---|
| 2004 | } |
---|
| 2005 | FlushBuf(); |
---|
| 2006 | |
---|
| 2007 | NewLines(1); |
---|
| 2008 | WriteComment("End INLINED RCONST"); |
---|
| 2009 | NewLines(1); |
---|
| 2010 | |
---|
| 2011 | for( i = 0; i < EqnNr; i++) { |
---|
| 2012 | if( kr[i].type == EXPRESION ) |
---|
| 2013 | Assign( Elm( RCONST, i ), Elm( KR, kr[i].val.st ) ); |
---|
| 2014 | if( kr[i].type == PHOTO ) |
---|
| 2015 | Assign( Elm( RCONST, i ), Elm( KR, kr[i].val.st ) ); |
---|
| 2016 | /* mz_rs_20050117+ */ |
---|
| 2017 | if ( kr[i].type == NUMBER ) { |
---|
| 2018 | F90_Inline("! RCONST(%d) = constant rate coefficient", i+1); |
---|
| 2019 | /* WriteComment("Constant rate coefficient (value inlined in the code):"); */ |
---|
| 2020 | /* Assign( Elm( RCONST, i ), Const( kr[i].val.f ) ); */ |
---|
| 2021 | } |
---|
| 2022 | /* mz_rs_20050117- */ |
---|
| 2023 | } |
---|
| 2024 | |
---|
| 2025 | MATLAB_Inline(" RCONST = RCONST(:);"); |
---|
| 2026 | |
---|
| 2027 | FunctionEnd( UPDATE_RCONST ); |
---|
| 2028 | FreeVariable( UPDATE_RCONST ); |
---|
| 2029 | } |
---|
| 2030 | |
---|
| 2031 | |
---|
| 2032 | |
---|
| 2033 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2034 | void GenerateUpdatePhoto() |
---|
| 2035 | { |
---|
| 2036 | int i; |
---|
| 2037 | int UPDATE_PHOTO; |
---|
| 2038 | |
---|
| 2039 | UseFile( rateFile ); |
---|
| 2040 | |
---|
| 2041 | UPDATE_PHOTO = DefFnc( "Update_PHOTO", 0, "function to update photolytical rate constants"); |
---|
| 2042 | |
---|
| 2043 | FunctionBegin( UPDATE_PHOTO ); |
---|
| 2044 | F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); |
---|
| 2045 | F90_Inline(" USE %s_Global", rootFileName); |
---|
| 2046 | MATLAB_Inline("global SUN TEMP RCONST"); |
---|
| 2047 | |
---|
| 2048 | NewLines(1); |
---|
| 2049 | |
---|
| 2050 | for( i = 0; i < EqnNr; i++) { |
---|
| 2051 | if( kr[i].type == PHOTO ) |
---|
| 2052 | Assign( Elm( RCONST, i ), Elm( KR, kr[i].val.st ) ); |
---|
| 2053 | } |
---|
| 2054 | |
---|
| 2055 | FunctionEnd( UPDATE_PHOTO ); |
---|
| 2056 | FreeVariable( UPDATE_PHOTO ); |
---|
| 2057 | } |
---|
| 2058 | |
---|
| 2059 | |
---|
| 2060 | |
---|
| 2061 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2062 | void GenerateIntegrator() |
---|
| 2063 | { |
---|
| 2064 | int TIN, TOUT, INTEGRATE; |
---|
| 2065 | |
---|
| 2066 | UseFile( integratorFile ); |
---|
| 2067 | |
---|
| 2068 | TIN = DefElm( "TIN", real, "Start Time for Integration"); |
---|
| 2069 | TOUT = DefElm( "TOUT", real, "End Time for Integration"); |
---|
| 2070 | INTEGRATE = DefFnc( "INTEGRATE", 2, "Integrator routine"); |
---|
| 2071 | CommentFunctionBegin( INTEGRATE, TIN, TOUT ); |
---|
| 2072 | |
---|
| 2073 | if( strchr( integrator, '/' ) ) |
---|
| 2074 | IncludeCode( integrator ); |
---|
| 2075 | else |
---|
| 2076 | IncludeCode( "%s/int/%s", Home, integrator ); |
---|
| 2077 | |
---|
| 2078 | CommentFunctionEnd( INTEGRATE ); |
---|
| 2079 | FreeVariable( INTEGRATE ); |
---|
| 2080 | } |
---|
| 2081 | |
---|
| 2082 | |
---|
| 2083 | |
---|
| 2084 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2085 | void GenerateDriver() |
---|
| 2086 | { |
---|
| 2087 | int MAIN; |
---|
| 2088 | |
---|
| 2089 | UseFile( driverFile ); |
---|
| 2090 | |
---|
| 2091 | MAIN = DefFnc( "MAIN", 0, "Main program - driver routine"); |
---|
| 2092 | CommentFunctionBegin( MAIN ); |
---|
| 2093 | |
---|
| 2094 | if( strchr( driver, '/' ) ) |
---|
| 2095 | IncludeCode( driver ); |
---|
| 2096 | else |
---|
| 2097 | IncludeCode( "%s/drv/%s", Home, driver ); |
---|
| 2098 | |
---|
| 2099 | CommentFunctionEnd( MAIN ); |
---|
| 2100 | FreeVariable( MAIN ); |
---|
| 2101 | } |
---|
| 2102 | |
---|
| 2103 | |
---|
| 2104 | |
---|
| 2105 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2106 | void GenerateUtil() |
---|
| 2107 | { |
---|
| 2108 | int UTIL; |
---|
| 2109 | |
---|
| 2110 | /* if (useLang == MATLAB_LANG) return; */ |
---|
| 2111 | |
---|
| 2112 | UseFile( utilFile ); |
---|
| 2113 | NewLines(1); |
---|
| 2114 | WriteComment("User INLINED Utility Functions"); |
---|
| 2115 | |
---|
| 2116 | switch( useLang ) { |
---|
| 2117 | case C_LANG: bprintf( InlineCode[ C_UTIL ].code ); |
---|
| 2118 | break; |
---|
| 2119 | case F77_LANG: bprintf( InlineCode[ F77_UTIL ].code ); |
---|
| 2120 | break; |
---|
| 2121 | case F90_LANG: bprintf( InlineCode[ F90_UTIL ].code ); |
---|
| 2122 | break; |
---|
| 2123 | case MATLAB_LANG:bprintf( InlineCode[ MATLAB_UTIL ].code ); |
---|
| 2124 | break; |
---|
| 2125 | } |
---|
| 2126 | FlushBuf(); |
---|
| 2127 | |
---|
| 2128 | NewLines(1); |
---|
| 2129 | WriteComment("End INLINED Utility Functions"); |
---|
| 2130 | NewLines(1); |
---|
| 2131 | |
---|
| 2132 | WriteComment("Utility Functions from KPP_HOME/util/util"); |
---|
| 2133 | UTIL = DefFnc( "UTIL", 0, "Utility functions"); |
---|
| 2134 | CommentFunctionBegin( UTIL); |
---|
| 2135 | |
---|
| 2136 | IncludeCode( "%s/util/util", Home ); |
---|
| 2137 | |
---|
| 2138 | if ((useLang == F90_LANG) && (useEqntags==1)) { |
---|
| 2139 | IncludeCode( "%s/util/tag2num", Home ); |
---|
| 2140 | } |
---|
| 2141 | |
---|
| 2142 | WriteComment("End Utility Functions from KPP_HOME/util/util"); |
---|
| 2143 | CommentFunctionEnd( UTIL ); |
---|
| 2144 | FreeVariable( UTIL ); |
---|
| 2145 | } |
---|
| 2146 | |
---|
| 2147 | |
---|
| 2148 | |
---|
| 2149 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2150 | void GenerateParamHeader() |
---|
| 2151 | { |
---|
| 2152 | int spc; |
---|
| 2153 | int i; |
---|
| 2154 | char name[20]; |
---|
| 2155 | int offs; |
---|
| 2156 | int mxyz; |
---|
| 2157 | |
---|
| 2158 | int j,dummy_species; |
---|
| 2159 | |
---|
| 2160 | /* ----------> First declaration of constants */ |
---|
| 2161 | UseFile( param_headerFile ); |
---|
| 2162 | |
---|
| 2163 | NewLines(1); |
---|
| 2164 | DeclareConstant( NSPEC, ascii( max(SpcNr, 1) ) ); |
---|
| 2165 | DeclareConstant( NVAR, ascii( max(VarNr, 1) ) ); |
---|
| 2166 | DeclareConstant( NVARACT, ascii( max(VarActiveNr, 1) ) ); |
---|
| 2167 | DeclareConstant( NFIX, ascii( max(FixNr, 1) ) ); |
---|
| 2168 | DeclareConstant( NREACT, ascii( max(EqnNr, 1) ) ); |
---|
| 2169 | DeclareConstant( NVARST, ascii( VarStartNr ) ); |
---|
| 2170 | DeclareConstant( NFIXST, ascii( FixStartNr ) ); |
---|
| 2171 | DeclareConstant( NONZERO, ascii( max(Jac_NZ, 1) ) ); |
---|
| 2172 | DeclareConstant( LU_NONZERO, ascii( max(LU_Jac_NZ, 1) ) ); |
---|
| 2173 | DeclareConstant( CNVAR, ascii( VarNr+1 ) ); |
---|
| 2174 | if ( useStoicmat ) { |
---|
| 2175 | DeclareConstant( CNEQN, ascii( EqnNr+1 ) ); |
---|
| 2176 | } |
---|
| 2177 | if ( useHessian ) { |
---|
| 2178 | DeclareConstant( NHESS, ascii( max(Hess_NZ, 1) ) ); |
---|
| 2179 | } |
---|
| 2180 | |
---|
| 2181 | DeclareConstant( NLOOKAT, ascii( nlookat ) ); |
---|
| 2182 | DeclareConstant( NMONITOR, ascii( nmoni ) ); |
---|
| 2183 | DeclareConstant( NMASS, ascii( nmass ) ); |
---|
| 2184 | |
---|
| 2185 | /* DeclareConstant( PI, "3.14159265358979" ); */ |
---|
| 2186 | |
---|
| 2187 | NewLines(1); |
---|
| 2188 | WriteComment("Index declaration for variable species in C and VAR"); |
---|
| 2189 | WriteComment(" VAR(ind_spc) = C(ind_spc)"); |
---|
| 2190 | NewLines(1); |
---|
| 2191 | for( i = 0; i < VarNr; i++) { |
---|
| 2192 | sprintf( name, "ind_%s", SpeciesTable[ Code[i] ].name ); |
---|
| 2193 | spc = DefConst( name, INT, 0 ); |
---|
| 2194 | DeclareConstant( spc, ascii( Index(i) ) ); |
---|
| 2195 | FreeVariable( spc ); |
---|
| 2196 | } |
---|
| 2197 | |
---|
| 2198 | NewLines(1); |
---|
| 2199 | WriteComment("Index declaration for fixed species in C"); |
---|
| 2200 | WriteComment(" C(ind_spc)"); |
---|
| 2201 | NewLines(1); |
---|
| 2202 | for( i = 0; i < FixNr; i++) { |
---|
| 2203 | sprintf( name, "ind_%s", SpeciesTable[ Code[i + VarNr] ].name ); |
---|
| 2204 | spc = DefConst( name, INT, 0 ); |
---|
| 2205 | DeclareConstant( spc, ascii( Index(i+VarNr) ) ); |
---|
| 2206 | FreeVariable( spc ); |
---|
| 2207 | } |
---|
| 2208 | |
---|
| 2209 | if (useDummyindex==1) { |
---|
| 2210 | NewLines(1); |
---|
| 2211 | WriteComment("Index declaration for dummy species"); |
---|
| 2212 | NewLines(1); |
---|
| 2213 | for( i = 0; i < MAX_SPECIES; i++) { |
---|
| 2214 | if (SpeciesTable[i].type == 0) continue; |
---|
| 2215 | dummy_species = 1; |
---|
| 2216 | for( j = 0; j < MAX_SPECIES; j++) |
---|
| 2217 | if (Code[j] == i) dummy_species = 0; |
---|
| 2218 | if (dummy_species) { |
---|
| 2219 | sprintf( name, "ind_%s", SpeciesTable[i].name ); |
---|
| 2220 | spc = DefConst( name, INT, 0 ); |
---|
| 2221 | DeclareConstant( spc, ascii( 0 ) ); |
---|
| 2222 | FreeVariable( spc ); |
---|
| 2223 | } |
---|
| 2224 | } |
---|
| 2225 | } |
---|
| 2226 | |
---|
| 2227 | NewLines(1); |
---|
| 2228 | WriteComment("Index declaration for fixed species in FIX"); |
---|
| 2229 | WriteComment(" FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)"); |
---|
| 2230 | NewLines(1); |
---|
| 2231 | for( i = 0; i < FixNr; i++) { |
---|
| 2232 | sprintf( name, "indf_%s", SpeciesTable[ Code[i + VarNr] ].name ); |
---|
| 2233 | spc = DefConst( name, INT, 0 ); |
---|
| 2234 | DeclareConstant( spc, ascii( Index(i) ) ); |
---|
| 2235 | FreeVariable( spc ); |
---|
| 2236 | } |
---|
| 2237 | } |
---|
| 2238 | |
---|
| 2239 | |
---|
| 2240 | |
---|
| 2241 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2242 | void GenerateGlobalHeader() |
---|
| 2243 | { |
---|
| 2244 | int spc; |
---|
| 2245 | int i; |
---|
| 2246 | char name[20]; |
---|
| 2247 | int offs; |
---|
| 2248 | int mxyz; |
---|
| 2249 | |
---|
| 2250 | UseFile( global_dataFile ); |
---|
| 2251 | |
---|
| 2252 | CommonName = "GDATA"; |
---|
| 2253 | |
---|
| 2254 | NewLines(1); |
---|
| 2255 | WriteComment("Declaration of global variables"); |
---|
| 2256 | NewLines(1); |
---|
| 2257 | |
---|
| 2258 | /* ExternDeclare( C_DEFAULT ); */ |
---|
| 2259 | |
---|
| 2260 | ExternDeclare( C ); |
---|
| 2261 | |
---|
| 2262 | if( useLang == F77_LANG ) { |
---|
| 2263 | |
---|
| 2264 | Declare( VAR ); |
---|
| 2265 | Declare( FIX ); |
---|
| 2266 | WriteComment("VAR, FIX are chunks of array C"); |
---|
| 2267 | F77_Inline(" EQUIVALENCE( %s(%d),%s(1) )", |
---|
| 2268 | varTable[C]->name, 1, varTable[VAR]->name ); |
---|
| 2269 | if ( FixNr > 0 ) { /* mz_rs_20050121 */ |
---|
| 2270 | F77_Inline(" EQUIVALENCE( %s(%d),%s(1) )", |
---|
| 2271 | varTable[C]->name, VarNr+1, varTable[FIX]->name ); |
---|
| 2272 | } |
---|
| 2273 | } |
---|
| 2274 | |
---|
| 2275 | if( useLang == F90_LANG ) { |
---|
| 2276 | ExternDeclare( VAR ); |
---|
| 2277 | ExternDeclare( FIX ); |
---|
| 2278 | WriteComment("VAR, FIX are chunks of array C"); |
---|
| 2279 | F90_Inline(" EQUIVALENCE( %s(%d),%s(1) )", |
---|
| 2280 | varTable[C]->name, 1, varTable[VAR]->name ); |
---|
| 2281 | if ( FixNr > 0 ) { /* mz_rs_20050121 */ |
---|
| 2282 | F90_Inline(" EQUIVALENCE( %s(%d),%s(1) )", |
---|
| 2283 | varTable[C]->name, VarNr+1, varTable[FIX]->name ); |
---|
| 2284 | } |
---|
| 2285 | } |
---|
| 2286 | |
---|
| 2287 | if( useLang == MATLAB_LANG ) { |
---|
| 2288 | ExternDeclare( VAR ); |
---|
| 2289 | ExternDeclare( FIX ); |
---|
| 2290 | } |
---|
| 2291 | |
---|
| 2292 | C_Inline(" extern %s * %s;", C_types[real], varTable[VAR]->name ); |
---|
| 2293 | C_Inline(" extern %s * %s;", C_types[real], varTable[FIX]->name ); |
---|
| 2294 | |
---|
| 2295 | |
---|
| 2296 | ExternDeclare( RCONST ); |
---|
| 2297 | ExternDeclare( TIME ); |
---|
| 2298 | ExternDeclare( SUN ); |
---|
| 2299 | ExternDeclare( TEMP ); |
---|
| 2300 | ExternDeclare( RTOLS ); |
---|
| 2301 | ExternDeclare( TSTART ); |
---|
| 2302 | ExternDeclare( TEND ); |
---|
| 2303 | ExternDeclare( DT ); |
---|
| 2304 | ExternDeclare( ATOL ); |
---|
| 2305 | ExternDeclare( RTOL ); |
---|
| 2306 | ExternDeclare( STEPMIN ); |
---|
| 2307 | ExternDeclare( STEPMAX ); |
---|
| 2308 | ExternDeclare( CFACTOR ); |
---|
| 2309 | if (useStochastic) |
---|
| 2310 | ExternDeclare( VOLUME ); |
---|
| 2311 | |
---|
| 2312 | CommonName = "INTGDATA"; |
---|
| 2313 | if ( useHessian ) { |
---|
| 2314 | ExternDeclare( DDMTYPE ); |
---|
| 2315 | } |
---|
| 2316 | |
---|
| 2317 | |
---|
| 2318 | if ( (useLang == C_LANG) || (useLang == F77_LANG) ) { |
---|
| 2319 | CommonName = "INTGDATA"; |
---|
| 2320 | ExternDeclare( LOOKAT ); |
---|
| 2321 | ExternDeclare( MONITOR ); |
---|
| 2322 | CommonName = "CHARGDATA"; |
---|
| 2323 | ExternDeclare( SPC_NAMES ); |
---|
| 2324 | ExternDeclare( SMASS ); |
---|
| 2325 | ExternDeclare( EQN_NAMES ); |
---|
| 2326 | ExternDeclare( EQN_TAGS ); |
---|
| 2327 | } |
---|
| 2328 | |
---|
| 2329 | NewLines(1); |
---|
| 2330 | WriteComment("INLINED global variable declarations"); |
---|
| 2331 | |
---|
| 2332 | switch( useLang ) { |
---|
| 2333 | case C_LANG: bprintf( InlineCode[ C_GLOBAL ].code ); |
---|
| 2334 | break; |
---|
| 2335 | case F77_LANG: bprintf( InlineCode[ F77_GLOBAL ].code ); |
---|
| 2336 | break; |
---|
| 2337 | case F90_LANG: bprintf( InlineCode[ F90_GLOBAL ].code ); |
---|
| 2338 | break; |
---|
| 2339 | case MATLAB_LANG: bprintf( InlineCode[ MATLAB_GLOBAL ].code ); |
---|
| 2340 | break; |
---|
| 2341 | } |
---|
| 2342 | FlushBuf(); |
---|
| 2343 | |
---|
| 2344 | NewLines(1); |
---|
| 2345 | WriteComment("INLINED global variable declarations"); |
---|
| 2346 | NewLines(1); |
---|
| 2347 | } |
---|
| 2348 | |
---|
| 2349 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2350 | void WriteSpec( int i, int j ) |
---|
| 2351 | { |
---|
| 2352 | char buf[100]; |
---|
| 2353 | |
---|
| 2354 | if( Reactive[j] ) |
---|
| 2355 | sprintf( buf, "%s (r)", SpeciesTable[ Code[j] ].name ); |
---|
| 2356 | else |
---|
| 2357 | sprintf( buf, "%s (n)", SpeciesTable[ Code[j] ].name ); |
---|
| 2358 | WriteAll("%3d = %-10s", 1 + i, buf ); |
---|
| 2359 | FlushBuf(); |
---|
| 2360 | } |
---|
| 2361 | |
---|
| 2362 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2363 | int EqnStr( int eq, char * buf, float** mat ) |
---|
| 2364 | { |
---|
| 2365 | int spc, first; |
---|
| 2366 | |
---|
| 2367 | /* bugfix if stoichiometric factor is not an integer */ |
---|
| 2368 | int n; |
---|
| 2369 | char s[40]; |
---|
| 2370 | |
---|
| 2371 | first = 1; |
---|
| 2372 | *buf = 0; |
---|
| 2373 | for( spc = 0; spc < SpcNr; spc++ ) |
---|
| 2374 | if( mat[spc][eq] != 0 ) { |
---|
| 2375 | if( ((mat[spc][eq] == 1)||(mat[spc][eq] == -1)) ) { |
---|
| 2376 | sprintf(s, ""); |
---|
| 2377 | } else { |
---|
| 2378 | /* real */ |
---|
| 2379 | /* mz_rs_20050130+ */ |
---|
| 2380 | /* sprintf(s, "%g", mat[spc][eq]); */ |
---|
| 2381 | /* remove the minus sign with fabs(), it will be re-inserted later */ |
---|
| 2382 | sprintf(s, "%g", mat[spc][eq]?mat[spc][eq]:(-mat[spc][eq])); |
---|
| 2383 | /* mz_rs_20050130- */ |
---|
| 2384 | /* remove trailing zeroes */ |
---|
| 2385 | for (n= strlen(s) - 1; n >= 0; n--) |
---|
| 2386 | if (s[n] != '0') break; |
---|
| 2387 | s[n + 1]= '\0'; |
---|
| 2388 | sprintf(s, "%s ", s); |
---|
| 2389 | } |
---|
| 2390 | |
---|
| 2391 | if( first ) { |
---|
| 2392 | if( mat[spc][eq] > 0 ) sprintf(buf, "%s%s", buf, s); |
---|
| 2393 | else sprintf(buf, "%s- %s", buf, s); |
---|
| 2394 | first = 0; |
---|
| 2395 | } else { |
---|
| 2396 | if( mat[spc][eq] > 0 ) sprintf(buf, "%s + %s", buf, s); |
---|
| 2397 | else sprintf(buf, "%s - %s", buf, s); |
---|
| 2398 | } |
---|
| 2399 | sprintf(buf, "%s%s", buf, SpeciesTable[ Code[spc] ].name); |
---|
| 2400 | if (strlen(buf)>MAX_EQNLEN/2) { /* truncate if eqn string too long */ |
---|
| 2401 | sprintf(buf, "%s ... etc.",buf); |
---|
| 2402 | break; |
---|
| 2403 | } |
---|
| 2404 | } |
---|
| 2405 | |
---|
| 2406 | return strlen(buf); |
---|
| 2407 | } |
---|
| 2408 | |
---|
| 2409 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2410 | int EqnString( int eq, char * buf ) |
---|
| 2411 | { |
---|
| 2412 | static int lhs = 0; |
---|
| 2413 | static int rhs = 0; |
---|
| 2414 | |
---|
| 2415 | int i, l; |
---|
| 2416 | char lhsbuf[MAX_EQNLEN], rhsbuf[MAX_EQNLEN]; |
---|
| 2417 | |
---|
| 2418 | if(lhs == 0) for( i = 0; i < EqnNr; i++ ) { |
---|
| 2419 | l = EqnStr( i, lhsbuf, Stoich_Left); |
---|
| 2420 | lhs = (lhs > l) ? lhs : l; |
---|
| 2421 | } |
---|
| 2422 | |
---|
| 2423 | if(rhs == 0) for( i = 0; i < EqnNr; i++ ) { |
---|
| 2424 | l = EqnStr( i, rhsbuf, Stoich_Right); |
---|
| 2425 | rhs = (rhs > l) ? lhs : l; |
---|
| 2426 | } |
---|
| 2427 | |
---|
| 2428 | |
---|
| 2429 | EqnStr( eq, lhsbuf, Stoich_Left); |
---|
| 2430 | EqnStr( eq, rhsbuf, Stoich_Right); |
---|
| 2431 | |
---|
| 2432 | sprintf(buf, "%*s --> %-*s", lhs, lhsbuf, rhs, rhsbuf); |
---|
| 2433 | return strlen(buf); |
---|
| 2434 | } |
---|
| 2435 | |
---|
| 2436 | |
---|
| 2437 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2438 | void GenerateMap() |
---|
| 2439 | { |
---|
| 2440 | int i; |
---|
| 2441 | int dn; |
---|
| 2442 | |
---|
| 2443 | UseFile( mapFile ); |
---|
| 2444 | |
---|
| 2445 | WriteAll("### Options -------------------------------------------\n"); |
---|
| 2446 | NewLines(1); |
---|
| 2447 | if( useAggregate ) WriteAll("FUNCTION - AGGREGATE\n"); |
---|
| 2448 | else WriteAll("FUNCTION - SPLIT\n"); |
---|
| 2449 | switch ( useJacobian ) { |
---|
| 2450 | case JAC_OFF: WriteAll("JACOBIAN - OFF\n"); break; |
---|
| 2451 | case JAC_FULL: WriteAll("JACOBIAN - FULL\n"); break; |
---|
| 2452 | case JAC_LU_ROW: WriteAll("JACOBIAN - SPARSE W/ ACCOUNT FOR LU DECOMPOSITION FILL-IN\n"); break; |
---|
| 2453 | case JAC_ROW: WriteAll("JACOBIAN - SPARSE\n"); break; |
---|
| 2454 | } |
---|
| 2455 | if( useDouble ) WriteAll("DOUBLE - ON\n"); |
---|
| 2456 | else WriteAll("DOUBLE - OFF\n"); |
---|
| 2457 | if( useReorder ) WriteAll("REORDER - ON\n"); |
---|
| 2458 | else WriteAll("REORDER - OFF\n"); |
---|
| 2459 | NewLines(1); |
---|
| 2460 | |
---|
| 2461 | WriteAll("### Parameters ----------------------------------------\n"); |
---|
| 2462 | NewLines(1); |
---|
| 2463 | |
---|
| 2464 | VarStartNr = Index(0); |
---|
| 2465 | FixStartNr = Index(VarNr); |
---|
| 2466 | |
---|
| 2467 | DeclareConstant( NSPEC, ascii( SpcNr ) ); |
---|
| 2468 | DeclareConstant( NVAR, ascii( max( VarNr, 1 ) ) ); |
---|
| 2469 | DeclareConstant( NVARACT, ascii( max( VarActiveNr, 1 ) ) ); |
---|
| 2470 | DeclareConstant( NFIX, ascii( max( FixNr, 1 ) ) ); |
---|
| 2471 | DeclareConstant( NREACT, ascii( EqnNr ) ); |
---|
| 2472 | DeclareConstant( NVARST, ascii( VarStartNr ) ); |
---|
| 2473 | DeclareConstant( NFIXST, ascii( FixStartNr ) ); |
---|
| 2474 | |
---|
| 2475 | NewLines(1); |
---|
| 2476 | WriteAll("### Species -------------------------------------------\n"); |
---|
| 2477 | |
---|
| 2478 | NewLines(1); |
---|
| 2479 | WriteAll("Variable species\n"); |
---|
| 2480 | |
---|
| 2481 | dn = VarNr/3 + 1; |
---|
| 2482 | for( i = 0; i < dn; i++ ) { |
---|
| 2483 | if( i < VarNr ) WriteSpec( i, i ); |
---|
| 2484 | i += dn; if( i < VarNr ) WriteSpec( i, i ); |
---|
| 2485 | i += dn; if( i < VarNr ) WriteSpec( i, i ); |
---|
| 2486 | i -= 2*dn; WriteAll("\n"); |
---|
| 2487 | } |
---|
| 2488 | |
---|
| 2489 | |
---|
| 2490 | NewLines(1); |
---|
| 2491 | WriteAll("Fixed species\n"); |
---|
| 2492 | |
---|
| 2493 | dn = FixNr/3 + 1; |
---|
| 2494 | for( i = 0; i < dn; i++ ) { |
---|
| 2495 | if( i < FixNr ) WriteSpec( i, i + VarNr ); |
---|
| 2496 | i += dn; if( i < FixNr ) WriteSpec( i, i + VarNr ); |
---|
| 2497 | i += dn; if( i < FixNr ) WriteSpec( i, i + VarNr ); |
---|
| 2498 | i -= 2*dn; WriteAll("\n"); |
---|
| 2499 | } |
---|
| 2500 | |
---|
| 2501 | NewLines(1); |
---|
| 2502 | WriteAll("### Subroutines ---------------------------------------\n"); |
---|
| 2503 | NewLines(1); |
---|
| 2504 | FlushBuf(); |
---|
| 2505 | } |
---|
| 2506 | |
---|
| 2507 | |
---|
| 2508 | |
---|
| 2509 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2510 | void GenerateInitialize() |
---|
| 2511 | { |
---|
| 2512 | int i; |
---|
| 2513 | int I, X; |
---|
| 2514 | int INITVAL; |
---|
| 2515 | |
---|
| 2516 | if ( (useLang==C_LANG)||(useLang==F77_LANG)||(useLang==F90_LANG) ) |
---|
| 2517 | UseFile( initFile ); |
---|
| 2518 | |
---|
| 2519 | INITVAL = DefFnc( "Initialize", 0, "function to initialize concentrations"); |
---|
| 2520 | FunctionBegin( INITVAL ); |
---|
| 2521 | F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); |
---|
| 2522 | F90_Inline(" USE %s_Global\n", rootFileName); |
---|
| 2523 | MATLAB_Inline("global CFACTOR VAR FIX NVAR NFIX", rootFileName); |
---|
| 2524 | |
---|
| 2525 | I = DefElm( "i", INT, 0); |
---|
| 2526 | X = DefElm( "x", real, 0); |
---|
| 2527 | Declare( I ); |
---|
| 2528 | Declare( X ); |
---|
| 2529 | |
---|
| 2530 | NewLines(1); |
---|
| 2531 | WriteAssign( varTable[CFACTOR]->name , ascid( (double)cfactor ) ); |
---|
| 2532 | NewLines(1); |
---|
| 2533 | |
---|
| 2534 | Assign( Elm( X ), Mul( Elm( IV, varDefault ), Elm( CFACTOR ) ) ); |
---|
| 2535 | C_Inline(" for( i = 0; i < NVAR; i++ )" ); |
---|
| 2536 | F77_Inline(" DO i = 1, NVAR" ); |
---|
| 2537 | F90_Inline(" DO i = 1, NVAR" ); |
---|
| 2538 | MATLAB_Inline(" for i = 1:NVAR" ); |
---|
| 2539 | ident++; |
---|
| 2540 | Assign( Elm( VAR, -I ), Elm( X ) ); |
---|
| 2541 | ident--; |
---|
| 2542 | F77_Inline(" END DO" ); |
---|
| 2543 | F90_Inline(" END DO" ); |
---|
| 2544 | MATLAB_Inline(" end" ); |
---|
| 2545 | |
---|
| 2546 | |
---|
| 2547 | NewLines(1); |
---|
| 2548 | Assign( Elm( X ), Mul( Elm( IV, fixDefault ), Elm( CFACTOR ) ) ); |
---|
| 2549 | C_Inline(" for( i = 0; i < NFIX; i++ )" ); |
---|
| 2550 | F77_Inline(" DO i = 1, NFIX" ); |
---|
| 2551 | F90_Inline(" DO i = 1, NFIX" ); |
---|
| 2552 | MATLAB_Inline(" for i = 1:NFIX" ); |
---|
| 2553 | ident++; |
---|
| 2554 | Assign( Elm( FIX, -I ), Elm( X ) ); |
---|
| 2555 | ident--; |
---|
| 2556 | F77_Inline(" END DO" ); |
---|
| 2557 | F90_Inline(" END DO" ); |
---|
| 2558 | MATLAB_Inline(" end" ); |
---|
| 2559 | |
---|
| 2560 | |
---|
| 2561 | NewLines(1); |
---|
| 2562 | |
---|
| 2563 | for( i = 0; i < VarNr; i++) { |
---|
| 2564 | if( *SpeciesTable[ Code[i] ].ival == 0 ) continue; |
---|
| 2565 | Assign( Elm( VAR, i ), Mul( |
---|
| 2566 | Elm( IV, SpeciesTable[ Code[i] ].ival ), |
---|
| 2567 | Elm( CFACTOR ) ) ); |
---|
| 2568 | } |
---|
| 2569 | |
---|
| 2570 | |
---|
| 2571 | for( i = 0; i < FixNr; i++) { |
---|
| 2572 | if( *SpeciesTable[ Code[i + VarNr] ].ival == 0 ) continue; |
---|
| 2573 | Assign( Elm( FIX, i ), Mul( |
---|
| 2574 | Elm( IV, SpeciesTable[ Code[i + VarNr] ].ival ), |
---|
| 2575 | Elm( CFACTOR ) ) ); |
---|
| 2576 | } |
---|
| 2577 | |
---|
| 2578 | /* NewLines(1); |
---|
| 2579 | C_Inline(" for( i = 0; i < NSPEC; i++ )" ); |
---|
| 2580 | F77_Inline(" do i = 1, NSPEC" ); |
---|
| 2581 | ident++; |
---|
| 2582 | Assign( Elm( C_DEFAULT, -I ), Elm( C, -I ) ); |
---|
| 2583 | ident--; |
---|
| 2584 | F77_Inline(" end do" ); |
---|
| 2585 | */ |
---|
| 2586 | |
---|
| 2587 | /* mz_rs_20050117+ */ |
---|
| 2588 | WriteComment("constant rate coefficients"); |
---|
| 2589 | for( i = 0; i < EqnNr; i++) { |
---|
| 2590 | if ( kr[i].type == NUMBER ) |
---|
| 2591 | Assign( Elm( RCONST, i ), Const( kr[i].val.f ) ); |
---|
| 2592 | } |
---|
| 2593 | WriteComment("END constant rate coefficients"); |
---|
| 2594 | /* mz_rs_20050117- */ |
---|
| 2595 | |
---|
| 2596 | NewLines(1); |
---|
| 2597 | WriteComment("INLINED initializations"); |
---|
| 2598 | |
---|
| 2599 | switch( useLang ) { |
---|
| 2600 | case C_LANG: bprintf( InlineCode[ C_INIT ].code ); |
---|
| 2601 | break; |
---|
| 2602 | case F77_LANG: bprintf( InlineCode[ F77_INIT ].code ); |
---|
| 2603 | break; |
---|
| 2604 | case F90_LANG: bprintf( InlineCode[ F90_INIT ].code ); |
---|
| 2605 | break; |
---|
| 2606 | case MATLAB_LANG: bprintf( InlineCode[ MATLAB_INIT ].code ); |
---|
| 2607 | break; |
---|
| 2608 | } |
---|
| 2609 | FlushBuf(); |
---|
| 2610 | |
---|
| 2611 | NewLines(1); |
---|
| 2612 | WriteComment("End INLINED initializations"); |
---|
| 2613 | NewLines(1); |
---|
| 2614 | |
---|
| 2615 | MATLAB_Inline(" VAR = VAR(:);\n FIX = FIX(:);\n" ); |
---|
| 2616 | |
---|
| 2617 | FreeVariable( X ); |
---|
| 2618 | FreeVariable( I ); |
---|
| 2619 | FunctionEnd( INITVAL ); |
---|
| 2620 | FreeVariable( INITVAL ); |
---|
| 2621 | } |
---|
| 2622 | |
---|
| 2623 | |
---|
| 2624 | |
---|
| 2625 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2626 | void GenerateShuffle_user2kpp() |
---|
| 2627 | { |
---|
| 2628 | int i,k,l; |
---|
| 2629 | int Shuffle_user2kpp; |
---|
| 2630 | |
---|
| 2631 | UseFile( utilFile ); |
---|
| 2632 | |
---|
| 2633 | Shuffle_user2kpp = DefFnc( "Shuffle_user2kpp", 2, "function to copy concentrations from USER to KPP"); |
---|
| 2634 | FunctionBegin( Shuffle_user2kpp, V_USER, V ); |
---|
| 2635 | |
---|
| 2636 | k = 0;l = 0; |
---|
| 2637 | for( i = 1; i < SpcNr; i++) { |
---|
| 2638 | if( ReverseCode[i] < 0 ) { |
---|
| 2639 | if( SpeciesTable[i].type == VAR_SPC ) k++; |
---|
| 2640 | continue; |
---|
| 2641 | } |
---|
| 2642 | switch( SpeciesTable[i].type ) { |
---|
| 2643 | case VAR_SPC: |
---|
| 2644 | if( k < initNr ) { |
---|
| 2645 | Assign( Elm( V, ReverseCode[i] ), Elm( V_USER, k++ ) ); |
---|
| 2646 | break; |
---|
| 2647 | } |
---|
| 2648 | case FIX_SPC: |
---|
| 2649 | case DUMMY_SPC: |
---|
| 2650 | default: break; |
---|
| 2651 | } |
---|
| 2652 | } |
---|
| 2653 | |
---|
| 2654 | FunctionEnd( Shuffle_user2kpp ); |
---|
| 2655 | FreeVariable( Shuffle_user2kpp ); |
---|
| 2656 | } |
---|
| 2657 | |
---|
| 2658 | |
---|
| 2659 | |
---|
| 2660 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2661 | void GenerateShuffle_kpp2user() |
---|
| 2662 | { |
---|
| 2663 | int i,k,l; |
---|
| 2664 | int Shuffle_kpp2user; |
---|
| 2665 | |
---|
| 2666 | UseFile( utilFile ); |
---|
| 2667 | |
---|
| 2668 | Shuffle_kpp2user = DefFnc( "Shuffle_kpp2user", 2, "function to restore concentrations from KPP to USER"); |
---|
| 2669 | FunctionBegin( Shuffle_kpp2user, V, V_USER ); |
---|
| 2670 | |
---|
| 2671 | k = 0; l = 0; |
---|
| 2672 | for( i = 0; i < SpcNr; i++) { |
---|
| 2673 | if( ReverseCode[i] < 0 ) { |
---|
| 2674 | if( SpeciesTable[i].type == VAR_SPC ) k++; |
---|
| 2675 | continue; |
---|
| 2676 | } |
---|
| 2677 | switch( SpeciesTable[i].type ) { |
---|
| 2678 | case VAR_SPC: |
---|
| 2679 | if( k < initNr ) |
---|
| 2680 | Assign( Elm( V_USER, k++ ), Elm( V, ReverseCode[i] ) ); |
---|
| 2681 | break; |
---|
| 2682 | case FIX_SPC: |
---|
| 2683 | case DUMMY_SPC: |
---|
| 2684 | default: break; |
---|
| 2685 | } |
---|
| 2686 | } |
---|
| 2687 | |
---|
| 2688 | FunctionEnd( Shuffle_kpp2user ); |
---|
| 2689 | FreeVariable( Shuffle_kpp2user ); |
---|
| 2690 | } |
---|
| 2691 | |
---|
| 2692 | |
---|
| 2693 | |
---|
| 2694 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2695 | void GenerateGetMass() |
---|
| 2696 | { |
---|
| 2697 | int i; |
---|
| 2698 | int atm, spc; |
---|
| 2699 | int GETMASS, MASS; |
---|
| 2700 | SPECIES_DEF *sp; |
---|
| 2701 | int numass; |
---|
| 2702 | |
---|
| 2703 | UseFile( utilFile ); |
---|
| 2704 | |
---|
| 2705 | nmass = 0; |
---|
| 2706 | for( atm = 0; atm < AtomNr; atm++ ) |
---|
| 2707 | if( AtomTable[atm].masscheck ) nmass++; |
---|
| 2708 | if( nmass == 0 ) nmass = 1; |
---|
| 2709 | |
---|
| 2710 | MASS = DefvElm( "Mass", real, nmass, "value of mass balance" ); |
---|
| 2711 | GETMASS = DefFnc( "GetMass", 2, "compute total mass of selected atoms"); |
---|
| 2712 | FunctionBegin( GETMASS, CL, MASS); |
---|
| 2713 | |
---|
| 2714 | numass = 0; |
---|
| 2715 | for( atm = 0; atm < AtomNr; atm++ ) { |
---|
| 2716 | if( AtomTable[atm].masscheck ) { |
---|
| 2717 | sum = Const( 0 ); |
---|
| 2718 | for( spc = 0; spc < SpcNr; spc++ ) { |
---|
| 2719 | sp = &SpeciesTable[ Code[spc] ]; |
---|
| 2720 | for( i = 0; i < sp->nratoms; i++ ) { |
---|
| 2721 | if( sp->atoms[i].code == atm ) { |
---|
| 2722 | sum = Add( sum, Mul( Const( sp->atoms[i].nr ), |
---|
| 2723 | Elm( CL, spc ) ) ); |
---|
| 2724 | } |
---|
| 2725 | } |
---|
| 2726 | } |
---|
| 2727 | Assign( Elm( MASS, numass ), sum ); |
---|
| 2728 | numass++; |
---|
| 2729 | } |
---|
| 2730 | } |
---|
| 2731 | |
---|
| 2732 | FunctionEnd( GETMASS ); |
---|
| 2733 | FreeVariable( GETMASS ); |
---|
| 2734 | } |
---|
| 2735 | |
---|
| 2736 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2737 | void GenerateMakefile() |
---|
| 2738 | { |
---|
| 2739 | char buf[100]; |
---|
| 2740 | |
---|
| 2741 | if ( useLang == MATLAB_LANG ) return; |
---|
| 2742 | |
---|
| 2743 | sprintf( buf, "Makefile_%s", rootFileName ); |
---|
| 2744 | makeFile = fopen(buf, "w"); |
---|
| 2745 | if( makeFile == 0 ) { |
---|
| 2746 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2747 | } |
---|
| 2748 | |
---|
| 2749 | UseFile( makeFile ); |
---|
| 2750 | |
---|
| 2751 | IncludeCode( "%s/util/Makefile", Home ); |
---|
| 2752 | |
---|
| 2753 | } |
---|
| 2754 | |
---|
| 2755 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2756 | void GenerateMex() |
---|
| 2757 | { |
---|
| 2758 | char buf[100], suffix[5]; |
---|
| 2759 | |
---|
| 2760 | if (useLang == MATLAB_LANG) return; |
---|
| 2761 | if (useMex == 0) return; |
---|
| 2762 | |
---|
| 2763 | switch( useLang ) { |
---|
| 2764 | case F77_LANG: sprintf( suffix, "f"); |
---|
| 2765 | break; |
---|
| 2766 | case F90_LANG: sprintf( suffix, "f90"); |
---|
| 2767 | break; |
---|
| 2768 | case C_LANG: sprintf( suffix, "c"); |
---|
| 2769 | break; |
---|
| 2770 | default: printf("\nCannot create mex files for language %d\n", useLang); |
---|
| 2771 | exit(1); |
---|
| 2772 | break; |
---|
| 2773 | } |
---|
| 2774 | |
---|
| 2775 | sprintf( buf, "%s_mex_Fun.%s", rootFileName, suffix ); |
---|
| 2776 | mex_funFile = fopen(buf, "w"); |
---|
| 2777 | if( mex_funFile == 0 ) { |
---|
| 2778 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2779 | } |
---|
| 2780 | UseFile( mex_funFile ); |
---|
| 2781 | IncludeCode( "%s/util/Mex_Fun", Home ); |
---|
| 2782 | |
---|
| 2783 | if (useJacSparse) { |
---|
| 2784 | sprintf( buf, "%s_mex_Jac_SP.%s", rootFileName, suffix ); |
---|
| 2785 | mex_jacFile = fopen(buf, "w"); |
---|
| 2786 | if( mex_jacFile == 0 ) { |
---|
| 2787 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2788 | } |
---|
| 2789 | UseFile( mex_jacFile ); |
---|
| 2790 | IncludeCode( "%s/util/Mex_Jac_SP", Home ); |
---|
| 2791 | } |
---|
| 2792 | |
---|
| 2793 | if (useHessian) { |
---|
| 2794 | sprintf( buf, "%s_mex_Hessian.%s", rootFileName, suffix ); |
---|
| 2795 | mex_hessFile = fopen(buf, "w"); |
---|
| 2796 | if( mex_hessFile == 0 ) { |
---|
| 2797 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2798 | } |
---|
| 2799 | UseFile( mex_hessFile ); |
---|
| 2800 | IncludeCode( "%s/util/Mex_Hessian", Home ); |
---|
| 2801 | } |
---|
| 2802 | |
---|
| 2803 | } |
---|
| 2804 | |
---|
| 2805 | |
---|
| 2806 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2807 | void GenerateMatlabTemplates() |
---|
| 2808 | { |
---|
| 2809 | char buf[200], suffix[5]; |
---|
| 2810 | |
---|
| 2811 | if (useLang != MATLAB_LANG) return; |
---|
| 2812 | |
---|
| 2813 | |
---|
| 2814 | sprintf( buf, "%s_Fun_Chem.m", rootFileName ); |
---|
| 2815 | mex_funFile = fopen(buf, "w"); |
---|
| 2816 | if( mex_funFile == 0 ) { |
---|
| 2817 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2818 | } |
---|
| 2819 | UseFile( mex_funFile ); |
---|
| 2820 | IncludeCode( "%s/util/Template_Fun_Chem", Home ); |
---|
| 2821 | |
---|
| 2822 | sprintf( buf, "%s_Update_SUN.m", rootFileName ); |
---|
| 2823 | mex_funFile = fopen(buf, "w"); |
---|
| 2824 | if( mex_funFile == 0 ) { |
---|
| 2825 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2826 | } |
---|
| 2827 | UseFile( mex_funFile ); |
---|
| 2828 | IncludeCode( "%s/util/UpdateSun", Home ); |
---|
| 2829 | |
---|
| 2830 | if (useJacSparse) { |
---|
| 2831 | sprintf( buf, "%s_Jac_Chem.m", rootFileName ); |
---|
| 2832 | mex_jacFile = fopen(buf, "w"); |
---|
| 2833 | if( mex_jacFile == 0 ) { |
---|
| 2834 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2835 | } |
---|
| 2836 | UseFile( mex_jacFile ); |
---|
| 2837 | IncludeCode( "%s/util/Template_Jac_Chem", Home ); |
---|
| 2838 | } |
---|
| 2839 | |
---|
| 2840 | if (useHessian) { |
---|
| 2841 | } |
---|
| 2842 | |
---|
| 2843 | } |
---|
| 2844 | |
---|
| 2845 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2846 | void GenerateF90Modules(char where) |
---|
| 2847 | { |
---|
| 2848 | char buf[200]; |
---|
| 2849 | |
---|
| 2850 | if (useLang != F90_LANG) return; |
---|
| 2851 | |
---|
| 2852 | switch (where) { |
---|
| 2853 | case 'h': |
---|
| 2854 | |
---|
| 2855 | sprintf( buf, "%s_Precision.f90", rootFileName ); |
---|
| 2856 | sparse_dataFile = fopen(buf, "w"); |
---|
| 2857 | if( sparse_dataFile == 0 ) { |
---|
| 2858 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 2859 | } |
---|
| 2860 | UseFile( sparse_dataFile ); |
---|
| 2861 | F90_Inline("\nMODULE %s_Precision\n", rootFileName ); |
---|
| 2862 | F90_Inline("!"); |
---|
| 2863 | F90_Inline("! Definition of different levels of accuracy"); |
---|
| 2864 | F90_Inline("! for REAL variables using KIND parameterization"); |
---|
| 2865 | F90_Inline("!"); |
---|
| 2866 | F90_Inline("! KPP SP - Single precision kind"); |
---|
| 2867 | F90_Inline(" INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6,30)"); |
---|
| 2868 | F90_Inline("! KPP DP - Double precision kind"); |
---|
| 2869 | F90_Inline(" INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14,300)"); |
---|
| 2870 | F90_Inline("! KPP QP - Quadruple precision kind"); |
---|
| 2871 | F90_Inline(" INTEGER, PARAMETER :: qp = SELECTED_REAL_KIND(18,400)"); |
---|
| 2872 | F90_Inline("\nEND MODULE %s_Precision\n\n", rootFileName ); |
---|
| 2873 | |
---|
| 2874 | UseFile( initFile ); |
---|
| 2875 | F90_Inline("MODULE %s_Initialize\n", rootFileName ); |
---|
| 2876 | F90_Inline(" USE %s_Parameters, ONLY: dp, NVAR, NFIX", rootFileName); |
---|
| 2877 | F90_Inline(" IMPLICIT NONE\n", rootFileName ); |
---|
| 2878 | F90_Inline("CONTAINS\n\n"); |
---|
| 2879 | |
---|
| 2880 | UseFile( param_headerFile ); |
---|
| 2881 | F90_Inline("MODULE %s_Parameters\n", rootFileName ); |
---|
| 2882 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2883 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2884 | |
---|
| 2885 | UseFile( global_dataFile ); |
---|
| 2886 | F90_Inline("MODULE %s_Global\n", rootFileName ); |
---|
| 2887 | if ( useDeclareValues ) |
---|
| 2888 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2889 | else |
---|
| 2890 | F90_Inline(" USE %s_Parameters, ONLY: dp, NSPEC, NVAR, NFIX, NREACT", rootFileName); |
---|
| 2891 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2892 | |
---|
| 2893 | UseFile( functionFile ); |
---|
| 2894 | F90_Inline("MODULE %s_Function\n", rootFileName ); |
---|
| 2895 | if ( useDeclareValues ) |
---|
| 2896 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2897 | else |
---|
| 2898 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2899 | F90_Inline(" IMPLICIT NONE\n", rootFileName ); |
---|
| 2900 | Declare( A ); /* mz_rs_20050117 */ |
---|
| 2901 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2902 | |
---|
| 2903 | UseFile( rateFile ); |
---|
| 2904 | F90_Inline("MODULE %s_Rates\n", rootFileName ); |
---|
| 2905 | if ( useDeclareValues ) |
---|
| 2906 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2907 | else |
---|
| 2908 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2909 | F90_Inline(" USE %s_Global", rootFileName ); |
---|
| 2910 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2911 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2912 | |
---|
| 2913 | if ( useStochastic ) { |
---|
| 2914 | UseFile(stochasticFile); |
---|
| 2915 | F90_Inline("MODULE %s_Stochastic\n", rootFileName); |
---|
| 2916 | if ( useDeclareValues ) |
---|
| 2917 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2918 | else |
---|
| 2919 | F90_Inline(" USE %s_Parameters, ONLY: NVAR, NFIX, NREACT", rootFileName ); |
---|
| 2920 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2921 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2922 | } |
---|
| 2923 | |
---|
| 2924 | if ( useJacSparse ) { |
---|
| 2925 | UseFile(sparse_jacFile); |
---|
| 2926 | F90_Inline("MODULE %s_JacobianSP\n", rootFileName); |
---|
| 2927 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2928 | } |
---|
| 2929 | |
---|
| 2930 | UseFile( jacobianFile ); |
---|
| 2931 | F90_Inline("MODULE %s_Jacobian\n", rootFileName ); |
---|
| 2932 | if ( useDeclareValues ) |
---|
| 2933 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2934 | else |
---|
| 2935 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2936 | if ( useJacSparse ) |
---|
| 2937 | F90_Inline(" USE %s_JacobianSP\n", rootFileName); |
---|
| 2938 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2939 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2940 | |
---|
| 2941 | if ( useStoicmat ) { |
---|
| 2942 | UseFile(sparse_stoicmFile); |
---|
| 2943 | F90_Inline("MODULE %s_StoichiomSP\n", rootFileName); |
---|
| 2944 | F90_Inline(" USE %s_Precision", rootFileName); |
---|
| 2945 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2946 | |
---|
| 2947 | UseFile( stoichiomFile ); |
---|
| 2948 | F90_Inline("MODULE %s_Stoichiom\n", rootFileName); |
---|
| 2949 | if ( useDeclareValues ) |
---|
| 2950 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2951 | else |
---|
| 2952 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2953 | F90_Inline(" USE %s_StoichiomSP\n", rootFileName); |
---|
| 2954 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2955 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2956 | } |
---|
| 2957 | |
---|
| 2958 | if ( useHessian ) { |
---|
| 2959 | UseFile(sparse_hessFile); |
---|
| 2960 | F90_Inline("MODULE %s_HessianSP\n", rootFileName); |
---|
| 2961 | /* F90_Inline(" USE %s_Precision", rootFileName ); */ /* mz_rs_20050321 */ |
---|
| 2962 | F90_Inline(" PUBLIC\n SAVE\n"); |
---|
| 2963 | |
---|
| 2964 | UseFile( hessianFile ); |
---|
| 2965 | F90_Inline("MODULE %s_Hessian\n", rootFileName); |
---|
| 2966 | if ( useDeclareValues ) |
---|
| 2967 | F90_Inline(" USE %s_Precision", rootFileName ); |
---|
| 2968 | else |
---|
| 2969 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2970 | F90_Inline(" USE %s_HessianSP\n", rootFileName); |
---|
| 2971 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2972 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2973 | } |
---|
| 2974 | |
---|
| 2975 | UseFile( monitorFile ); |
---|
| 2976 | F90_Inline("MODULE %s_Monitor", rootFileName); |
---|
| 2977 | |
---|
| 2978 | UseFile( linalgFile ); |
---|
| 2979 | F90_Inline("MODULE %s_LinearAlgebra\n", rootFileName); |
---|
| 2980 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2981 | /* mz_rs_20050511+ if( useJacSparse ) added */ |
---|
| 2982 | if ( useJacSparse ) |
---|
| 2983 | F90_Inline(" USE %s_JacobianSP\n", rootFileName); |
---|
| 2984 | /* mz_rs_20050511- */ |
---|
| 2985 | /* mz_rs_20050321+ */ |
---|
| 2986 | /* if (useHessian) */ |
---|
| 2987 | /* F90_Inline(" USE %s_HessianSP\n", rootFileName); */ |
---|
| 2988 | /* mz_rs_20050321- */ |
---|
| 2989 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2990 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2991 | |
---|
| 2992 | UseFile( utilFile ); |
---|
| 2993 | F90_Inline("MODULE %s_Util\n", rootFileName); |
---|
| 2994 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 2995 | F90_Inline(" IMPLICIT NONE", rootFileName ); |
---|
| 2996 | F90_Inline("\nCONTAINS\n\n"); |
---|
| 2997 | |
---|
| 2998 | /* Here we define the model module which aggregates everything */ |
---|
| 2999 | /* put module rootFileName_Model into separate file */ |
---|
| 3000 | /* (reusing "sparse_dataFile" as done above for _Precision file) */ |
---|
| 3001 | sprintf( buf, "%s_Model.f90", rootFileName ); |
---|
| 3002 | sparse_dataFile = fopen(buf, "w"); |
---|
| 3003 | if( sparse_dataFile == 0 ) { |
---|
| 3004 | FatalError(3,"%s: Can't create file", buf ); |
---|
| 3005 | } |
---|
| 3006 | UseFile( sparse_dataFile ); |
---|
| 3007 | F90_Inline("MODULE %s_Model\n", rootFileName); |
---|
| 3008 | F90_Inline("!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); |
---|
| 3009 | F90_Inline("! Completely defines the model %s", rootFileName); |
---|
| 3010 | F90_Inline("! by using all the associated modules"); |
---|
| 3011 | F90_Inline("!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); |
---|
| 3012 | F90_Inline("\n USE %s_Precision", rootFileName ); |
---|
| 3013 | F90_Inline(" USE %s_Parameters", rootFileName ); |
---|
| 3014 | F90_Inline(" USE %s_Global", rootFileName ); |
---|
| 3015 | F90_Inline(" USE %s_Function", rootFileName ); |
---|
| 3016 | F90_Inline(" USE %s_Integrator", rootFileName ); |
---|
| 3017 | F90_Inline(" USE %s_Rates", rootFileName ); |
---|
| 3018 | if ( useStochastic ) |
---|
| 3019 | F90_Inline(" USE %s_Stochastic", rootFileName ); |
---|
| 3020 | if ( useJacobian ) |
---|
| 3021 | F90_Inline(" USE %s_Jacobian", rootFileName ); |
---|
| 3022 | if ( useHessian ) |
---|
| 3023 | F90_Inline(" USE %s_Hessian", rootFileName); |
---|
| 3024 | if ( useStoicmat ) |
---|
| 3025 | F90_Inline(" USE %s_Stoichiom", rootFileName); |
---|
| 3026 | F90_Inline(" USE %s_LinearAlgebra", rootFileName); |
---|
| 3027 | F90_Inline(" USE %s_Monitor", rootFileName); |
---|
| 3028 | F90_Inline(" USE %s_Util", rootFileName); |
---|
| 3029 | F90_Inline("\nEND MODULE %s_Model\n", rootFileName); |
---|
| 3030 | |
---|
| 3031 | /* mz_rs_20050518+ */ |
---|
| 3032 | /* UseFile( driverFile ); */ |
---|
| 3033 | /* WriteDelim(); */ |
---|
| 3034 | /* mz_rs_20050518- */ |
---|
| 3035 | |
---|
| 3036 | break; |
---|
| 3037 | |
---|
| 3038 | case 't': |
---|
| 3039 | |
---|
| 3040 | /* mz_rs_20050117+ */ |
---|
| 3041 | UseFile( initFile ); |
---|
| 3042 | F90_Inline("\nEND MODULE %s_Initialize\n", rootFileName ); |
---|
| 3043 | /* mz_rs_20050117- */ |
---|
| 3044 | |
---|
| 3045 | UseFile( param_headerFile ); |
---|
| 3046 | F90_Inline("\nEND MODULE %s_Parameters\n", rootFileName ); |
---|
| 3047 | |
---|
| 3048 | UseFile( global_dataFile ); |
---|
| 3049 | F90_Inline("\nEND MODULE %s_Global\n", rootFileName ); |
---|
| 3050 | |
---|
| 3051 | UseFile( functionFile ); |
---|
| 3052 | F90_Inline("\nEND MODULE %s_Function\n", rootFileName ); |
---|
| 3053 | |
---|
| 3054 | UseFile( rateFile ); |
---|
| 3055 | F90_Inline("\nEND MODULE %s_Rates\n", rootFileName ); |
---|
| 3056 | |
---|
| 3057 | if ( useStochastic ) { |
---|
| 3058 | UseFile(stochasticFile); |
---|
| 3059 | F90_Inline("\nEND MODULE %s_Stochastic\n", rootFileName); |
---|
| 3060 | } |
---|
| 3061 | |
---|
| 3062 | if ( useJacSparse ) { |
---|
| 3063 | UseFile(sparse_jacFile); |
---|
| 3064 | F90_Inline("\nEND MODULE %s_JacobianSP\n", rootFileName); |
---|
| 3065 | } |
---|
| 3066 | |
---|
| 3067 | UseFile( jacobianFile ); |
---|
| 3068 | F90_Inline("\nEND MODULE %s_Jacobian\n", rootFileName ); |
---|
| 3069 | |
---|
| 3070 | if ( useStoicmat ) { |
---|
| 3071 | UseFile(sparse_stoicmFile); |
---|
| 3072 | F90_Inline("\nEND MODULE %s_StoichiomSP\n", rootFileName); |
---|
| 3073 | |
---|
| 3074 | UseFile( stoichiomFile ); |
---|
| 3075 | F90_Inline("\nEND MODULE %s_Stoichiom\n", rootFileName); |
---|
| 3076 | } |
---|
| 3077 | |
---|
| 3078 | if ( useHessian ) { |
---|
| 3079 | UseFile(sparse_hessFile); |
---|
| 3080 | F90_Inline("\nEND MODULE %s_HessianSP\n", rootFileName); |
---|
| 3081 | |
---|
| 3082 | UseFile( hessianFile ); |
---|
| 3083 | F90_Inline("\nEND MODULE %s_Hessian\n", rootFileName ); |
---|
| 3084 | } |
---|
| 3085 | |
---|
| 3086 | UseFile(monitorFile); |
---|
| 3087 | F90_Inline("\nEND MODULE %s_Monitor", rootFileName); |
---|
| 3088 | |
---|
| 3089 | UseFile( linalgFile ); |
---|
| 3090 | F90_Inline("\nEND MODULE %s_LinearAlgebra\n", rootFileName); |
---|
| 3091 | |
---|
| 3092 | UseFile( utilFile ); |
---|
| 3093 | F90_Inline("\nEND MODULE %s_Util\n", rootFileName); |
---|
| 3094 | |
---|
| 3095 | break; |
---|
| 3096 | |
---|
| 3097 | default: |
---|
| 3098 | printf("\n Unrecognized option '%s' in GenerateF90Modules\n", where); |
---|
| 3099 | break; |
---|
| 3100 | } |
---|
| 3101 | } |
---|
| 3102 | |
---|
| 3103 | |
---|
| 3104 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 3105 | void Generate() |
---|
| 3106 | { |
---|
| 3107 | int i, j; |
---|
| 3108 | int n; |
---|
| 3109 | |
---|
| 3110 | VarStartNr = 0; |
---|
| 3111 | FixStartNr = VarNr; |
---|
| 3112 | |
---|
| 3113 | real = useDouble ? DOUBLE : REAL; |
---|
| 3114 | |
---|
| 3115 | n = MAX_OUTBUF; |
---|
| 3116 | for( i = 1; i < INLINE_OPT; i++ ) |
---|
| 3117 | if( InlineCode[i].maxlen > n ) |
---|
| 3118 | n = InlineCode[i].maxlen; |
---|
| 3119 | |
---|
| 3120 | outBuf = (char*)malloc( n ); |
---|
| 3121 | outBuffer = outBuf; |
---|
| 3122 | |
---|
| 3123 | switch( useLang ) { |
---|
| 3124 | case F77_LANG: Use_F( rootFileName ); |
---|
| 3125 | break; |
---|
| 3126 | case F90_LANG: Use_F90( rootFileName ); |
---|
| 3127 | break; |
---|
| 3128 | case C_LANG: Use_C( rootFileName ); |
---|
| 3129 | break; |
---|
| 3130 | case MATLAB_LANG: Use_MATLAB( rootFileName ); |
---|
| 3131 | break; |
---|
| 3132 | default: printf("\n Language no '%s' unknown\n",useLang ); |
---|
| 3133 | } |
---|
| 3134 | printf("\nKPP is initializing the code generation."); |
---|
| 3135 | InitGen(); |
---|
| 3136 | |
---|
| 3137 | if ( useLang == F90_LANG ) |
---|
| 3138 | GenerateF90Modules('h'); |
---|
| 3139 | |
---|
| 3140 | GenerateMap(); |
---|
| 3141 | |
---|
| 3142 | /* if( (useLang == F77_LANG)||(useLang == F90_LANG)||(useLang == C_LANG) ) |
---|
| 3143 | {*/ |
---|
| 3144 | printf("\nKPP is generating the monitor data:"); |
---|
| 3145 | printf("\n - %s_Monitor",rootFileName); |
---|
| 3146 | GenerateMonitorData(); |
---|
| 3147 | /* }*/ |
---|
| 3148 | |
---|
| 3149 | printf("\nKPP is generating the utility data:"); |
---|
| 3150 | printf("\n - %s_Util",rootFileName); |
---|
| 3151 | GenerateUtil(); |
---|
| 3152 | |
---|
| 3153 | printf("\nKPP is generating the global declarations:"); |
---|
| 3154 | printf("\n - %s_Main",rootFileName); |
---|
| 3155 | GenerateGData(); |
---|
| 3156 | |
---|
| 3157 | |
---|
| 3158 | printf("\nKPP is generating the ODE function:"); |
---|
| 3159 | printf("\n - %s_Function",rootFileName); |
---|
| 3160 | GenerateFun(); |
---|
| 3161 | |
---|
| 3162 | if ( useStochastic ) { |
---|
| 3163 | printf("\nKPP is generating the Stochastic description:"); |
---|
| 3164 | printf("\n - %s_Function",rootFileName); |
---|
| 3165 | GenerateStochastic(); |
---|
| 3166 | } |
---|
| 3167 | |
---|
| 3168 | if ( useJacobian ) { |
---|
| 3169 | printf("\nKPP is generating the ODE Jacobian:"); |
---|
| 3170 | printf("\n - %s_Jacobian\n - %s_JacobianSP",rootFileName,rootFileName); |
---|
| 3171 | GenerateJacobianSparseData(); |
---|
| 3172 | GenerateJac(); |
---|
| 3173 | if ( (useLang == F77_LANG)||(useLang == F90_LANG)||(useLang == C_LANG) ) { |
---|
| 3174 | GenerateJacVect(); |
---|
| 3175 | GenerateJacTRVect(); |
---|
| 3176 | if( useJacSparse ) { |
---|
| 3177 | printf("\nKPP is generating the linear algebra routines:"); |
---|
| 3178 | printf("\n - %s_LinearAlgebra",rootFileName); |
---|
| 3179 | GenerateSparseUtil(); |
---|
| 3180 | GenerateSolve(); |
---|
| 3181 | GenerateTRSolve(); |
---|
| 3182 | } |
---|
| 3183 | } |
---|
| 3184 | } |
---|
| 3185 | |
---|
| 3186 | GenerateBlas(); |
---|
| 3187 | |
---|
| 3188 | if( useHessian ) { |
---|
| 3189 | printf("\nKPP is generating the Hessian:"); |
---|
| 3190 | printf("\n - %s_Hessian\n - %s_HessianSP",rootFileName,rootFileName); |
---|
| 3191 | GenerateHessian(); |
---|
| 3192 | GenerateHessianSparseData(); |
---|
| 3193 | } |
---|
| 3194 | |
---|
| 3195 | printf("\nKPP is generating the utility functions:"); |
---|
| 3196 | printf("\n - %s_Util",rootFileName); |
---|
| 3197 | |
---|
| 3198 | GenerateInitialize(); |
---|
| 3199 | |
---|
| 3200 | GenerateShuffle_user2kpp(); |
---|
| 3201 | GenerateShuffle_kpp2user(); |
---|
| 3202 | |
---|
| 3203 | printf("\nKPP is generating the rate laws:"); |
---|
| 3204 | printf("\n - %s_Rates",rootFileName); |
---|
| 3205 | |
---|
| 3206 | GenerateRateLaws(); |
---|
| 3207 | GenerateUpdateSun(); |
---|
| 3208 | GenerateUpdateRconst(); |
---|
| 3209 | GenerateUpdatePhoto(); |
---|
| 3210 | GenerateGetMass(); |
---|
| 3211 | |
---|
| 3212 | |
---|
| 3213 | printf("\nKPP is generating the parameters:"); |
---|
| 3214 | printf("\n - %s_Parameters",rootFileName); |
---|
| 3215 | |
---|
| 3216 | GenerateParamHeader(); |
---|
| 3217 | |
---|
| 3218 | printf("\nKPP is generating the global data:"); |
---|
| 3219 | printf("\n - %s_Global",rootFileName); |
---|
| 3220 | |
---|
| 3221 | GenerateGlobalHeader(); |
---|
| 3222 | |
---|
| 3223 | if ( (useLang == F77_LANG)||(useLang == C_LANG)||(useLang == MATLAB_LANG) ) { |
---|
| 3224 | printf("\nKPP is generating the sparsity data:"); |
---|
| 3225 | if( useJacSparse ) { |
---|
| 3226 | GenerateJacobianSparseHeader(); |
---|
| 3227 | printf("\n - %s_JacobianSP",rootFileName); |
---|
| 3228 | } |
---|
| 3229 | if( useHessian ) { |
---|
| 3230 | GenerateHessianSparseHeader(); |
---|
| 3231 | printf("\n - %s_HessianSP",rootFileName); |
---|
| 3232 | } |
---|
| 3233 | } |
---|
| 3234 | |
---|
| 3235 | if ( useStoicmat ) { |
---|
| 3236 | printf("\nKPP is generating the stoichiometric description files:"); |
---|
| 3237 | printf("\n - %s_Stoichiom\n - %s_StoichiomSP",rootFileName,rootFileName); |
---|
| 3238 | GenerateReactantProd(); |
---|
| 3239 | GenerateJacReactantProd(); |
---|
| 3240 | GenerateStoicmSparseData(); |
---|
| 3241 | if ( (useLang == F77_LANG)||(useLang == C_LANG)||(useLang == MATLAB_LANG) ) |
---|
| 3242 | GenerateStoicmSparseHeader(); |
---|
| 3243 | GenerateDFunDRcoeff(); |
---|
| 3244 | GenerateDJacDRcoeff(); |
---|
| 3245 | } |
---|
| 3246 | |
---|
| 3247 | printf("\nKPP is generating the driver from %s.f90:", driver); |
---|
| 3248 | printf("\n - %s_Main",rootFileName); |
---|
| 3249 | |
---|
| 3250 | if ( (useLang == F77_LANG)||(useLang == F90_LANG)||(useLang == C_LANG) ) |
---|
| 3251 | GenerateIntegrator(); |
---|
| 3252 | |
---|
| 3253 | /* mz_rs_20050518+ no driver file if driver = none */ |
---|
| 3254 | if( strcmp( driver, "none" ) != 0 ) |
---|
| 3255 | GenerateDriver(); |
---|
| 3256 | /* mz_rs_20050518- */ |
---|
| 3257 | |
---|
| 3258 | if ( (useLang == F77_LANG)||(useLang == F90_LANG)||(useLang == C_LANG) ) |
---|
| 3259 | GenerateMakefile(); |
---|
| 3260 | |
---|
| 3261 | if ( useLang == F90_LANG ) |
---|
| 3262 | GenerateF90Modules('t'); |
---|
| 3263 | |
---|
| 3264 | if ( useLang == MATLAB_LANG ) |
---|
| 3265 | GenerateMatlabTemplates(); |
---|
| 3266 | |
---|
| 3267 | if ( (useLang == F77_LANG)||(useLang == F90_LANG)||(useLang == C_LANG) ) |
---|
| 3268 | GenerateMex(); |
---|
| 3269 | |
---|
| 3270 | /* mz_rs_20050117+ */ |
---|
| 3271 | if( initFile ) fclose( initFile ); |
---|
| 3272 | /* mz_rs_20050117- */ |
---|
| 3273 | if( driverFile ) fclose( driverFile ); |
---|
| 3274 | if( functionFile ) fclose( functionFile ); |
---|
| 3275 | if( global_dataFile ) fclose( global_dataFile ); |
---|
| 3276 | if( hessianFile ) fclose( hessianFile ); |
---|
| 3277 | if( integratorFile ) fclose( integratorFile ); |
---|
| 3278 | if( jacobianFile ) fclose( jacobianFile ); |
---|
| 3279 | if( linalgFile ) fclose( linalgFile ); |
---|
| 3280 | if( mapFile ) fclose( mapFile ); |
---|
| 3281 | if( makeFile ) fclose( makeFile ); |
---|
| 3282 | if( monitorFile ) fclose( monitorFile ); |
---|
| 3283 | if( mex_funFile ) fclose( mex_funFile ); |
---|
| 3284 | if( mex_jacFile ) fclose( mex_jacFile ); |
---|
| 3285 | if( mex_hessFile ) fclose( mex_hessFile ); |
---|
| 3286 | if( param_headerFile ) fclose( param_headerFile ); |
---|
| 3287 | if( rateFile ) fclose( rateFile ); |
---|
| 3288 | if( sparse_dataFile ) fclose( sparse_dataFile ); |
---|
| 3289 | if( sparse_jacFile ) fclose( sparse_jacFile ); |
---|
| 3290 | if( sparse_hessFile ) fclose( sparse_hessFile ); |
---|
| 3291 | if( sparse_stoicmFile ) fclose( sparse_stoicmFile ); |
---|
| 3292 | if( stoichiomFile ) fclose( stoichiomFile ); |
---|
| 3293 | if( utilFile ) fclose( utilFile ); |
---|
| 3294 | if( stochasticFile ) fclose( stochasticFile ); |
---|
| 3295 | |
---|
| 3296 | } |
---|
| 3297 | |
---|
| 3298 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 3299 | int* AllocIntegerVector(int n, char* message) |
---|
| 3300 | { |
---|
| 3301 | int* vec; |
---|
| 3302 | if ( ( vec=(int*)calloc(n,sizeof(int)) ) == NULL ) |
---|
| 3303 | FatalError(-30,"%s: Cannot allocate vector.",message); |
---|
| 3304 | return vec; |
---|
| 3305 | } |
---|
| 3306 | |
---|
| 3307 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 3308 | /* Allocates a matrix of integers */ |
---|
| 3309 | int** AllocIntegerMatrix(int m, int n, char* message) |
---|
| 3310 | { |
---|
| 3311 | int** mat; |
---|
| 3312 | int i; |
---|
| 3313 | if ( (mat = (int**)calloc(m,sizeof(int*)))==NULL ) { |
---|
| 3314 | FatalError(-30,"%s: Cannot allocate matrix.", message); |
---|
| 3315 | } |
---|
| 3316 | for (i=0; i<m; i++) |
---|
| 3317 | if ( (mat[i] = (int*)calloc(n,sizeof(int)))==NULL ) { |
---|
| 3318 | FatalError(-30,"%s: Cannot allocate matrix[%d].", message, i); |
---|
| 3319 | } |
---|
| 3320 | return mat; |
---|
| 3321 | } |
---|
| 3322 | |
---|
| 3323 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 3324 | /* Frees the memory allocated by AllocIntegerMatrix */ |
---|
| 3325 | void FreeIntegerMatrix(int** mat, int m, int n) |
---|
| 3326 | { |
---|
| 3327 | int i; |
---|
| 3328 | for (i=0; i<m; i++) |
---|
| 3329 | free(mat[i]); |
---|
| 3330 | free(mat); |
---|
| 3331 | } |
---|
| 3332 | |
---|
| 3333 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 3334 | /* i = C-type index; returns language-appropriate index */ |
---|
| 3335 | int Index( int i ) |
---|
| 3336 | { |
---|
| 3337 | switch( useLang ) { |
---|
| 3338 | case C_LANG: |
---|
| 3339 | return i; |
---|
| 3340 | case F77_LANG: |
---|
| 3341 | return i+1; |
---|
| 3342 | case F90_LANG: |
---|
| 3343 | return i+1; |
---|
| 3344 | case MATLAB_LANG: |
---|
| 3345 | return i+1; |
---|
| 3346 | default: printf("\n Unknown language no %d\n",useLang); |
---|
| 3347 | exit(1); |
---|
| 3348 | } |
---|
| 3349 | } |
---|