[2696] | 1 | /****************************************************************************** |
---|
| 2 | |
---|
| 3 | KPP - The Kinetic PreProcessor |
---|
| 4 | Builds simulation code for chemical kinetic systems |
---|
| 5 | |
---|
[4843] | 6 | Copyright (C) -2021 996 Valeriu Damian and Adrian Sandu |
---|
| 7 | Copyright (C) -2021 005 Adrian Sandu |
---|
[2696] | 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,MAX_SPECIES); |
---|
| 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 | } |
---|