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