source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/src/gen.c @ 4843

Last change on this file since 4843 was 4843, checked in by raasch, 16 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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