Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 4725 was 4481, checked in by maronga, 5 years ago

Bugfix for copyright updates in document_changes; copyright update applied to all files
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) -2020 996 Valeriu Damian and Adrian Sandu
7  Copyright (C) -2020 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.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |