chem/KPP/kpp/kpp-2.1/src.org/kpp.c

   1 /******************************************************************************
   2
   3   KPP - The Kinetic PreProcessor
   4         Builds simulation code for chemical kinetic systems
   5
   6   Copyright (C) 1995-1996 Valeriu Damian and Adrian Sandu
   7   Copyright (C) 1997-2005 Adrian Sandu
   8
   9   KPP is free software; you can redistribute it and/or modify it under the
  10   terms of the GNU General Public License as published by the Free Software
  11   Foundation (http://www.gnu.org/copyleft/gpl.html); either version 2 of the
  12   License, or (at your option) any later version.
  13
  14   KPP is distributed in the hope that it will be useful, but WITHOUT ANY
  15   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  17   details.
  18
  19   You should have received a copy of the GNU General Public License along
  20   with this program; if not, consult http://www.gnu.org/copyleft/gpl.html or
  21   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  22   Boston, MA  02111-1307,  USA.
  23
  24   Adrian Sandu
  25   Computer Science Department
  26   Virginia Polytechnic Institute and State University
  27   Blacksburg, VA 24060
  28   E-mail: sandu@cs.vt.edu
  29
  30 ******************************************************************************/
  31
  32
  33 #include <stdlib.h>
  34 #include <string.h>
  35 #include "gdata.h"
  36 #include "scan.h"
  37
  38 char *eqFileName;
  39 char *rootFileName = "ff";
  40 char Home[ MAX_PATH ] = "";
  41
  42 short int linStru[ MAX_SPECIES ];
  43 short int colStru[ MAX_SPECIES ];
  44 short int bestStru[ MAX_SPECIES ];
  45 short int *Stru;
  46
  47 enum stru_criteria { UNSORT, LINSORT, COLSORT, BESTSORT };
  48
  49 void EqCopy( EQ_VECT e1, EQ_VECT e2 )
  50 {
  51 int i;
  52
  53   for( i = 0; i < EqnNr; i++ ) e2[i] = e1[i];
  54 }
  55
  56 int NoSort( const void *p1, const void *p2 )
  57 {
  58   return -1;
  59 }
  60
  61 int CodeCmp( const void *p1, const void *p2 )
  62 {
  63 CODE *c1, *c2;
  64
  65   c1 = (CODE*)p1;
  66   c2 = (CODE*)p2;
  67
  68   if ( *c1 < *c2 ) return -1;
  69   if ( *c1 > *c2 ) return 1;
  70   return 0;
  71 }
  72
  73 int CodeRCmp( const void *p1, const void *p2 )
  74 {
  75 int rc1, rc2;
  76 CODE *c1, *c2;
  77
  78   c1 = (CODE*)p1;
  79   c2 = (CODE*)p2;
  80
  81   rc1 = Reactive[ ReverseCode[ *c1 ] ];
  82   rc2 = Reactive[ ReverseCode[ *c2 ] ];
  83   if ( rc1 > rc2 ) return -1;
  84   if ( rc1 < rc2 ) return 1;
  85   if ( *c1 < *c2 ) return -1;
  86   if ( *c1 > *c2 ) return 1;
  87   return 0;
  88 }
  89
  90 int CodeSSCmp( const void *p1, const void *p2 )
  91 {
  92   return -CodeRCmp(p1,p2);
  93 }
  94
  95 int CodeSCmp( const void *p1, const void *p2 )
  96 {
  97 CODE *c1, *c2;
  98 short int sc1, sc2;
  99
 100   c1 = (CODE*)p1;
 101   c2 = (CODE*)p2;
 102
 103   sc1 = Stru[ ReverseCode[ *c1 ] ];
 104   sc2 = Stru[ ReverseCode[ *c2 ] ];
 105
 106   if ( sc1 > sc2 ) return 1;
 107   if ( sc1 < sc2 ) return -1;
 108   if ( *c1 < *c2 ) return 1;
 109   if ( *c1 > *c2 ) return -1;
 110   return 0;
 111 }
 112
 113 void UpdateStructJ()
 114 {
 115 int i,j,k;
 116
 117   for ( i=0; i<VarNr; i++ )
 118     for ( j=0; j<VarNr; j++ )
 119       structJ[i][j]=(i==j)?1:0;
 120
 121   for (i = 0; i < VarNr; i++)
 122     for (j = 0; j < VarNr; j++)
 123       for (k = 0; k < EqnNr; k++)
 124         if ( Stoich[i][k]*((Stoich_Left[j][k])?1:0) != 0.0 )
 125           structJ[i][j]=1;
 126
 127   for ( i=0; i<VarNr; i++ )
 128     for ( j=0; j<VarNr; j++ )
 129       LUstructJ[i][j]=structJ[i][j];
 130
 131 }
 132
 133 int ComputeLUStructJ()
 134 {
 135 int i,j,k;
 136 int nu,nl;
 137
 138   for (j = 0; j < VarNr-1; j++) {
 139     for (i = j+1; i < VarNr; i++) {
 140       if( LUstructJ[i][j] ) {
 141         for (k = j; k < VarNr; k++)
 142          /*   LUstructJ[i][k] += LUstructJ[j][k]; */
 143            if ( LUstructJ[j][k] != 0 )
 144               LUstructJ[i][k] = 1;
 145       }
 146     }
 147   }
 148
 149
 150   nu = 0; nl = 0;
 151   for (i = 0; i < VarNr; i++)
 152     for (j = 0; j < VarNr; j++)
 153       if( LUstructJ[i][j] ) {
 154         if(i > j) nl++;
 155         if(i <= j) nu++;
 156       }
 157
 158   return nu+nl;
 159 }
 160
 161 int LUnonZero()
 162 {
 163 CODE v[MAX_SPECIES];
 164 CODE *var;
 165 int i,j,k;
 166 int nu,nl;
 167
 168   var = v;
 169   if( Stru != bestStru ) {
 170     for( i=0; i<VarNr; i++ )
 171       var[i] = Code[i];
 172     qsort( (void*)var, VarNr, sizeof(CODE), CodeSCmp );
 173   } else {
 174     var = bestStru;
 175   }
 176
 177   for (i = 0; i < VarNr; i++)
 178     for (j = 0; j < VarNr; j++)
 179       LUstructJ[i][j] = structJ[ ReverseCode[var[i]] ][ ReverseCode[var[j]] ];
 180
 181   return ComputeLUStructJ();
 182 }
 183
 184 void LinColSparsity()
 185 {
 186 int i,j,k;
 187 int nlin, ncol;
 188 FILE * fff;
 189
 190   for ( i=0; i<VarNr; i++ )
 191     for ( j=0; j<VarNr; j++ )
 192       structJ[i][j]=(i==j)?1:0;
 193
 194   for (i = 0; i < VarNr; i++)
 195     for (j = 0; j < VarNr; j++)
 196       for (k = 0; k < EqnNr; k++)
 197         if ( Stoich[i][k]*((Stoich_Left[j][k])?1:0) != 0.0 )
 198           structJ[i][j]=1;
 199
 200   for ( i=0; i<VarNr; i++ ) {
 201     linStru[i] = 0;
 202     for (j = 0; j < VarNr; j++)
 203       linStru[i] += structJ[i][j];
 204   }
 205
 206   for ( i=0; i<VarNr; i++ ) {
 207     colStru[i] = 0;
 208     for (j = 0; j < VarNr; j++)
 209       colStru[i] += structJ[j][i];
 210     colStru[i] *= linStru[i];
 211   }
 212
 213   Stru = linStru;
 214   nlin = LUnonZero();
 215   Stru = colStru;
 216   ncol = LUnonZero();
 217   if( nlin <= ncol ) {
 218     Stru = linStru;
 219     LUnonZero();
 220   }
 221 }
 222
 223 void BestSparsity()
 224 {
 225 int i,j,k;
 226 int cnz, lnz;
 227 int best, crt;
 228 int best_i;
 229 int tmp;
 230 int s;
 231
 232   UpdateStructJ();
 233
 234   for ( i=0; i<VarNr; i++ )
 235     bestStru[i] = Code[i];
 236
 237   for ( s=0; s<VarNr-1; s++ ) {
 238     best = MAX_SPECIES*MAX_SPECIES; best_i = 0;
 239     for ( i=s; i<VarNr; i++ ) {
 240       cnz = 0;lnz = 0;
 241       for (j = s; j < VarNr; j++) {
 242         cnz += (LUstructJ[i][j]?1:0);
 243         lnz += (LUstructJ[j][i]?1:0);
 244       }
 245       crt = (cnz-1)*(lnz-1);
 246       if( crt < best ) {
 247         best = crt;
 248         best_i = i;
 249       }
 250     }
 251     for ( i=0; i<VarNr; i++ ) {
 252       tmp = LUstructJ[s][i];
 253       LUstructJ[s][i] = LUstructJ[best_i][i];
 254       LUstructJ[best_i][i] = tmp;
 255     }
 256     for ( i=0; i<VarNr; i++ ) {
 257       tmp = LUstructJ[i][s];
 258       LUstructJ[i][s] = LUstructJ[i][best_i];
 259       LUstructJ[i][best_i] = tmp;
 260     }
 261     tmp = bestStru[s];
 262     bestStru[s] = bestStru[best_i];
 263     bestStru[best_i] = tmp;
 264
 265     for (i = s+1; i < VarNr; i++) {
 266       if( LUstructJ[i][s] ) {
 267         for (k = s; k < VarNr; k++)
 268           LUstructJ[i][k] += LUstructJ[s][k];
 269       }
 270     }
 271   }
 272
 273   Stru = bestStru;
 274 }
 275
 276 void ReorderSpecies( int criteria )
 277 {
 278 CODE *var;
 279 CODE *fix;
 280 CODE *dummy;
 281 EQ_VECT *tmpStoich_Left;
 282 EQ_VECT *tmpStoich_Right;
 283 EQ_VECT *tmpStoich;
 284 CODE *tmpCode;
 285 CODE *tmpReact;
 286 int i, k;
 287 int new;
 288 int (*cmpVar)(const void *, const void *);
 289 int (*cmpFix)(const void *, const void *);
 290 int dummyNr;
 291
 292   cmpVar = CodeCmp;
 293   cmpFix = CodeCmp;
 294
 295   switch( criteria ) {
 296     case UNSORT:   cmpVar = useJacobian ? CodeRCmp : CodeCmp;
 297                    break;
 298     case LINSORT:  cmpVar = useJacobian ? CodeSCmp : CodeCmp;
 299                    Stru = linStru;
 300                    break;
 301     case COLSORT:  cmpVar = useJacobian ? CodeSCmp : CodeCmp;
 302                    Stru = colStru;
 303                    break;
 304     case BESTSORT: cmpVar = useJacobian ? NoSort : CodeCmp;
 305                    break;
 306   }
 307
 308   VarNr = 0;
 309   VarActiveNr = 0;
 310   FixNr = 0;
 311   dummyNr = 0;
 312
 313   var = (CODE*)malloc( SpcNr * sizeof(CODE) );
 314   fix = (CODE*)malloc( SpcNr * sizeof(CODE) );
 315   dummy = (CODE*)malloc( 5 * sizeof(CODE) );
 316   tmpStoich_Left = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
 317   tmpStoich_Right = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
 318   tmpStoich = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
 319   tmpCode = (CODE*)malloc( SpcNr * sizeof(CODE) );
 320   tmpReact = (CODE*)malloc( SpcNr * sizeof(CODE) );
 321
 322   for( i = 0; i < SpcNr; i++ ) {
 323     switch( SpeciesTable[ Code[i] ].type ) {
 324       case VAR_SPC:  var[ VarNr++ ] = Code[ i ];
 325                      break;
 326       case FIX_SPC:  fix[ FixNr++ ] = Code[ i ];
 327                      break;
 328       case DUMMY_SPC:dummy[ dummyNr++ ] = Code[ i ];
 329                      break;
 330     }
 331   }
 332
 333   if( Stru != bestStru ) {
 334     qsort( (void*)var, VarNr, sizeof(CODE), cmpVar );
 335   } else {
 336     for( i = 0; i < SpcNr; i++ )
 337       var[i] = bestStru[i];
 338   }
 339   qsort( (void*)fix, FixNr, sizeof(CODE), cmpFix );
 340
 341   for( i = 0; i < SpcNr; i++ ) {
 342     EqCopy( Stoich_Left[i], tmpStoich_Left[i] );
 343     EqCopy( Stoich_Right[i], tmpStoich_Right[i] );
 344     EqCopy( Stoich[i], tmpStoich[i] );
 345     tmpCode[i] = Code[i];
 346     tmpReact[i] = Reactive[i];
 347   }
 348
 349   SpcNr -= dummyNr;
 350   dummyNr = 0;
 351
 352   k = 0;
 353   for( i = 0; i < VarNr; i++ ) {
 354     new = ReverseCode[ var[i] ];
 355     EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
 356     EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
 357     EqCopy( tmpStoich[ new ], Stoich[ k ] );
 358     Code[ k ] = tmpCode[ new ];
 359     Reactive[ k ] = tmpReact[ new ];
 360     if( Reactive[ k ] ) VarActiveNr++;
 361     k++;
 362   }
 363   for( i = 0; i < FixNr; i++ ) {
 364     new = ReverseCode[ fix[i] ];
 365     EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
 366     EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
 367     EqCopy( tmpStoich[ new ], Stoich[ k ] );
 368     Code[ k ] = tmpCode[ new ];
 369     Reactive[ k ] = tmpReact[ new ];
 370     k++;
 371   }
 372   for( i = 0; i < dummyNr; i++ ) {
 373     new = ReverseCode[ dummy[i] ];
 374     EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
 375     EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
 376     EqCopy( tmpStoich[ new ], Stoich[ k ] );
 377     Code[ k ] = tmpCode[ new ];
 378     Reactive[ k ] = tmpReact[ new ];
 379     k++;
 380   }
 381
 382
 383   for( i = 0; i < SpcNr+dummyNr; i++ )
 384     ReverseCode[ Code[i] ] = i;
 385
 386   free( tmpReact );
 387   free( tmpCode );
 388   free( tmpStoich );
 389   free( tmpStoich_Right );
 390   free( tmpStoich_Left );
 391   free( dummy );
 392   free( fix );
 393   free( var );
 394
 395   fflush(stdout);
 396 }
 397
 398 /* Allocate Internal Arrays */
 399 void  AllocInternalArrays( void )
 400 {
 401 int i;
 402
 403 if ( (Stoich_Left =(float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL )
 404     FatalError(-30,"Cannot allocate Stoich_Left.\n");
 405
 406 for (i=0; i<MAX_SPECIES; i++)
 407     if ( (Stoich_Left[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
 408         FatalError(-30,"Cannot allocate Stoich_Left[%d]",i);
 409     }
 410
 411 if ( (Stoich_Right = (float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL )
 412     FatalError(-30,"Cannot allocate Stoich_Right.\n");
 413
 414 for (i=0; i<MAX_SPECIES; i++)
 415     if ( (Stoich_Right[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
 416         FatalError(-30,"Cannot allocate Stoich_Right[%d].",i);
 417     }
 418
 419 if ( (Stoich = (float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL )
 420     FatalError(-30,"Cannot allocate Stoich.\n");
 421
 422 for (i=0; i<MAX_SPECIES; i++)
 423     if ( (Stoich[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
 424         FatalError(-30,"Cannot allocate Stoich[%d].",i);
 425     }
 426
 427 }
 428
 429
 430 /* Allocate Structure Arrays */
 431 void  AllocStructArrays( void )
 432 {
 433 int i;
 434
 435
 436 if ( (structB = (int**)calloc(EqnNr,sizeof(int*)))==NULL )
 437     FatalError(-30, "Cannot allocate structB.");
 438
 439 for (i=0; i<EqnNr; i++)
 440     if ( (structB[i] =(int*) calloc(SpcNr,sizeof(int)))==NULL )
 441         FatalError(-30, "Cannot allocate structB[%d].\n",i);
 442
 443 if ( (structJ = (int**)calloc(SpcNr,sizeof(int*)))==NULL )
 444     FatalError(-30, "Cannot allocate structJ.");
 445
 446 for (i=0; i<SpcNr; i++)
 447     if ( (structJ[i] =(int*) calloc(SpcNr,sizeof(int)))==NULL )
 448         FatalError(-30, "Cannot allocate structJ[%d].\n",i);
 449
 450 if ( (LUstructJ = (int**)calloc(SpcNr,sizeof(int*)))==NULL )
 451     FatalError(-30, "Cannot allocate LUstructJ.");
 452
 453 for (i=0; i<SpcNr; i++)
 454     if ( (LUstructJ[i] = (int*)calloc(SpcNr,sizeof(int)))==NULL )
 455         FatalError(-30, "Cannot allocate LUstructJ[%d].\n",i);
 456
 457 }
 458
 459 /*******************************************************************/
 460 int Postprocess( char * root )
 461 {
 462 char buf[ 200 ];
 463 char cmd[500];
 464 char cmdexe[500];
 465 static char tmpfile[] = "kppfile.tmp";
 466 FILE * fp;
 467
 468   if ( useLang == MATLAB_LANG ) {
 469  /*  Add rate function definitions as internal functions to the Update_RCONST file*/
 470     sprintf( buf, "cat %s_Update_RCONST.m %s_Rates.m > tmp; mv tmp %s_Update_RCONST.m;",
 471         root, root, root );
 472     system( buf );
 473   }
 474
 475 /*    Postprocessing to replace parameter names by values in the declarations
 476   strcpy( cmd, "sed " );
 477   sprintf( cmd, "%s -e 's/(NVAR)/(%d)/g'", cmd, VarNr );
 478   sprintf( cmd, "%s -e 's/(NFIX)/(%d)/g'", cmd, FixNr );
 479   sprintf( cmd, "%s -e 's/(NSPEC)/(%d)/g'", cmd,SpcNr );
 480   sprintf( cmd, "%s -e 's/(NREACT)/(%d)/g'", cmd, EqnNr );
 481   sprintf( cmd, "%s -e 's/(NONZERO)/(%d)/g'", cmd, Jac_NZ );
 482   sprintf( cmd, "%s -e 's/(LU_NONZERO)/(%d)/g'", cmd, LU_Jac_NZ );
 483   sprintf( cmd, "%s -e 's/(NHESS)/(%)/g'", cmd, Hess_NZ );
 484
 485   sprintf( buf, "%s_Function", rootFileName );
 486   switch( useLang ) {
 487     case F77_LANG: sprintf( buf, "%s.f", buf );
 488                  break;
 489     case F90_LANG: sprintf( buf, "%s.f90", buf );
 490                  break;
 491     case C_LANG: sprintf( buf, "%s.c", buf );
 492                  break;
 493     case MATLAB_LANG: sprintf( buf, "%s.m", buf );
 494                  break;
 495     default: printf("\n Language '%d' not implemented!\n",useLang);
 496                  exit(1);
 497   }
 498   sprintf( cmdexe, "%s %s > %s; mv %s %s;", cmd, buf, tmpfile, tmpfile, buf );
 499   printf("\n\nCMDEXE='%s'\n",cmdexe);
 500   system( cmdexe );
 501 */
 502 }
 503
 504 /*******************************************************************/
 505 int main( int argc, char * argv[] )
 506 {
 507 int status;
 508 char name[ 200 ];
 509 char *p;
 510 int i,j;
 511
 512   AllocInternalArrays();
 513
 514   p = getenv("KPP_HOME");
 515   if( p ) strcpy( Home, p );
 516
 517   switch( argc ) {
 518     case 3: eqFileName = argv[1];
 519             rootFileName = argv[2];
 520             break;
 521     case 2: eqFileName = argv[1];
 522             strcpy( name, eqFileName );
 523             p = name + strlen(name);
 524             while( p > name ) {
 525               if( *p == '.') {
 526                 *p = '\0';
 527                 break;
 528               }
 529               p--;
 530             }
 531             rootFileName = name;
 532             break;
 533     default: FatalError(1,"\nUsage :"
 534                           "\n        kpp <equations file> [output file]\n");
 535   }
 536
 537   printf("\nThis is KPP-%s.\n", KPP_VERSION);
 538
 539   printf("\nKPP is parsing the equation file.");
 540   status = ParseEquationFile( argv[1] );
 541
 542   if( status ) FatalError(2,"%d errors and %d warnings encountered.",
 543                            nError, nWarning );
 544   /* Allocate some internal data structures */
 545   AllocStructArrays();
 546
 547   printf("\nKPP is computing Jacobian sparsity structure.");
 548   ReorderSpecies( UNSORT );
 549   if (useReorder==1){
 550     BestSparsity();
 551     ReorderSpecies( BESTSORT );
 552     }
 553   UpdateStructJ();
 554   ComputeLUStructJ();
 555
 556   if( initNr == -1 ) initNr = VarNr;
 557
 558
 559   printf("\nKPP is starting the code generation.");
 560   Generate( rootFileName );
 561
 562   printf("\nKPP is starting the code post-processing.");
 563   Postprocess( rootFileName );
 564
 565   printf("\n\nKPP has succesfully created the model \"%s\".\n\n",rootFileName);
 566
 567   if( nError ) exit(4);
 568   if( nWarning ) exit(5);
 569
 570   exit(0);
 571 }