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