123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626 |
- /* Utility definitions for code generated
- * from Modelica translator.
- *
- * Author: Hilding Elmqvist, Dynasim AB
- *
- * Copyright (C) 1998-2001 Dynasim AB.
- * All rights reserved.
- *
- * Version: 1.0, 1998-08-18
- * 1.1, 1998-10-21, handling of states
- * 1.2, 1998-10-29, handling of events
- * 1.3, 1999-01-13, matrix features
- * 1.4, 1999-02-22, support for solving nonlinear equations
- */
- #if !defined(RT)
- #include <assert.h>
- #endif
- #ifndef DYMOLA_MOUTIL_C
- #include <math.h>
- #define DYMOLA_MOUTIL_C
- #include <math.h>
- #if !defined(DYMOLA_DSPACE) && !defined(NO_FILE)
- #include "amat.h"
- #endif
- #if defined(DymolaStoreProtectedVariables_)
- static int DYNStoreProtectedVariables=1;
- #else
- static int DYNStoreProtectedVariables=0;
- #endif
- #if !defined(DSE_STRUCT)
- #if defined(DS_EMBEDDED)
- #define DSE_STRUCT sts->
- #else
- #define DSE_STRUCT
- #endif
- #endif
- #ifndef DYNX
- #if _MSC_VER>=1700
- #define DYNX(s,i) (*(s+i))
- #else
- #define DYNX(s,i) s[i]
- #endif
- #endif
- #define Modelica
- DYMOLA_STATIC double* delayID;
- DYMOLA_STATIC double Buffersize;
- #include "dsblock.h"
- #ifndef DYNassert_static
- #define DYNassert_static(e) do { enum { DYN_assert_static_=1/(e) }; } while(0)
- #endif
- struct NonAlias_ {int tag;const char*n;const char*d;double val;double mi;double ma;double w;int ix;};
- #undef SpuriousEvents_
- #define SpuriousEvents_ 0
- /* ---------------------------------------------------- */
- /* Macros to define variable names. */
- #define SParameter(i) ((const char*)(DYNX(cuser_,i)))
- #define Parameter(i) DYNX(DP_,i)
- #define Variable(i) DYNX(W_,i)
- #define State(i) DYNX(X_,i)
- #define Derivative(i) DYNX(F_,i)
- #define Input(i) DYNX(U_,i)
- #define Output(i) DYNX(Y_,i)
- #define Aux(i) DYNX(Aux_,i)
- #define SAux(i) DYNX(DYNAuxStr_,i+DYNStrInit(did_))
- /* ---------------------------------------------------- */
- /* Macros for equations */
- #define IF (
- #define THEN ) ? (
- #define ELSE ) :
- #define AND &&
- #define OR ||
- #define NSAND &
- #define NSOR |
- #define NOT !
- #define EventEPS 1E-8
- #define Crossing(x, y, nr, op) \
- (Qenable_[nr]=true, \
- QZ_[2*nr] = (x) - (y) - EventEPS, \
- QZ_[2*nr+1] = (x) - (y) + EventEPS
- static
- #if defined(_MSC_VER) && _MSC_VER>=1200
- __inline
- #elif __GNUC__
- __inline
- #endif
- double ScaledCross(double x,double y,double nom,double*r) {
- double v=x-y;
- *r=(x>nom) ? x : (-x>nom)? -x : nom;
- return v;
- }
- #define Greater(x, sx, y, sy, nr) HandleEvent2( (ScaledCross((x),(y),1.0,Qscaled_+nr)),sx,">",sy,(nr+1), 0, 0)
- #define Less(x, sx, y, sy, nr) HandleEvent2( (ScaledCross((x),(y),1.0,Qscaled_+nr)) ,sx, "<", sy,(nr+1), 1, 0)
- #define GreaterEqual(x,sx,y,sy,nr) HandleEvent2( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, ">=", sy,(nr+1), 1, 1)
- #define LessEqual(x,sx,y,sy,nr) HandleEvent2( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, "<=", sy,(nr+1), 0, 1)
- #define GreaterMinor(x, sx, y, sy, nr) HandleEvent2Minor( (ScaledCross((x),(y),1.0,Qscaled_+nr)),sx,">",sy,(nr+1), 0, 0)
- #define LessMinor(x, sx, y, sy, nr) HandleEvent2Minor( (ScaledCross((x),(y),1.0,Qscaled_+nr)) ,sx, "<", sy,(nr+1), 1, 0)
- #define GreaterEqualMinor(x,sx,y,sy,nr) HandleEvent2Minor( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, ">=", sy,(nr+1), 1, 1)
- #define LessEqualMinor(x,sx,y,sy,nr) HandleEvent2Minor( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, "<=", sy,(nr+1), 0, 1)
- #define GreaterS(x, sx, y, sy, nr) HandleEvent2S( (ScaledCross((x),(y),1.0,Qscaled_+nr)) ,sx,">",sy,(nr+1), 0, 0)
- #define LessS(x, sx, y, sy, nr) HandleEvent2S( (ScaledCross((x),(y),1.0,Qscaled_+nr)) ,sx, "<", sy,(nr+1), 1, 0)
- #define GreaterEqualS(x,sx,y,sy,nr) HandleEvent2S( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, ">=", sy,(nr+1), 1, 1)
- #define LessEqualS(x,sx,y,sy,nr) HandleEvent2S( (ScaledCross((x),(y),1.0,Qscaled_+nr)), sx, "<=", sy,(nr+1), 0, 1)
- #define absEvent(x,sx,nr) (Qscaled_[nr]=1,HandleEvent2((x),sx," > ","0",(nr+1),0,0)?QRel_[nr]:-QRel_[nr])
- #define signEvent(x,sx,nr1,nr2) (Qscaled_[nr1]=Qscaled_[nr2]=1,HandleEvent2((x),sx," > ","0",(nr1+1),0,0)?1: \
- HandleEvent2(QRel_[nr1],sx," < ","0",(nr2+1),1,0)?-1:0)
- #define simplesignEvent(x,sx,nr1) (Qscaled_[nr1]=1,HandleEvent2((x),sx," > ","0",(nr1+1),0,0)?1:-1)
- #define Equal(x, y, nr) \
- ((x) == (y))
- #define NotEqual(x, y, nr) \
- ((x) != (y))
- #define pre(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;AnyEvent_=true;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- AnyEvent_=true;\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- #define preCont(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else {\
- DYNX(QPre_,nr)=x;\
- }
- #define preW(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- #define preD(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;AnyEvent_=true;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- AnyEvent_=true;\
- AnyDEvent_=true;\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- #define preContD(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else {\
- if (DYNX(QPre_,nr)!=x) AnyDEvent_=true;\
- DYNX(QPre_,nr)=x;\
- }
- #define preWD(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- AnyDEvent_=true;\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- #define preI(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;AnyEvent_=true;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- AnyEvent_=true;\
- AnyIEvent_=true;\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- #define preContI(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else {\
- if (DYNX(QPre_,nr)!=x) AnyIEvent_=true;\
- DYNX(QPre_,nr)=x;\
- }
- #define preWI(x, sx, start, nr)\
- if (Init_ && (DymolaOneIteration_ == 0 || !final_ && DymolaOneIteration_ == 2) ) { DYNX(QPre_,nr)=start;}\
- else if (DYNX(QPre_,nr)!=x && !(!Init_ && DymolaOneIteration_ == 5 && !final_)) {\
- AnyIEvent_=true;\
- if (PrintEvent&(1<<1)) VariableChanged(sx,x,DYNTime);\
- if (!Init_ && final_==1 && (DymolaOneIteration_==2 || DymolaOneIteration_==3 || DymolaOneIteration_==5)) ; else DYNX(QPre_,nr)=x;\
- }
- static double spliceFunction2(double x,double deltax) {
- if (deltax<0) {
- x=-x;
- deltax=-deltax;
- }
- if (x<=-deltax) {
- return 0;
- } else if (x>=deltax) {
- return 1;
- } else {
- double piHalf=asin(1.0);
- return (1+tanh(tan((x/deltax)*piHalf)))/2;
- }
- }
- static double spliceFunction2der(double x,double deltax,double dx,double ddeltax) {
- if (deltax<0) {
- x=-x;
- deltax=-deltax;
- dx=-dx;
- ddeltax=-ddeltax;
- }
- if (x<=-deltax || x>=deltax) {
- return 0;
- } else {
- double piHalf=asin(1.0);
- double x1=(x/deltax);
- double y=x1*piHalf;
- double z=cosh(tan(y))*cos(y);
- return ((dx-x1*ddeltax)*(piHalf)/(2*deltax))/(z*z);
- }
- }
- #define spliceFunction(pos, neg, x, deltax, id) \
- (((double*)DYNhelp)[id]=spliceFunction2((x),(deltax)), \
- ((double*)DYNhelp)[id]*(((double*)DYNhelp)[id]>0 ? (pos) : 0.0)+(1-((double*)DYNhelp)[id])*(((double*)DYNhelp)[id]<1 ? (neg) : 0.0))
- #define spliceFunction_der(pos, neg, x, deltax, pos_der, neg_der, x_der, deltax_der, id) \
- (((double*)DYNhelp)[id]=spliceFunction2((x),(deltax)), \
- ((double*)DYNhelp)[id]*(((double*)DYNhelp)[id]>0 ? (pos_der) : 0.0)+(1-((double*)DYNhelp)[id])*(((double*)DYNhelp)[id]<1 ? (neg_der) : 0.0)\
- + (((double*)DYNhelp)[id]>0 && ((double*)DYNhelp)[id]<1 ? spliceFunction2der((x),(deltax),(x_der),(deltax_der))*((pos)-(neg)): 0.0))
- #define PRE(x, nr) DYNX(QPre_,nr)
- #define Remember_(x, nr) \
- x=(solverHandleEq_ && EqRemember1Time_>-1e30 && EqRemember1Time_<EqRemember2Time_) ?\
- DYNX(EqRemember1_,nr)+Dymola_min(4,(DYNTime-EqRemember1Time_)/(EqRemember2Time_-EqRemember1Time_))*(DYNX(EqRemember2_,nr)-DYNX(EqRemember1_,nr)):\
- DYNX(EqRemember2_,nr)
- #define RememberSimple_(x, nr) \
- DYNX(EqRemember2_,nr)
- #define DoRemember_(x,val,nr) DYNX(EqRemember2_,nr)=(initializationPhase_==1) ? (val) : (x)
- #define event(x) \
- (x)
- #define noEvent(x) \
- (x)
- #define reinit(x, expr) \
- (x = expr, AnyREvent_=1, 0);
- #define reinitNew(x, expr, hv, iv) \
- (((DymolaOneIteration_==0) ? (x = expr):(hv=expr,iv=1)), AnyREvent_=1, 0);
- #define reinitNewArray(x, expr, hv, iv) \
- (((DymolaOneIteration_==0) ? (RealAssign(x,expr)):(RealAssign(hv,expr),iv=1)), AnyREvent_=1, 0);
- #define delay(expr, DelayTime, MaxDelayTime, nr, str) \
- (delayID[nr]= (delayID[nr]==0) ? 1+delayIni(Dymola_abs(DelayTime), Buffersize, DYNTime, expr) : delayID[nr], \
- delayFunction(delayID[nr]-1, Dymola_abs(DelayTime), MaxDelayTime, DYNTime, expr, str))
- #define delayDC(expr, DelayTime, nr, nrC) \
- delayDCFunction(expr, DelayTime, delayID+nr, Buffersize, DYNTime, DYNEvent, QRel_+nrC, Qp_+nrC, Qn_+nrC,\
- Qenable_+1+nrC, QZ_+2*nrC,&AnyEvent_)
- #define delayDD(expr, DelayTime, nr, nrT) \
- delayDDFunction(expr, DelayTime, delayID+nr, Buffersize, DYNTime, DYNEvent, QTimed_+nrT)
- /* ---------------------------------------------------- */
- #define NonLinearSystemOfEquations(Jac, res, x, n, numericJacobian, nr, base, hv, nhv, iv, niv) \
- { \
- DeclareDidRealMatrix1(Jac, n, n, hv); \
- Real *Jac##__DWork=(&hv)+n*n; \
- int Jac##__LDWork=(20+((n+15)*n)/2);\
- int *Jac##__IWork=&iv;\
- static const int QBase=base;\
- struct MyJmpBuf saved_jmp_buf_env; \
- DeclareDidRealVector1(res, n, (&hv)[n*n+(20+((n+15)*n)/2)+0*n]); \
- DeclareDidRealVector1(x, n, (&hv)[n*n+(20+((n+15)*n)/2)+1*n]); \
- DeclareDidRealVector1(xInitial, n, (&hv)[n*n+(20+((n+15)*n)/2)+2*n]); \
- DeclareDidRealVector1(xJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+3*n]); \
- DeclareDidRealVector1(resJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+4*n]); \
- DeclareDidRealVector1(resJac2, n, (&hv)[n*n+(20+((n+15)*n)/2)+5*n]); \
- double lambda_=-2;\
- int *storedXP_=(&iv)+20+n;\
- int firstTime=1;\
- int* storedResetCounterP=(&iv)+21+n;\
- int veryFirstTime=ResetCounter_!=*(storedResetCounterP);\
- int veryFirstTime2=0;\
- int changePmatrix_=0;\
- int changePmatrixOld_=0;\
- int changePloop_=0;\
- int QNnl;int Qinfo,QNLnr,QiOpt,QInfRev;\
- DeclareDidRealMatrix2(Jac, n, n, hv);\
- DeclareDidRealVector2(res, n, (&hv)[n*n+(20+((n+15)*n)/2)+0*n]); \
- DeclareDidRealVector2(x, n, (&hv)[n*n+(20+((n+15)*n)/2)+1*n]); \
- DeclareDidRealVector2(xInitial, n, (&hv)[n*n+(20+((n+15)*n)/2)+2*n]); \
- DeclareDidRealVector2(xJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+3*n]); \
- DeclareDidRealVector2(resJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+4*n]); \
- DeclareDidRealVector2(resJac2, n, (&hv)[n*n+(20+((n+15)*n)/2)+5*n]); \
- DYNassert_static(nhv==n*n+(20+((n+15)*n)/2)+6*n && niv==22+n);\
- if (veryFirstTime) {int i;(*storedXP_)=0;for(i=0;i<n;i++) x##__internal[i]=0,res##__internal[i]=0;for(i=0;i<20+n;++i) Jac##__IWork[i]=0;for(i=0;i<20+((n+15)*n)/2;++i) Jac##__DWork[i]=0;}\
- again##nr: firstTime=1; againsecond##nr:\
- if (DymolaUserHomotopy) {lambda_=1-DymolaHomotopyLambda;if (DymolaHomotopyLambda==0) veryFirstTime=1;}\
- veryFirstTime2=(*(storedResetCounterP)==0);\
- *(storedResetCounterP)=ResetCounter_;\
- saved_jmp_buf_env = DYNGlobal_jmp_buf_env;\
- QNnl = n; \
- if (DYNEvent || ! DAEsolver_) { \
- QNLnr = nr; \
- QiOpt = numericJacobian+1; \
- QInfRev = -11+(numericJacobian?1:-DymolaNonLinearIterations_); \
- if ((!DymolaNewJacobian_ && !veryFirstTime) || !firstTime) QInfRev=-20-DymolaNonLinearIterations_-10*(1-DymolaAllowNewJacobianAnyway_);\
- if (inJacobian_ && !(*storedXP_)) {RealAssign(xJac,x); RealAssign(resJac,res); (*storedXP_)=1;\
- } else if ((*storedXP_)) {RealAssign(x,xJac);(*storedXP_)=inJacobian_;}\
- } else { \
- int i_;\
- QInfRev = 0;\
- for (i_=1; i_<=QNnl; i_++) \
- SetVector(x, i_, ((double*)(DSE_STRUCT X_))[i_+QBase]); \
- }
- #define NonLinearSystemOfEquationsNH(Jac, res, x, n, numericJacobian, nr, base, nin, hv, nhv, iv, niv) \
- { \
- DeclareDidRealMatrix1(Jac, n, n, hv); \
- Real *Jac##__DWork=(&hv)+n*n; \
- int Jac##__LDWork=(20+((n+15)*n)/2);\
- int *Jac##__IWork=&iv;\
- static const int QBase=base;\
- struct MyJmpBuf saved_jmp_buf_env; \
- DeclareDidRealVector1(res, n, (&hv)[n*n+(20+((n+15)*n)/2)+0*n]); \
- DeclareDidRealVector1(x, n, (&hv)[n*n+(20+((n+15)*n)/2)+1*n]); \
- DeclareDidRealVector1(xInitial, n, (&hv)[n*n+(20+((n+15)*n)/2)+2*n]); \
- DeclareDidRealVector1(xJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+3*n]); \
- DeclareDidRealVector1(resJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+4*n]); \
- DeclareDidRealVector1(resJac2, n, (&hv)[n*n+(20+((n+15)*n)/2)+5*n]); \
- Real *hom1=(&hv)+n*n+(20+((n+15)*n)/2)+6*n;\
- Real *hom2=(&hv)+n*n+(20+((n+15)*n)/2)+6*n+nin;\
- Real homlocal,homlocaldelta,homlocalfail;\
- double lambda_=-2;\
- int *storedXP_=(&iv)+20+n;\
- int firstTime=1;\
- int* storedResetCounterP=(&iv)+21+n;\
- int veryFirstTime=ResetCounter_!=*(storedResetCounterP);\
- int veryFirstTime2=0;\
- int changePmatrix_=0;\
- int changePmatrixOld_=0;\
- int changePloop_=0;\
- int QNnl;int Qinfo,QNLnr,QiOpt,QInfRev;\
- DeclareDidRealMatrix2(Jac, n, n, hv);\
- DeclareDidRealVector2(res, n, (&hv)[n*n+(20+((n+15)*n)/2)+0*n]); \
- DeclareDidRealVector2(x, n, (&hv)[n*n+(20+((n+15)*n)/2)+1*n]); \
- DeclareDidRealVector2(xInitial, n, (&hv)[n*n+(20+((n+15)*n)/2)+2*n]); \
- DeclareDidRealVector2(xJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+3*n]); \
- DeclareDidRealVector2(resJac, n, (&hv)[n*n+(20+((n+15)*n)/2)+4*n]); \
- DeclareDidRealVector2(resJac2, n, (&hv)[n*n+(20+((n+15)*n)/2)+5*n]); \
- DYNassert_static(nhv==n*n+(20+((n+15)*n)/2)+2*nin+6*n && niv==22+n);\
- homlocal=-1;homlocaldelta=0;homlocalfail=0;\
- if (veryFirstTime) {int i;(*storedXP_)=0;for(i=0;i<n;i++) x##__internal[i]=0,res##__internal[i]=0;for(i=0;i<20+n;++i) Jac##__IWork[i]=0;for(i=0;i<20+((n+15)*n)/2;++i) Jac##__DWork[i]=0;}\
- again##nr: firstTime=1; againsecond##nr:\
- if (DymolaUserHomotopy) {lambda_=1-DymolaHomotopyLambda;if (DymolaHomotopyLambda==0) veryFirstTime=1;}\
- veryFirstTime2=(*(storedResetCounterP)==0);\
- *(storedResetCounterP)=ResetCounter_;\
- saved_jmp_buf_env = DYNGlobal_jmp_buf_env;\
- QNnl = n; \
- if (DYNEvent || ! DAEsolver_) { \
- QNLnr = nr; \
- QiOpt = numericJacobian+1; \
- QInfRev = -11+(numericJacobian?1:-DymolaNonLinearIterations_); \
- if ((!DymolaNewJacobian_ && !veryFirstTime) || !firstTime) QInfRev=-20-DymolaNonLinearIterations_-10*(1-DymolaAllowNewJacobianAnyway_);\
- if (inJacobian_ && !(*storedXP_)) {RealAssign(xJac,x); RealAssign(resJac,res); (*storedXP_)=1;\
- } else if ((*storedXP_)) {RealAssign(x,xJac);(*storedXP_)=inJacobian_;}\
- } else { \
- int i_;\
- QInfRev = 0;\
- for (i_=1; i_<=QNnl; i_++) \
- SetVector(x, i_, ((double*)X_)[i_+QBase]); \
- }
- #define DummyNonLinear(Jac,residue,m,n,Pvar,var,stateP, hv, nhv, iv, niv) \
- {\
- LIBDS_API int CheckP(const double*,double*,int,int,int*,const double*);\
- LIBDS_API void ComputeP(const double*,double*,int,int,int,int*,double*);\
- DeclareDidRealMatrix1(JacCopy_,n,m,hv);\
- int i_,j_;\
- int*Pivots_=&iv;\
- DeclareDidRealMatrix2(JacCopy_,n,m,hv);\
- DYNassert_static(nhv==n*m && niv==3*n);\
- if (Init_ || (DYNEvent && changePmatrixOld_)) {\
- int info;\
- int uno_=1,n_=n;\
- RealArray varP_=var;\
- ComputeP(&(GetMatrix(Jac,1,1)),&(GetMatrix(JacCopy_,1,1)),m,m,n,Pivots_,Pvar);\
- {int m1_=n-m;double d1_=1.0,d2_=0.0;\
- dgemv_("T",&n_,&m1_,&d1_,Pvar,&n_,varP_.data,&uno_,\
- &d2_,stateP,&uno_);}\
- for(i_=m+1;i_<=n;i_++) {double d_=-(stateP)[i_-(m+1)];\
- for(j_=1;j_<=n;j_++) d_+=(Pvar)[(i_-(m+1))*n+(j_-1)]*varP_.data[j_-1];\
- SetVector(residue,i_,d_);}\
- for(i_=m+1;i_<=n;i_++) for(j_=1;j_<=n;j_++) SetMatrix(Jac,i_,j_,(Pvar)[(i_-(m+1))*n+(j_-1)]);\
- Release();\
- changePmatrixOld_=0;\
- if (PrintEvent&(1<<13)) ReportDummySelection(varnames_,Pvar,DYNTime,m,n);\
- AnyREvent_=true;\
- } else if (CheckP(&(GetMatrix(Jac,1,1)),&(GetMatrix(JacCopy_,1,1)),m,n,Pivots_,Pvar)) changePmatrix_=1;\
- if ((PrintEvent&(1<<14))&&terminal()) ReportDummySelection(varnames_,Pvar,DYNTime,m,n);\
- }
- #define SetInit(x, i, j, var, val) \
- if (firstTime) \
- var = val; \
- if (firstTime && !(DYNEvent || ! DAEsolver_)) \
- SetMatrix(x, i, j, var); \
- else \
- var = GetMatrix(x, i, j)
- #define SetInitVector(x, i, var, val) \
- if (DymolaUserHomotopy && DymolaHomotopyLambda>0) {} else\
- if (firstTime) \
- var = val; \
- if ((DAEsolver_ && !DYNEvent) || !firstTime || inJacobian_) \
- var = GetVector(x, i); \
- else \
- SetVector(x, i, var); \
- #define SetInitVectorNH(x, i, var, val) \
- if (homlocal>0 || DymolaUserHomotopy && DymolaHomotopyLambda>0) var=GetVector(x, i); else\
- if (firstTime) \
- var = val; \
- if ((DAEsolver_ && !DYNEvent) || !firstTime || inJacobian_) \
- var = GetVector(x, i); \
- else \
- SetVector(x, i, var); \
- #define Residues \
- firstTime=0;\
- while (true) {\
- if (setjmp(DYNGlobal_jmp_buf_env.TheBuffer)!=0) {\
- PopAllMarks();\
- GlobalError_ = 0;\
- if (QInfRev>0) {\
- QInfRev+=100;\
- if (!(PrintEvent&(1<<10))) DymosimMessage("Non-linear solver will attempt to handle this problem.\n");\
- }\
- else {\
- if (!(PrintEvent&(1<<10))) DymosimMessage("First evaluation failed for non-linear solver.\n");\
- if (solverHandleEq_) EqRemember2Time_=-1e40;\
- if (!((*(DSE_STRUCT QiErr)||MixedFailFlag_)&&(QInfRev==0))) QInfRev=-1;\
- }\
- } else {\
- #define Jacobian(Jac) \
- if (QInfRev == 2 || QInfRev<=-11 && QInfRev>-20 || (DAEsolver_ && !DYNEvent && !inJacobian_)) {
- #if 1
- LIBDS_API void dymnon7_(int*qinfrev, int qiopt, int qnnl, double*qsol,double* qres,double* qjac,
- double qtol, double qtoldesired,
- int*qinfo, double*qd, int lqd, int*idum, int printevent, int inJacobian_, int qnlnr, double time,
- int*qnlfunc, int*qnljac, int qnlmax, const char*const*varnames,const double*nom,int*ierr___,const char*tag);
- #define DymNonlinear2(QInfRev,QiOpt,QNnl,QSol,QRes,QJac,Qtol,Qtol2,Qinfo,QD,LQD,DumI, \
- PrintEvent,QNLnr,Time,QNLfunc,QNLjac,QNLmax,ierr_,tag) {int i1_,i2_;\
- i1_ = LQD; \
- i2_ = QNLmax; \
- if (QNLnr > QNLmax) { \
- DymosimMessage("Internal error: The capacity for solving nonlinear systems of equations, QNLmax, needs to be increased"); \
- DYN_GOTOLEAVE; \
- } \
- dymnon7_(&QInfRev,QiOpt,QNnl,QSol,QRes,QJac,Qtol,Qtol2,&Qinfo,QD,i1_,DumI,\
- PrintEvent,inJacobian_,QNLnr,Time,QNLfunc,QNLjac,i2_,varnames_,nominal_,ierr_,tag);}
- #else
- LIBDS_API void dymnon6_(int*qinfrev, int qiopt, int qnnl, double*qsol,double* qres,double* qjac,
- double qtol, double qtoldesired,
- int*qinfo, double*qd, int lqd, int*idum, int printevent, int inJacobian_, int qnlnr, double time,
- int*qnlfunc, int*qnljac, int qnlmax, const char*const*varnames,const double*nom,int*ierr___);
- #define DymNonlinear2(QInfRev,QiOpt,QNnl,QSol,QRes,QJac,Qtol,Qtol2,Qinfo,QD,LQD,DumI, \
- PrintEvent,QNLnr,Time,QNLfunc,QNLjac,QNLmax,ierr_,tag) {int i1_,i2_;\
- i1_ = LQD; \
- i2_ = QNLmax; \
- if (QNLnr > QNLmax) { \
- DymosimMessage("Internal error: The capacity for solving nonlinear systems of equations, QNLmax, needs to be increased"); \
- DYN_GOTOLEAVE; \
- } \
- dymnon6_(&QInfRev,QiOpt,QNnl,QSol,QRes,QJac,Qtol,Qtol2,&Qinfo,QD,i1_,DumI,\
- PrintEvent,inJacobian_,QNLnr,Time,QNLfunc,QNLjac,i2_,varnames_,nominal_,ierr_);}
- #endif
- /*
- \
- if (QInfRev==0 && Qinfo==4) goto leave;
- */
- #define OverdeterminedNonLinear(Jac,Jac2,residue,residue2,m,n,m1,m2,var,hv,nhv) \
- {\
- LIBDS_API void CreateJacFull(const double*,double*,int,int,int,const int*);\
- LIBDS_API void ComputePfull(const double*,double*,int,int,int,int,int,int*,int*,double*,const double*,const double*);\
- LIBDS_API int CheckPfull(const double*,double*,int,int,int,int,int,int*,const int*);\
- DeclareDidRealMatrix1(JacCopy_,m,n,hv);\
- int i_,j_;\
- DeclareDidRealMatrix2(JacCopy_,m,n,hv);\
- DYNassert_static(nhv==n*m);\
- if (Init_ || (DYNEvent && changePmatrixOld_)) {\
- ComputePfull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(JacCopy_,1,1)),m,m,n,m1,m2,IWork_+(m+n),IWork_,0,0,0);\
- changePmatrixOld_=0;\
- } else if (CheckPfull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(JacCopy_,1,1)),m,m,n,m1,m2,IWork_+(m+n),IWork_)) changePmatrix_=1;\
- CreateJacFull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(Jac,1,1)),m,m,n,IWork_);\
- } \
- } {\
- LIBDS_API void NonLinearRoutine(int,int,double*,const int*,const double*,const double*,const double*,const double*);\
- NonLinearRoutine(m,n,&(GetVector(residue,1)),IWork_,&(GetVector(residue2,1)),0,0,0);\
- #define OverdeterminedDummyNonLinear(Jac,Jac2,residue,residue2,l,m,n,m1,m2,Pvar,var,stateP,hv,nhv) \
- {\
- LIBDS_API void CreateJacFull(const double*,double*,int,int,int,const int*);\
- LIBDS_API void ComputePfull(const double*,double*,int,int,int,int,int,int*,int*,double*,const double*,const double*);\
- LIBDS_API int CheckPfull(const double*,double*,int,int,int,int,int,int*,const int*);\
- DeclareDidRealMatrix1(JacCopy_,m,n,hv);\
- int i_,j_;\
- DeclareDidRealMatrix2(JacCopy_,m,n,hv);\
- DYNassert_static(nhv==n*m);\
- if (Init_ || (DYNEvent && changePmatrixOld_)) {\
- RealArray varP_=var;\
- ComputePfull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(JacCopy_,1,1)),l,m,n,m1,m2,IWork_+(l+n),IWork_,Pvar,varP_.data,stateP);\
- Release();\
- changePmatrixOld_=0;\
- if (PrintEvent&(1<<13)) ReportDummySelection(varnames_,Pvar,DYNTime,m,n);\
- AnyREvent_=true;\
- } else if (CheckPfull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(JacCopy_,1,1)),l,m,n,m1,m2,IWork_+(l+n),IWork_)) changePmatrix_=1;\
- CreateJacFull(&(GetMatrix(Jac2,1,1)),&(GetMatrix(Jac,1,1)),l,m,n,IWork_);\
- } \
- } {\
- LIBDS_API void NonLinearRoutine(int,int,double*,const int*,const double*,const double*,const double*,const double*);\
- RealArray varP_=var;\
- NonLinearRoutine(m,n,&(GetVector(residue,1)),IWork_,&(GetVector(residue2,1)),Pvar,varP_.data,stateP);\
- Release();
- #define SolveNonLinearSystemOfEquationsInner(Jac, res, x, tag) \
- } \
- } DYNGlobal_jmp_buf_env=saved_jmp_buf_env;\
- if (!(QInfRev != 0 || ((DYNEvent && !Init_ && changePmatrix_)?(changePmatrixOld_=changePmatrix_,QInfRev=-1,changePmatrix_=0,changePloop_++<5):0)))\
- break;\
- if (DYNEvent || ! DAEsolver_) { \
- double Qtol1=Qtol;\
- double Qtol2=DymolaAimForHighAccuracy_ ? Dymola_min(EPS_*0.5,0.01*Qtol) : Qtol;\
- if (inJacobian_ && (*storedXP_)) RealSubtractAssign(res,res,resJac);\
- DymNonlinear2(QInfRev,QiOpt,QNnl,&(GetVector(x,1)),&(GetVector(res,1)),\
- &(GetMatrix(Jac,1,1)),Qtol1,Qtol2,Qinfo, \
- Jac##__DWork,Jac##__LDWork,Jac##__IWork,PrintEvent,QNLnr,DYNTime,QNLfunc,QNLjac,QNLmax,QiErr,tag); \
- if (*QiErr != 0) {*QiErr=0;QInfRev=0;DymosimError("Inner failed");}
- #define SolveNonLinearSystemOfEquations(Jac, res, x, tag) \
- } \
- } DYNGlobal_jmp_buf_env=saved_jmp_buf_env;\
- if (QInfRev == 0) {\
- if (DYNEvent && !Init_ && changePmatrix_) {\
- changePmatrixOld_=changePmatrix_;\
- QInfRev=-1;\
- changePmatrix_=0;\
- if (changePloop_++>=5) break;\
- } else break;\
- }\
- if (DYNEvent || ! DAEsolver_) { \
- double Qtol1=Qtol;\
- double Qtol2=DymolaAimForHighAccuracy_ ? Dymola_min(EPS_*0.5,0.01*Qtol) : Qtol;\
- if (inJacobian_ && (*storedXP_)) RealSubtractAssign(res,res,resJac);\
- DymNonlinear2(QInfRev,QiOpt,QNnl,&(GetVector(x,1)),&(GetVector(res,1)),\
- &(GetMatrix(Jac,1,1)),Qtol1,Qtol2,Qinfo, \
- Jac##__DWork,Jac##__LDWork,Jac##__IWork,PrintEvent,QNLnr,DYNTime,QNLfunc,QNLjac,QNLmax,QiErr,tag); \
- if (*QiErr != 0) DYN_GOTOLEAVE
- #define SolveNonLinearSystemOfEquationsNH(Jac, res, x, nr, tag) \
- } \
- } DYNGlobal_jmp_buf_env=saved_jmp_buf_env;\
- if (QInfRev == 0) {\
- if (DYNEvent && !Init_ && changePmatrix_) {\
- changePmatrixOld_=changePmatrix_;\
- QInfRev=-1;\
- changePmatrix_=0;\
- if (changePloop_++>=5) break;\
- } else break;\
- }\
- if (DYNEvent || ! DAEsolver_) { \
- double Qtol1=Qtol;\
- double Qtol2=DymolaAimForHighAccuracy_ ? Dymola_min(EPS_*0.5,0.01*Qtol) : Qtol;\
- if (inJacobian_ && (*storedXP_)) RealSubtractAssign(res,res,resJac);\
- DymNonlinear2(QInfRev,QiOpt,QNnl,&(GetVector(x,1)),&(GetVector(res,1)),\
- &(GetMatrix(Jac,1,1)),Qtol1,Qtol2,Qinfo, \
- Jac##__DWork,Jac##__LDWork,Jac##__IWork,PrintEvent,QNLnr,DYNTime,QNLfunc,QNLjac,QNLmax,QiErr,tag); \
- if (*QiErr && DYNEvent && homlocal==-1) { *QiErr=0;DymosimMessage("Trying homotopy");homlocal=0;homlocaldelta=0.125;goto again##nr;}\
- if (*QiErr && DYNEvent && homlocal>0 && homlocaldelta>1e-4) {*QiErr=0;DymosimMessage("Homotopy reducing");homlocal-=homlocaldelta;homlocalfail=homlocal;homlocaldelta/=2;homlocal+=homlocaldelta;goto again##nr;}\
- if (*QiErr != 0) DYN_GOTOLEAVE
- #define NonLinearSystemSave(a, nr) \
- if (homlocal==-1) {if (DYNEvent) hom1[nr]=a; else hom2[nr]=a;}\
- else a=hom1[nr]*homlocal+hom2[nr]*(1-homlocal);
- #define SolveNonLinearSystemOfEquationsInit(Jac, res, x, tag) \
- } \
- } DYNGlobal_jmp_buf_env=saved_jmp_buf_env;\
- if (!(QInfRev != 0 || ((DYNEvent && !Init_ && changePmatrix_)?(changePmatrixOld_=changePmatrix_,QInfRev=-1,changePmatrix_=0,changePloop_++<5):0)))\
- break;\
- if (DYNEvent || ! DAEsolver_) { \
- double Qtol1=Qtol;\
- double Qtol2=DymolaAimForHighAccuracy_ ? Dymola_min(EPS_*0.5,0.01*Qtol) : Qtol;\
- if (lambda_<0 && !inJacobian_) {Qtol1=Qtol2;}\
- if (lambda_==-2 && !inJacobian_ && (QInfRev==1 || QInfRev<-1)) {RealAssign(resJac,res);lambda_=-1;}\
- if (lambda_>0 && !DymolaUserHomotopy) {RealSubtractAssign(res,res,RealScale(resJac,lambda_));PopAllMarks();}\
- if (inJacobian_ && (*storedXP_)) RealSubtractAssign(res,res,resJac);\
- DymNonlinear2(QInfRev,QiOpt,QNnl,&(GetVector(x,1)),&(GetVector(res,1)),\
- &(GetMatrix(Jac,1,1)),Qtol1,Qtol2,Qinfo, \
- Jac##__DWork,Jac##__LDWork,Jac##__IWork,PrintEvent,QNLnr,DYNTime,QNLfunc,QNLjac,QNLmax,(DSE_STRUCT QiErr),tag); \
- if (*(DSE_STRUCT QiErr) != 0 && !(*(DSE_STRUCT QiErr)>=1 && *(DSE_STRUCT QiErr)<=4 && lambda_==-1)) goto leave
- #define SolveNonLinearSystemOfEquationsMixed(Jac, res, x, tag) \
- } \
- } DYNGlobal_jmp_buf_env=saved_jmp_buf_env;\
- if (!(QInfRev != 0 || ((DYNEvent && !Init_ && changePmatrix_)?(changePmatrixOld_=changePmatrix_,QInfRev=-1,changePmatrix_=0,changePloop_++<5):0)))\
- break;\
- if (DYNEvent || ! DAEsolver_) { \
- double Qtol1=Qtol;\
- double Qtol2=DymolaAimForHighAccuracy_ ? Dymola_min(EPS_*0.5,0.01*Qtol) : Qtol;\
- if (inJacobian_ && (*storedXP_)) RealSubtractAssign(res,res,resJac);\
- DymNonlinear2(QInfRev,QiOpt,QNnl,&(GetVector(x,1)),&(GetVector(res,1)),\
- &(GetMatrix(Jac,1,1)),Qtol1,Qtol2,Qinfo, \
- Jac##__DWork,Jac##__LDWork,Jac##__IWork,PrintEvent,QNLnr,DYNTime,QNLfunc,QNLjac,QNLmax,QiErr,tag); \
- if (*QiErr != 0) {MixedFailFlag_=*QiErr;*QiErr=0;QInfRev=0;}
- #define VerifyOverdeterminedResidual(res1, n1, res2, n2) \
- if (QInfRev == 0 && !(DYNEvent && !Init_ && changePmatrix_)) {\
- if (!inJacobian_) {\
- int i;double res2=0,res1=0;\
- for(i=1;i<=n1;++i) res1+=sqr(GetVector(res1, i));for(i=1;i<=n2;++i) res2+=sqr(GetVector(res2, i));\
- if (res2>res1*1000+Qtol*Qtol*10) {*QiErr=1;if (!(PrintEvent&(1<<10))) DymosimMessage("Too large residual for overdetermined system");DYN_GOTOLEAVE;}\
- }\
- }
- #define EndNonLinearSystemOfEquations(res, x) \
- } \
- triggerStepEvent_=triggerStepEvent_|| (changePmatrix_ && !DYNEvent);\
- } \
- if (DAEsolver_) {\
- int i_;\
- triggerStepEvent_=triggerStepEvent_|| (changePmatrix_ && !DYNEvent);\
- if (!DYNEvent) for (i_=1; i_<=QNnl; i_++) \
- ((double*)F_)[i_+QBase] = GetVector(res,i_); \
- else for (i_=1; i_<=QNnl; i_++) \
- ((double*)X_)[i_+QBase] = GetVector(x,i_); \
- }\
- }
- #define EndNonLinearSystemOfEquationsNH(res, x, nr) \
- } \
- triggerStepEvent_=triggerStepEvent_|| (changePmatrix_ && !DYNEvent);\
- } \
- if (DAEsolver_) {\
- int i_;\
- triggerStepEvent_=triggerStepEvent_|| (changePmatrix_ && !DYNEvent);\
- if (!DYNEvent) for (i_=1; i_<=QNnl; i_++) \
- ((double*)F_)[i_+QBase] = GetVector(res,i_); \
- else for (i_=1; i_<=QNnl; i_++) \
- ((double*)X_)[i_+QBase] = GetVector(x,i_); \
- }\
- if (homlocal!=-1 && homlocal!=1) {if (homlocaldelta<0.05 && homlocal>homlocalfail+4*homlocalfail) homlocaldelta*=2;homlocal+=homlocaldelta;if (homlocal>1) homlocal=1;goto again##nr;}\
- }
- #define EndNonLinearSystemOfEquationsInit(res, x,nr) \
- } \
- triggerStepEvent_=triggerStepEvent_|| (changePmatrix_ && !DYNEvent);\
- }\
- if (*(DSE_STRUCT QiErr)>=1 && *(DSE_STRUCT QiErr)<=4 && lambda_==-1) {*(DSE_STRUCT QiErr)=0;lambda_=1;if (!(PrintEvent&(1<<10))) DymosimMessage("Trying to solve non-linear system using global homotopy-method.");if (DymolaHomotopyUsed){DymolaUserHomotopy=1;GlobalError_=-997;goto leave;} goto again##nr;} \
- if (lambda_>0 && !DymolaHomotopyUsed) {lambda_-=0.1; if (lambda_-1e-3<0) lambda_=0;goto again##nr;}\
- }
- #define ModelicaPrintMatrixRows(Mat, nRows) \
- DymosimMessageDoubleMatrix(#Mat, &(GetMatrix(Mat,1,1)), nRows, \
- MatrixNrColumns(Mat), MatrixNrRows(Mat) )
- #define MixedSystemOfEquations(n,hv)\
- {\
- DeclareDidRealVector1(oldVector,n, hv); \
- int changeflag = 1;\
- int niter = 0;\
- double deltaTheta = 1e-3;\
- double deltaTheta2 = 1e-12;\
- DeclareDidRealVector2(oldVector,n, hv); \
- while (1) {\
- if (MixedFailFlag_&&DYNEvent&&!changeflag) {AnyEvent_=1;AnyDEvent_=1;MixedFailFlag_=0;AnyEventParameter_=1;break;}\
- if (!changeflag) break;\
- if (MixedFailFlag_&&(niter>0) && !(PrintEvent&(1<<10))) DymosimMessage("Fixed point iterating to handle problem.");\
- changeflag = 0;MixedFailFlag_ = 0;\
- ++niter;\
- if (niter > 13) {\
- if (PrintEvent&(1<<1)) DymosimMessageDouble("Fix point iteration did not converge at time :", DYNTime);\
- if (!AnyEvent_) MixedFailFlag_=1;\
- break;\
- }
- #define MixedInitialSystemOfEquations(n, hv)\
- {\
- DeclareDidRealVector1(oldVector,n, hv); \
- int changeflag = 1;\
- int niter = 0;\
- double deltaTheta = 1e-3;\
- double deltaTheta2 = 1e-12;\
- DeclareDidRealVector2(oldVector,n, hv); \
- while (changeflag) {\
- ++niter;\
- if (niter > 13) {\
- MixedFailFlag_=true;\
- if ((PrintEvent&(1<<1))||MixedFailFlag_&&!(PrintEvent&(1<<10))) DymosimMessageDouble("Fix point iteration did not converge at time :", DYNTime);\
- break;\
- } \
- if (MixedFailFlag_&& DymolaHomotopyUsed && !DymolaUserHomotopy) {if (!(PrintEvent&(1<<10))) DymosimMessage("Solving using global homotopy-method.");MixedFailFlag_=0;DymolaUserHomotopy=1;GlobalError_=-997;goto leave;}\
- if (MixedFailFlag_&&(niter>0) && !(PrintEvent&(1<<10))) DymosimMessage("Fixed point iterating to handle problem.");\
- changeflag = 0;MixedFailFlag_ = 0;
- #define EndMixedSystemOfEquations\
- }}}if (MixedFailFlag_) {*(DSE_STRUCT QiErr)=MixedFailFlag_;DYN_GOTOLEAVE;}}
- #define MixedModeStartBoolean \
- { int mixedModeFlag=(DYNEvent || !RootFinder_) && (niter>1);\
- int thetaIter=-1;int thetaIterFlag=1;double mixedModeTheta,mixedModeTheta2;\
- int thetaNr=0;\
- for(;thetaIterFlag!=0;thetaIter++) { \
- mixedModeTheta=mixedModeTheta2=1; \
- if (thetaIter>-1) {
- #define MixedModeInit(nr, hv) \
- } \
- if (mixedModeFlag) { \
- double *Qold_=&(hv);
- #define ThetaMixedCross(j,nr) \
- if (thetaIter==-1) DYNX(Qold_,j)=DYNX(Qenable_,nr+1) ? DYNX(QRel_,nr) : 0; \
- else if (DYNX(Qold_,j)*DYNX(QRel_,nr)<0) {\
- double th=DYNX(Qold_,j)/(DYNX(Qold_,j)-DYNX(QRel_,nr));\
- if (th<mixedModeTheta) {if(th<mixedModeTheta-deltaTheta2) mixedModeTheta2=mixedModeTheta;\
- mixedModeTheta=th;thetaNr=nr;}\
- else if (th<mixedModeTheta2-deltaTheta2 && th>mixedModeTheta+deltaTheta2) {mixedModeTheta2=th;}\
- }
- #define MixedModeEndBoolean \
- } \
- if (thetaIter==-1) continue; \
- if (thetaIterFlag==2) break; else thetaIterFlag=2;\
- if (mixedModeTheta<1-deltaTheta) {\
- changeflag=1;\
- mixedModeTheta=Dymola_min(mixedModeTheta+deltaTheta,(mixedModeTheta+mixedModeTheta2)/2);\
- if (PrintEvent&(1<<1)) DymosimMessageDouble("Shorter attempt for mixed equation system",mixedModeTheta);\
- } else {\
- mixedModeTheta=1;\
- thetaIterFlag=0;\
- }
- #define UpdateVariable(var, val)\
- {\
- int val1 = val;\
- if (!inJacobian_ && var != val1) {changeflag = 1;var = val1;if (PrintEvent&(1<<1)) VariableChanged(#var,val,DYNTime);}\
- }
- #define UpdateVariableNamed(var, varn, val)\
- {\
- int val1 = val;\
- if (!inJacobian_ && var != val1) {changeflag = 1;var = val1;if (PrintEvent&(1<<1)) VariableChanged(varn,val,DYNTime);}\
- }
- #define UpdateReal(var, nr) \
- if(mixedModeTheta<1) \
- var=GetVector(oldVector,nr)+(mixedModeTheta)*(var-GetVector(oldVector,nr));\
- SetVector(oldVector, nr, var);
- #if defined(DS_EMBEDDED)
- #define TranslatedEquations\
- dse_GetVariables(vars, sts);\
- sectioncondition = true;\
- {
- #define EndTranslatedEquations\
- PopAllMarks();\
- dse_SetVariables(vars, sts);\
- }
- #define ModelFunctionEnd() \
- return 0;finish: return 1;leave:return 2
- #else
- #define TranslatedEquations\
- sectioncondition = true;\
- {
- #define EndTranslatedEquations\
- PopAllMarks();\
- }
- #endif
- #define DefaultSection\
- }\
- sectioncondition = true;\
- {\
- #define InitialSection2\
- }\
- sectioncondition = initializationPhase_?initializationPhase_==2:(initial() && !continueSimulate_);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialSection3\
- }\
- sectioncondition = (initializationPhase_==1) || (initial() && !continueSimulate_);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialSection\
- }\
- sectioncondition = (initializationPhase_ || initialAlsoContinue()) && (DymolaOneIteration_<3);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialBoundSection\
- }\
- sectioncondition = (initializationPhase_ || initialAlsoContinue()) && (DymolaOneIteration_<3) && (!DymolaUserHomotopy || DymolaUserHomotopy && DymolaHomotopyLambda==0);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialStartSection2\
- }\
- sectioncondition = (initializationPhase_?initializationPhase_==2:(initial() && !continueSimulate_)) && (!DymolaUserHomotopy || DymolaHomotopyLambda<=0) && (DymolaOneIteration_<3);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialStartSection3\
- }\
- sectioncondition = ((initializationPhase_==1) || (initial() && !continueSimulate_)) && (!DymolaUserHomotopy || DymolaHomotopyLambda<=0) && (DymolaOneIteration_<3);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define InitialStartSection\
- }\
- sectioncondition = (initializationPhase_ || initialAlsoContinue()) && (DymolaOneIteration_<3) && (!DymolaUserHomotopy || DymolaHomotopyLambda<=0);\
- if (sectioncondition) {\
- *icall_ = 0;
- #define EndInitialSection\
- if (DymolaUserHomotopy) {if (DymolaHomotopyLambda<1) {if (GlobalError_==0)GlobalError_=-996;goto leave;} else DymolaUserHomotopy=0;DymolaHomotopyLambda=1;}
- #define BoundParameterSection \
- }\
- sectioncondition = NewParameters_ && (!DymolaUserHomotopy || DymolaUserHomotopy && DymolaHomotopyLambda==0);\
- if (sectioncondition) {
- DYMOLA_STATIC void initializeData(double *time, double* X_, double* XD_, double* U_,
- double* DP_, int IP_[], Dymola_bool LP_[],
- double* F_, double* Y_, double* W_, double QZ_[],
- double duser_[], int iuser_[], void*cuser_[],
- int final_,int Init,int m_Event,int solverCall);
- DYMOLA_STATIC void initializeDataNew(double *time, double* X_, double* XD_, double* U_,
- double* DP_, int IP_[], Dymola_bool LP_[],
- double* F_, double* Y_, double* W_, double QZ_[],
- double duser_[], int iuser_[], void*cuser_[],struct DYNInstanceData*,
- int final_,int solverCall);
- #define InitializeData(fin) {\
- initializeDataNew(time, X_, XD_, U_, \
- DP_, IP_, LP_, \
- F_, Y_, W_, QZ_, \
- duser_, iuser_, cuser_, did_,\
- fin,fin&&(did_->DymolaEventOptional_var&&solverHandleEq_&&(*idemand_==4))); \
- }
- #define InitializeDataSolver() {\
- initializeDataNew(time, X_, XD_, U_, \
- DP_, IP_, LP_, \
- F_, Y_, W_, QZ_, \
- duser_, iuser_, cuser_, did_,\
- 0,1); \
- }
- #define StartDataBlock DYMOLA_STATIC void initializeData(double *time, double* X_, double* XD_, double* U_, \
- double* DP_, int IP_[], Dymola_bool LP_[], \
- double* F_, double* Y_, double* W_, double QZ_[], \
- double duser_[], int iuser_[], void*cuser_[],\
- int final_,int Init,int m_Event,int solverCall) {tempData.Init_var=Init;tempData.Event_var=m_Event;\
- initializeDataNew(time,X_,XD_,U_,DP_,IP_,LP_,F_,Y_,W_,QZ_,duser_,iuser_,cuser_,&tempData,final_,solverCall);}\
- DYMOLA_STATIC void initializeDataNew(double *time, double* X_, double* XD_, double* U_, \
- double* DP_, int IP_[], Dymola_bool LP_[], \
- double* F_, double* Y_, double* W_, double QZ_[], \
- double duser_[], int iuser_[], void*cuser_[], struct DYNInstanceData*did_,\
- int final_,int solverCall) {
- #define StartEqBlock if (solverCall?1:(Init_||(DYNEvent&&final_)||(!solverHandleEq_&&!inJacobian_&&final_))) {\
- if ((solverCall||(DYNEvent&&solverHandleEq_)) && *time>EqRemember2Time_) {int ii_;EqRemember1Time_=EqRemember2Time_;\
- for(ii_=0;ii_<SizeEq_;ii_++) EqRemember1_[ii_]=EqRemember2_[ii_];}\
- EqRemember2Time_=*time;
- DYMOLA_STATIC void UpdateInitVars(double *time, double* X_, double* XD_, double* U_,
- double* DP_, int IP_[], Dymola_bool LP_[],
- double* F_, double* Y_, double* W_, double QZ_[], double duser_[], int iuser_[], void*cuser_[],struct DYNInstanceData*did_);
- #define OutputSection\
- }\
- sectioncondition = ((*idemand_ >= 1 && *icall_ < 1) || (DYNEvent || *idemand_==4)) && (DymolaOneIteration_!=-1);\
- if (sectioncondition) {\
- *icall_ = 1;
- #define DynamicsSection\
- }\
- sectioncondition = ((*idemand_ >= 2 && *icall_ < 2) || (DYNEvent || *idemand_==4)) && (DymolaOneIteration_!=-1);\
- if (sectioncondition) {\
- *icall_ = 2;
- #define AcceptedSection\
- }\
- sectioncondition = ((*idemand_ >= 3 && *icall_ < 3) || (DYNEvent || *idemand_==4)) && (DymolaOneIteration_!=-1);\
- if (sectioncondition) {\
- *icall_ = 3;
- #define AcceptedSection1\
- }\
- sectioncondition = ((*idemand_ >= 3 && *icall_ < 3) || (DYNEvent || *idemand_==4)) && (DymolaOneIteration_!=-1);\
- if (sectioncondition) {\
- *icall_ = 3;
- #define ODEJacobianSection\
- }\
- sectioncondition = (did_->QJacobian_var || did_->QJacobianSparse_var);\
- if (sectioncondition) {\
- if (did_->QJacobian_var) {int i_;for(i_=0;i_<did_->QJacobianN_var*did_->QJacobianN_var;++i_) did_->QJacobian_var[i_]=0.0;}
- #define ODEBCDJacobianSection\
- sectioncondition = (did_->QBJacobian_var && did_->QCJacobian_var && did_->QDJacobian_var) || (did_->QJacobianSparse_var && did_->QSparseABCD_var);\
- if (sectioncondition && (did_->QBJacobian_var && did_->QCJacobian_var && did_->QDJacobian_var)) {int i_;\
- for(i_=0;i_<did_->QJacobianN_var*did_->QJacobianNU_var;++i_) did_->QBJacobian_var[i_]=0.0;\
- for(i_=0;i_<did_->QJacobianNY_var*did_->QJacobianN_var;++i_) did_->QCJacobian_var[i_]=0.0;\
- for(i_=0;i_<did_->QJacobianNY_var*did_->QJacobianNU_var;++i_) did_->QDJacobian_var[i_]=0.0;\
- }\
- if (sectioncondition)
- #define ODEJacobian(iz, r, c, v) (did_->QJacobian_var?did_->QJacobian_var+(r-1)+(c-1)*did_->QJacobianN_var:did_->QJacobianSparse_var+iz)[0] = v; if (did_->QJacobianSparse_var) {did_->QJacobianSparseR_var[iz]=r-1;did_->QJacobianSparseC_var[iz]=c-1;}
- #define ODEBJacobian(iz, r, c, v) (did_->QBJacobian_var?did_->QBJacobian_var+(r-1)+(c-1)*did_->QJacobianN_var:did_->QJacobianSparse_var+iz)[0]=v; if (did_->QJacobianSparse_var && did_->QSparseABCD_var) {did_->QJacobianSparseR_var[iz]=r-1;did_->QJacobianSparseC_var[iz]=c-1+AnalyticJacobianSize_;}
- #define ODECJacobian(iz, r, c, v) (did_->QCJacobian_var?did_->QCJacobian_var+(r-1)+(c-1)*did_->QJacobianNY_var:did_->QJacobianSparse_var+iz)[0]=v; if (did_->QJacobianSparse_var && did_->QSparseABCD_var) {did_->QJacobianSparseR_var[iz]=r-1+AnalyticJacobianSize_;did_->QJacobianSparseC_var[iz]=c-1;}
- #define ODEDJacobian(iz, r, c, v) (did_->QDJacobian_var?did_->QDJacobian_var+(r-1)+(c-1)*did_->QJacobianNY_var:did_->QJacobianSparse_var+iz)[0]=v; if (did_->QJacobianSparse_var && did_->QSparseABCD_var) {did_->QJacobianSparseR_var[iz]=r-1+AnalyticJacobianSize_;did_->QJacobianSparseC_var[iz]=c-1+AnalyticJacobianSize_;}
- /* ---------------------------------------------------- */
- /* Macros to build variable name table. */
- struct DeclarePhase;
- #if defined(Matlab5) || defined(Matlab51) || defined(SimStruct) || defined(DYMOLA_DSPACE) || defined(DYM2DS) || defined(FMU_SOURCE_CODE_EXPORT)
- #if !defined(DynSimStruct)
- #define DynSimStruct 1
- #endif
- #endif
- #define PreNonAliasDef(nr) DYMOLA_STATIC void declare__##nr(double x0_[],double dp_[],double du_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase);
- #define FuncStartNonAlias(nr) DYMOLA_STATIC void declare__##nr(double x0_[],double dp_[],double du_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase) {
- #if defined(DynSimStruct) || !defined(DYMOSIM)
- /* Not in Dymosim => cannot use blo-routines */
- #if GenerateResultInNonDymosim==2 || (GenerateResultInNonDymosim==1 && !(defined(RT) || defined(DYMOLA_DSPACE)))
- #include "GenerateResultInNonDymosim.h"
- /* Special code for generating result files from Matlab
- #define GenerateResultInMatlab 2
- - do this always
- #define GenerateResultInMatlab 1
- - only when running inside Matlab and not on realtime platforms
- */
- #define DeclareSParameter(par, descpar, ipar, val, mi, ma, w, typ, con) {0+(typ<<3)+(con<<5),par,val,0,0.0,0.0,w,ipar},
- #define DeclareParameter(par, descpar, ipar, val, mi, ma, w,typ,con) {1+(typ<<3)+(con<<5),par,descpar,val,mi, ma,w,ipar},
- #define DeclareState(x, descx, ix, val, mi, ma,w, typ,con) {2+(typ<<3)+(con<<5),x,descx,val,mi, ma, w, ix},
- #define DeclareDerivative(dx, descdx, val, mi, ma,w, typ,con) {3+(typ<<3)+(con<<5),dx,descdx,val,mi, ma, w},
- #define DeclareOutput(var, descvar, iy, val, mi, ma, w,typ, con) {4+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareInput(var, descvar, iu, val, mi, ma, w,typ, con) {5+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareVariable(var,descvar, val, mi, ma, w, typ, con) {6+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareAlias(name, description, aliasname, sign, sigtype, sigindex) {7+(sigtype<<3),name,description,sign,0.0,0.0,0,sigindex},
- #define DeclareAlias2(name, description, aliasname, sign, sigtype, sigindex, typ) {7+(sigtype<<3)+(typ<<6),name,description,sign,0.0,0.0,0,sigindex},
- #define PreNonAliasNew(nr) {declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define PreNonAlias(nr) {PreNonAliasDef(nr);declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define StartNonAlias(nr) static struct NonAlias_ nonalias_##nr []={
- #define EndNonAlias(nr) {100}}; FuncStartNonAlias(nr) \
- struct NonAlias_*nonalias_=nonalias_##nr;\
- int i_,imax_;\
- switch(phase ? phase->temp.declarePhase_ : 0) {\
- case 0: /* Handle literal start-values and parameters.*/\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int t=nonalias_[i_].tag&7;\
- switch(t) {\
- case 0:if (cuser_) cuser_[nonalias_[i_].ix]=(void*)(nonalias_[i_].d);break;\
- case 1:if (dp_) dp_[nonalias_[i_].ix]=nonalias_[i_].val;break;\
- case 2:if (x0_) x0_[nonalias_[i_].ix]=nonalias_[i_].val;break;\
- default:break;\
- }}break;\
- case 11: /* Count items */\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int storF=(nonalias_[i_].tag>>6)&(512|1024);\
- int stor=DYNStoreProtectedVariables || (storF==0 || storF==(512|1024));\
- int t=nonalias_[i_].tag&7;\
- switch(t) {\
- case 1:if(stor)phase->temp.isPublicP_[setDefaultP_]|=1;++setDefaultP_;break;\
- case 2:if(stor)phase->temp.isPublicX_[setDefaultX_]|=1;++setDefaultX_;break;\
- case 3:if(stor)phase->temp.isPublicXD_[setDefaultDX_]|=1;++setDefaultDX_;break;\
- case 6:if(stor)phase->temp.isPublicW_[setDefaultW_]|=1;++setDefaultW_;break;\
- case 7:{\
- static const int mapToO[8]={0,2,5,4,0,6,3,1};\
- int kind=mapToO[(nonalias_[i_].tag>>3)&7];\
- if (kind==6 && stor) phase->temp.isPublicW_[nonalias_[i_].ix]|=2;\
- if (kind==2 && stor) phase->temp.isPublicX_[nonalias_[i_].ix]|=2;\
- if (kind==3 && stor) phase->temp.isPublicXD_[nonalias_[i_].ix]|=2;\
- if (kind==1 && stor) phase->temp.isPublicP_[nonalias_[i_].ix]|=2;\
- }\
- default:break;\
- }}break;\
- case 1: /* Count items */\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int t=nonalias_[i_].tag&7;\
- int storF=(nonalias_[i_].tag>>6)&(512|1024);\
- int stor=DYNStoreProtectedVariables||(storF==0 || storF==(512|1024));\
- if (t>=1 && t<=6) {\
- if (!DYNStoreProtectedVariables && ((t==6 && !phase->temp.isPublicW_[setDefaultW_++])||(t==2 && !phase->temp.isPublicX_[setDefaultX_++])||(t==3 && !phase->temp.isPublicXD_[setDefaultDX_++])||(t==1 && !phase->temp.isPublicP_[setDefaultP_++]))) {phase->temp.offsets_[t]++;continue;}\
- } else if(t==7) {\
- static const int mapToO[8]={0,2,5,4,0,6,3,1};\
- int kind=mapToO[(nonalias_[i_].tag>>3)&7];\
- if (!DYNStoreProtectedVariables && (kind==6||kind==2||kind==3||kind==1) && !stor) continue;\
- if (!DYNStoreProtectedVariables && (kind==6 && phase->temp.isPublicW_[nonalias_[i_].ix]==2)) {phase->temp.isPublicW_[nonalias_[i_].ix]=4;continue;}\
- if (!DYNStoreProtectedVariables && (kind==2 && phase->temp.isPublicX_[nonalias_[i_].ix]==2)) {phase->temp.isPublicX_[nonalias_[i_].ix]=4;continue;}\
- if (!DYNStoreProtectedVariables && (kind==3 && phase->temp.isPublicXD_[nonalias_[i_].ix]==2)) {phase->temp.isPublicXD_[nonalias_[i_].ix]=4;continue;}\
- if (!DYNStoreProtectedVariables && (kind==1 && phase->temp.isPublicP_[nonalias_[i_].ix]==2)) {phase->temp.isPublicP_[nonalias_[i_].ix]=4;continue;}\
- }\
- switch(t) {\
- case 1:*((phase->temp.offsets_[t])++)= - ++(phase->temp.countConst_);break;\
- case 2:*((phase->temp.offsets_[t])++)= ++(phase->temp.countNonConst_);break;\
- case 3:*((phase->temp.offsets_[t])++)=++(phase->temp.countNonConst_);break;\
- case 4:case 5:case 6:if (nonalias_[i_].tag&32) *(phase->temp.offsets_[t]++)= - ++(phase->temp.countConst_);else *((phase->temp.offsets_[t])++) = ++(phase->temp.countNonConst_);break;\
- case 7:++(phase->temp.countAlias_);break;\
- default:break;\
- }}break;\
- case 3:/* Add names, info */\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int t=nonalias_[i_].tag&7;\
- int data=-1,offset,sgn=1;\
- int storF=(nonalias_[i_].tag>>6)&(512|1024);\
- int stor=DYNStoreProtectedVariables||(storF==0 || storF==(512|1024));\
- if (t>=1 && t<=7) {\
- if (!DYNStoreProtectedVariables && (t==6||t==2||t==3||t==1) && !stor) {phase->temp.offsets_[t]++;continue;}\
- if (t==7) {\
- static const int mapToO[8]={0,2,5,4,0,6,3,1};\
- int kind=mapToO[(nonalias_[i_].tag>>3)&7];\
- assert(kind!=0);\
- if (!DYNStoreProtectedVariables && (kind==6||kind==2||kind==3||kind==1) && !stor) continue;\
- offset=phase->temp.baseoffsets_[kind][nonalias_[i_].ix];\
- if (nonalias_[i_].val<0) sgn=-1;\
- } else offset=*((phase->temp.offsets_[t])++);\
- if (offset>0) {data=2;}else{offset=-offset;data=1;}\
- }\
- if (data!=-1) {\
- phase->temp.trajNames_[phase->temp.infoCount_]=nonalias_[i_].n;\
- phase->temp.trajDesc_[phase->temp.infoCount_]=nonalias_[i_].d;\
- phase->temp.trajInfo_[phase->temp.infoCount_+phase->totalInfoCount*0]=data;\
- phase->temp.trajInfo_[phase->temp.infoCount_+phase->totalInfoCount*1]=sgn*(offset+1);\
- phase->temp.trajInfo_[phase->temp.infoCount_+phase->totalInfoCount*2]=0;\
- phase->temp.trajInfo_[phase->temp.infoCount_+phase->totalInfoCount*3]=-1;\
- ++(phase->temp.infoCount_);\
- }\
- }break;\
- case 5: /* Update constants*/\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int t=nonalias_[i_].tag&7;\
- if (!DYNStoreProtectedVariables && ((t==6 && !phase->temp.isPublicW_[setDefaultW_++])||(t==2 && !phase->temp.isPublicX_[setDefaultX_++])||(t==3 && !phase->temp.isPublicXD_[setDefaultDX_++])||(t==1 && !phase->temp.isPublicP_[setDefaultP_++]))) {++(phase->temp.offsets_[t]);++(phase->temp.pointers_[t]);continue;}\
- if (t>=1 && t<=6) {\
- if (*(phase->temp.offsets_[t])<0) \
- phase->temp.resultPart_[2* - *(phase->temp.offsets_[t])+1] =phase->temp.resultPart_[2* - *(phase->temp.offsets_[t])] = (float)(*(phase->temp.pointers_[t]));\
- ++(phase->temp.offsets_[t]);\
- ++(phase->temp.pointers_[t]);\
- }\
- }break;\
- case 7:/*Update non-constants */\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int t=nonalias_[i_].tag&7;\
- if (t>=1 && t<=6) {\
- if (!DYNStoreProtectedVariables && ((t==6 && !phase->temp.isPublicW_[setDefaultW_++])||(t==2 && !phase->temp.isPublicX_[setDefaultX_++])||(t==3 && !phase->temp.isPublicXD_[setDefaultDX_++])||(t==1 && !phase->temp.isPublicP_[setDefaultP_++]))) {++(phase->temp.offsets_[t]);++(phase->temp.pointers_[t]);continue;}\
- if (*(phase->temp.offsets_[t])>0) \
- phase->temp.resultPart_[*(phase->temp.offsets_[t])] = (float)(*(phase->temp.pointers_[t]));\
- ++(phase->temp.offsets_[t]);\
- ++(phase->temp.pointers_[t]);\
- }\
- }break;\
- default:break;\
- }
- #else
- #define DeclareParameter(par, descpar, ipar, val, mi, ma, w,typ,con) if (dp_) dp_[ipar]=val;
- #define DeclareSParameter(par, descpar, ipar, val, mi, ma, w,typ,con) if (cuser_) cuser_[ipar]=(const char*)(val);
- #define DeclareState(x, descx, ix, val, mi, ma,w, typ, con) if (x0_) x0_[ix]=val;
- #define DeclareDerivative(dx, descdx, val, mi, ma,w, typ, con)
- #define DeclareOutput(var, descvar, iy, val, mi, ma, w,typ, con)
- #define DeclareInput(var, descvar, iu, val, mi, ma, w,typ, con) if (du_) du_[iu]=val;
- #define DeclareVariable(var,descvar, val, mi, ma, w, typ, con)
- #define DeclareAlias(name, description, aliasname, sign, sigtype, sigindex)
- #define DeclareAlias2(name, description, aliasname, sign, sigtype, sigindex, typ)
- #define PreNonAliasNew(nr) {declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define PreNonAlias(nr) {PreNonAliasDef(nr);declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define StartNonAlias(nr) FuncStartNonAlias(nr)
- #define EndNonAlias(nr)
- #endif
- #define DeclareConstant(c)
- #define StartAlias(nr) void declarea__##nr(double x0_[],double dp_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase) {
- #define EndAlias(nr)
- #define PreAlias(nr) {void declarea__##nr(double x0_[],double dp_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase);declarea__##nr(x0_,dp_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #else
- #define DeclareSParameter(par,descpar,ipar,val,mi,ma,w,typ,con) {0+(typ<<3)+(con<<5),par,val,0,0.0,0.0,w,ipar},
- #define DeclareParameter(par, descpar, ipar, val, mi, ma, w,typ,con) {1+(typ<<3)+(con<<5),par,descpar,val,mi, ma,w,ipar},
- #define DeclareState(x, descx, ix, val, mi, ma,w, typ,con) {2+(typ<<3)+(con<<5),x,descx,val,mi, ma, w, ix},
- #define DeclareDerivative(dx, descdx, val, mi, ma,w, typ,con) {3+(typ<<3)+(con<<5),dx,descdx,val,mi, ma, w},
- #define DeclareOutput(var, descvar, iy, val, mi, ma, w,typ, con) {4+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareInput(var, descvar, iu, val, mi, ma, w,typ, con) {5+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareVariable(var,descvar, val, mi, ma, w, typ, con) {6+(typ<<3)+(con<<5),var,descvar,val,mi, ma, w},
- #define DeclareAlias(name, description, aliasname, sign, sigtype, sigindex) {7+(sigtype<<3),name,description,sign,0.0,0.0,0,sigindex},
- #define DeclareAlias2(name, description, aliasname, sign, sigtype, sigindex, typ) {7+(sigtype<<3)+(typ<<6),name,description,sign,0.0,0.0,0,sigindex},
- #define PreNonAliasDef(nr) DYMOLA_STATIC void declare__##nr(double x0_[],double dp_[],double du_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase);
- #define PreNonAliasNew(nr) {declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define PreNonAlias(nr) {PreNonAliasDef(nr);declare__##nr(x0_,dp_,du_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define StartNonAlias(nr) static struct NonAlias_ nonalias_##nr []={
- #define EndNonAlias(nr) {100}}; FuncStartNonAlias(nr) \
- LIBDS_API void blo3val(int,int,double,double,double,double,int);\
- int c1_,c2_,c3_;double d1_;int i_,imax_;\
- struct NonAlias_*nonalias_=nonalias_##nr;\
- c1_=c2_=c3_=0;d1_=0;\
- if (!setDefault_) {\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- switch(nonalias_[i_].tag&7) {\
- case 0:blo3psd_(nonalias_[i_].n,nonalias_[i_].d);cuser_[nonalias_[i_].ix]=(void*)(nonalias_[i_].d);break;\
- case 1:blo3pd_(nonalias_[i_].n,nonalias_[i_].d);dp_[nonalias_[i_].ix]=nonalias_[i_].val;if (nonalias_[i_].tag>=64) blo3pt_(nonalias_[i_].tag>>6);break;\
- case 2:blo3sx_(nonalias_[i_].n,nonalias_[i_].d);x0_[nonalias_[i_].ix]=nonalias_[i_].val;if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- case 3:blo3sd_(nonalias_[i_].n,nonalias_[i_].d);if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- case 4:blo3sy_(nonalias_[i_].n,nonalias_[i_].d);if (nonalias_[i_].tag&32) {blo3co_(nonalias_[i_].val);};if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- case 5:blo3su_(nonalias_[i_].n,nonalias_[i_].d);if (nonalias_[i_].tag&32) {blo3co_(nonalias_[i_].val);};if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- case 6:blo3sw_(nonalias_[i_].n,nonalias_[i_].d);if (nonalias_[i_].tag&32) {blo3co_(nonalias_[i_].val);};if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- case 7:blo3sa_(nonalias_[i_].n,nonalias_[i_].d,(nonalias_[i_].tag>>3)&7,nonalias_[i_].ix+1,(nonalias_[i_].val>0)?1:-1);if (nonalias_[i_].tag>=64) blo3st_(nonalias_[i_].tag>>6);break;\
- default:break;\
- }\
- }} else {\
- for(i_=0,imax_=sizeof(nonalias_##nr)/sizeof(nonalias_##nr [0])-1;i_<imax_;++i_){\
- int num=0;int sect=nonalias_[i_].tag&7;\
- if (sect!=7 && sect!=0) {switch(sect) {\
- case 1:num=++setDefaultP_;break;\
- case 2:num=++setDefaultX_;break;\
- case 3:num=++setDefaultDX_;break;\
- case 4:num=++setDefaultY_;break;\
- case 5:num=++setDefaultU_;break;\
- case 6:num=++setDefaultW_;break;\
- default:break;\
- }blo3val(sect,num,nonalias_[i_].val,nonalias_[i_].mi,nonalias_[i_].ma,nonalias_[i_].w,(nonalias_[i_].tag>>3)&3);}\
- }}\
- #define PreAlias(nr) {void declarea__##nr(double x0_[],double dp_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase);declarea__##nr(x0_,dp_,cuser_,QiErr,setDefault_,setDefaultX_,setDefaultU_,setDefaultY_,setDefaultP_,setDefaultDX_,setDefaultW_,phase);}}
- #define StartAlias(nr) static struct NonAlias_ alias_##nr []={
- #define EndAlias(nr) {100}}; void declarea__##nr(double x0_[],double dp_[],void*cuser_[],int*QiErr,int setDefault_,int setDefaultX_,int setDefaultU_,int setDefaultY_,int setDefaultP_,int setDefaultDX_,int setDefaultW_, struct DeclarePhase*phase) {\
- LIBDS_API void blo3val(int,int,double,double,double,double,int);\
- int c1_,c2_,c3_;double d1_;int i_,imax_;\
- struct NonAlias_*alias_=alias_##nr;\
- c1_=c2_=c3_=0;d1_=0;\
- if (!setDefault_) {\
- for(i_=0,imax_=sizeof(alias_##nr)/sizeof(alias_##nr [0])-1;i_<imax_;++i_){\
- switch(alias_[i_].tag&7) {\
- case 0:blo3psd_(nonalias_[i_].n,nonalias_[i_].d);cuser_[nonalias_[i_].ix]=(void*)(nonalias_[i_].d);break;\
- case 1:blo3pd_(alias_[i_].n,alias_[i_].d);dp_[alias_[i_].ix]=alias_[i_].val;if (alias_[i_].tag>=64) blo3pt_(alias_[i_].tag>>6);break;\
- case 2:blo3sx_(alias_[i_].n,alias_[i_].d);x0_[alias_[i_].ix]=alias_[i_].val;if (alias_[i_].tag>=64) blo3st_(alias_[i_].tag>>6);break;\
- case 3:blo3sd_(alias_[i_].n,alias_[i_].d);if (alias_[i_].tag>=64) blo3st_(alias_[i_].tag>>6);break;\
- case 4:blo3sy_(alias_[i_].n,alias_[i_].d);if (alias_[i_].tag&32) {blo3co_(alias_[i_].val);};if (alias_[i_].tag>=64) blo3st_(alias_[i_].tag>>6);break;\
- case 5:blo3su_(alias_[i_].n,alias_[i_].d);if (alias_[i_].tag&32) {blo3co_(alias_[i_].val);};if (alias_[i_].tag>=64) blo3st_(alias_[i_].tag>>6);break;\
- case 6:blo3sw_(alias_[i_].n,alias_[i_].d);if (alias_[i_].tag&32) {blo3co_(alias_[i_].val);};if (alias_[i_].tag>=64) blo3st_(alias_[i_].tag>>6);break;\
- case 7:blo3sa_(alias_[i_].n,alias_[i_].d,(alias_[i_].tag>>3)&7,alias_[i_].ix+1,(alias_[i_].val>0)?1:-1);break;\
- default:break;\
- }\
- }} else {\
- for(i_=0,imax_=sizeof(alias_##nr)/sizeof(alias_##nr [0])-1;i_<imax_;++i_){\
- int num=0;int sect=alias_[i_].tag&7;\
- if (sect!=7 && sect!=0) {switch(sect) {\
- case 1:num=++setDefaultP_;break;\
- case 2:num=++setDefaultX_;break;\
- case 3:num=++setDefaultDX_;break;\
- case 4:num=++setDefaultY_;break;\
- case 5:num=++setDefaultU_;break;\
- case 6:num=++setDefaultW_;break;\
- default:break;\
- }blo3val(sect,num,alias_[i_].val,alias_[i_].mi,alias_[i_].ma,alias_[i_].w,(alias_[i_].tag>>3)&3);}\
- }}\
- #endif
- #if defined(DYMOLA_DSPACE) || defined(DynSimStruct) && DymolaAutoRemoveAuxiliaries_
- #define AcceptedSection2\
- }sectioncondition=0; if (0) {*icall_=2;
- #elif DymolaAutoRemoveAuxiliaries_
- #define AcceptedSection2\
- } sectioncondition=sectioncondition && DymolaStoreAuxiliaries_; if (sectioncondition) {
- #else
- #define AcceptedSection2
- #endif
- #define CrossingSection\
- }\
- sectioncondition = (*idemand_ >= 4 && *icall_ < 4) || (DYNEvent || *idemand_==4);\
- if (sectioncondition) {\
- *icall_ = 4;
- #if !defined(DYMOSIM)
- #define SetInitVectorSimple(x, i, var, val) \
- if (veryFirstTime && lambda_<-1) \
- var = val; \
- if ((firstTime || ! DAEsolver_) && !inJacobian_) \
- SetVector(x, i, var); \
- else \
- var = GetVector(x, i)
- #define SetInitVectorSimpleNH(x, i, var, val) \
- if (veryFirstTime && lambda_<-1) \
- var = val; \
- if ((firstTime || ! DAEsolver_) && !inJacobian_) \
- SetVector(x, i, var); \
- else \
- var = GetVector(x, i)
- #else
- #define SetInitVectorSimpleOld(x, i, var, val) \
- if (veryFirstTime && (!DymolaUserHomotopy || DymolaUserHomotopy && DymolaHomotopyLambda==0)) var=val;\
- if ((firstTime || ! DAEsolver_) && !inJacobian_) \
- SetVector(x, i, var); \
- else \
- var = GetVector(x, i)
- #define SetInitVectorSimple(x, i, var, val) \
- if (veryFirstTime && (!DymolaUserHomotopy || DymolaUserHomotopy && DymolaHomotopyLambda==0)) {\
- if (veryFirstTime2) SetVector(xInitial, i, var);else var = GetVector(xInitial, i);\
- }\
- if (firstTime || ! DAEsolver_) \
- SetVector(x, i, var); \
- else \
- var = GetVector(x, i)
- #define SetInitVectorSimpleNH(x, i, var, val) \
- if (veryFirstTime && !(homlocal>0 || DymolaUserHomotopy && DymolaHomotopyLambda>0)) var=val;\
- if ((firstTime || ! DAEsolver_) && !inJacobian_) \
- SetVector(x, i, var); \
- else \
- var = GetVector(x, i)
- #endif
- #if defined(DYMOLA_DSPACE) || defined(NO_FILE)
- #define InitializeModes(datfile,Modes,nModes) {nModes=0;}
- #if 0
- /* Change to #if 1 to see modes in CPUClk */
- #define AddMode(Modes,nModes,Mode) { \
- int j; \
- unsigned int cnt = 0; \
- for(j = 1; j <= MatrixNrRows(Modes); j++) \
- cnt = (cnt << 1) | (GetVector(Mode,j) != 0); \
- CPUClk = cnt; \
- }
- #else
- #define AddMode(Modes,nModes,Mode)
- #endif
- #define OutputModes
- #else
- /* To handle modes */
- #define AddMode(Modes,nModes,Mode) {\
- int i,j;\
- int doAppend=(nModes<MaxModes);\
- for(i=1;doAppend&&(i<=nModes);i++) {\
- int allSame=1;\
- for(j=1;allSame&&(j<=MatrixNrRows(Modes));j++) \
- allSame=GetMatrix(Modes,j,i)==GetVector(Mode,j);\
- doAppend=!allSame;\
- }\
- if (doAppend) {\
- nModes++;\
- for(j=1;j<=MatrixNrRows(Modes);j++)\
- SetMatrix(Modes,j,nModes,GetVector(Mode,j)); \
- }\
- }
- #define InitializeModes(datfile,Modes,nModes) {\
- Amatrix bigModes; AclassRead c[1];\
- bigModes.name="bigModes"; bigModes.nrow=MatrixNrRows(Modes);\
- bigModes.ncol=-MatrixNrColumns(Modes); bigModes.type=doubleMatrix; bigModes.data.d=Modes.data;\
- c[0].matrix=&bigModes; c[0].req=1; c[0].rowDim=0; c[0].colDim=0; c[0].rowIndex=0; c[0].colIndex=0; c[0].onFile=0;\
- nModes=0;\
- if (amatReadAll(datfile,"NoClass",(char*)0,c,1)==0) {\
- nModes=bigModes.ncol;\
- } else {\
- DymosimMessage(amatError);\
- }\
- }
- #define OutputModes(datfile,mofile,newmodel,origmodel,pre,Modes,nModes,nModesOriginal) {\
- if (datfile) {\
- Amatrix bigModes;\
- bigModes.name="bigModes"; bigModes.nrow=MatrixNrRows(Modes);\
- bigModes.ncol=nModes; bigModes.type=doubleMatrix; bigModes.data.d=Modes.data;\
- amatWrite(datfile,amatBinary,bigModes);\
- }\
- if (mofile) {FILE *f=fopen(mofile,"w");\
- if (f) {\
- int i,j,jmax=MatrixNrRows(Modes);\
- fprintf(f,"// Automatically constructed mode variables\n\n");\
- fprintf(f,"model %s\n\n",newmodel);\
- fprintf(f," extends %s;\n\n",origmodel);\
- fprintf(f," %s\n\n",pre);\
- fprintf(f," Real bigModes[%d,%d]={",nModes,jmax);\
- for(i=1;i<=nModes;i++) {fprintf(f,"\n { ");for(j=1;j<=jmax;j++) \
- fprintf(f,"%s%g",(j==1)?"":", ",GetMatrix(Modes,j,i));fprintf(f,(i==nModes)?"}":"},");}\
- fprintf(f," };\n\n");\
- fprintf(f,"end %s;\n\n",newmodel);\
- }\
- fclose(f);}\
- if (nModes>nModesOriginal) DymosimMessage("Found additional modes, please recompile model");\
- }
- #include <bloutil.h>
- #ifdef Dymola_Include_MT
- /*
- A C-program for MT19937, with initialization improved 2002/1/26.
- Coded by Takuji Nishimura and Makoto Matsumoto.
- Before using, initialize the state by using init_genrand(seed)
- or init_by_array(init_key, key_length).
- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
- All rights reserved.
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
- 1. Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
- 2. Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
- 3. The names of its contributors may not be used to endorse or promote
- products derived from this software without specific prior written
- permission.
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- Any feedback is very welcome.
- http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
- email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
- Modified to include prefix; added code for using time to randomize values (and print it)
- and removed unused functions.
- */
- /* Period parameters */
- #define Dymola_MT_N 624
- #define Dymola_MT_M 397
- #define Dymola_MT_MATRIX_A 0x9908b0dfUL /* constant vector a */
- #define Dymola_MT_UPPER_MASK 0x80000000UL /* most significant w-r bits */
- #define Dymola_MT_LOWER_MASK 0x7fffffffUL /* least significant r bits */
- static unsigned int Dymola_MT_mt[Dymola_MT_N]; /* the array for the state vector */
- static int Dymola_MT_mti=Dymola_MT_N+1; /* mti==N+1 means mt[N] is not initialized */
- /* initializes mt[N] with a seed */
- static void Dymola_MT_init_genrand(unsigned int s)
- {
- Dymola_MT_mt[0]= s & 0xffffffffUL;
- for (Dymola_MT_mti=1; Dymola_MT_mti<Dymola_MT_N; Dymola_MT_mti++) {
- Dymola_MT_mt[Dymola_MT_mti] =
- (1812433253UL * (Dymola_MT_mt[Dymola_MT_mti-1] ^ (Dymola_MT_mt[Dymola_MT_mti-1] >> 30)) + Dymola_MT_mti);
- /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
- /* In the previous versions, MSBs of the seed affect */
- /* only MSBs of the array mt[]. */
- /* 2002/01/09 modified by Makoto Matsumoto */
- Dymola_MT_mt[Dymola_MT_mti] &= 0xffffffffUL;
- /* for >32 bit machines */
- }
- }
- /* initialize by an array with array-length */
- /* init_key is the array for initializing keys */
- /* key_length is its length */
- /* slight change for C++, 2004/2/26 */
- static void Dymola_MT_init_by_array(unsigned int init_key[], int key_length)
- {
- int i, j, k;
- Dymola_MT_init_genrand(19650218UL);
- i=1; j=0;
- k = (Dymola_MT_N>key_length ? Dymola_MT_N : key_length);
- for (; k; k--) {
- Dymola_MT_mt[i] = (Dymola_MT_mt[i] ^ ((Dymola_MT_mt[i-1] ^ (Dymola_MT_mt[i-1] >> 30)) * 1664525UL))
- + init_key[j] + j; /* non linear */
- Dymola_MT_mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++; j++;
- if (i>=Dymola_MT_N) { Dymola_MT_mt[0] = Dymola_MT_mt[Dymola_MT_N-1]; i=1; }
- if (j>=key_length) j=0;
- }
- for (k=Dymola_MT_N-1; k; k--) {
- Dymola_MT_mt[i] = (Dymola_MT_mt[i] ^ ((Dymola_MT_mt[i-1] ^ (Dymola_MT_mt[i-1] >> 30)) * 1566083941UL))
- - i; /* non linear */
- Dymola_MT_mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++;
- if (i>=Dymola_MT_N) { Dymola_MT_mt[0] = Dymola_MT_mt[Dymola_MT_N-1]; i=1; }
- }
- Dymola_MT_mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
- }
- #include <time.h>
- #ifdef _MSC_VER
- #include <process.h>
- #endif
- static void Dymola_MT_init_default(void) {
- char str[100];
- #if 0
- Dymola_MT_init_genrand(5489UL); /* a default initial seed is used */
- #else
- unsigned int s =(unsigned int)(time(0));
- #ifdef _MSC_VER
- s+=(unsigned int)(_getpid())*10235489UL;
- #endif
- Dymola_MT_init_genrand(s);
- DymosimMessageInt("Random number seed for MT19937:",s);
- #endif
- }
- /* generates a random number on [0,0xffffffff]-interval */
- static unsigned int Dymola_MT_genrand_int32(void)
- {
- unsigned int y;
- static unsigned int mag01[2]={0x0UL, Dymola_MT_MATRIX_A};
- /* mag01[x] = x * MATRIX_A for x=0,1 */
- if (Dymola_MT_mti >= Dymola_MT_N) { /* generate N words at one time */
- int kk;
- if (Dymola_MT_mti == Dymola_MT_N+1) /* if init_genrand() has not been called, */
- Dymola_MT_init_default(); /* a default initial seed is used */
- for (kk=0;kk<Dymola_MT_N-Dymola_MT_M;kk++) {
- y = (Dymola_MT_mt[kk]&Dymola_MT_UPPER_MASK)|(Dymola_MT_mt[kk+1]&Dymola_MT_LOWER_MASK);
- Dymola_MT_mt[kk] = Dymola_MT_mt[kk+Dymola_MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
- }
- for (;kk<Dymola_MT_N-1;kk++) {
- y = (Dymola_MT_mt[kk]&Dymola_MT_UPPER_MASK)|(Dymola_MT_mt[kk+1]&Dymola_MT_LOWER_MASK);
- Dymola_MT_mt[kk] = Dymola_MT_mt[kk+(Dymola_MT_M-Dymola_MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
- }
- y = (Dymola_MT_mt[Dymola_MT_N-1]&Dymola_MT_UPPER_MASK)|(Dymola_MT_mt[0]&Dymola_MT_LOWER_MASK);
- Dymola_MT_mt[Dymola_MT_N-1] = Dymola_MT_mt[Dymola_MT_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
- Dymola_MT_mti = 0;
- }
-
- y = Dymola_MT_mt[Dymola_MT_mti++];
- /* Tempering */
- y ^= (y >> 11);
- y ^= (y << 7) & 0x9d2c5680UL;
- y ^= (y << 15) & 0xefc60000UL;
- y ^= (y >> 18);
- return y;
- }
- /* generates a random number on [0,1) with 53-bit resolution*/
- static double Dymola_MT_genrand_res53(void)
- {
- unsigned int a=Dymola_MT_genrand_int32()>>5, b=Dymola_MT_genrand_int32()>>6;
- return(a*67108864.0+b)*(1.0/9007199254740992.0);
- }
- /* These real versions are due to Isaku Wada, 2002/01/09 added */
- /* End of MT19937*/
- #endif
- #ifndef DYMOLA_TIMES
- #define DYMOLA_TIMES
- struct DymolaTimes {
- int num;
- double maxim;
- double minim;
- double total;
- int numAccum;
- double maximAccum;
- double minimAccum;
- double totalAccum;
- const char*name;
- };
- #endif
- #if defined(DYMOLA_DSPACE) || defined(NO_FILE) || defined(FMU_SOURCE_CODE_EXPORT) || defined(LABCAR)
- #define DymolaStartTimer(x)
- #define DymolaEndTimer(x)
- #define DymolaEndTimerName(x,n)
- #define DymolaResetTimer(x) 0
- #else
- double DymolaTimerCounterSince(double*,int);
- extern void DymolaOnInterrupt(void);
- extern int DymolaNoInterrupt(void);
- static void dymolaEndTimer(int num,struct DymolaTimes*dt,double Time,double*start,const char*name,double*tz,int tz_len,int*tz_cnt) {
- double diff=DymolaTimerCounterSince(start,1);
- if (num==0) {
- DymolaOnInterrupt();
- }
- if (!dt) return; /* Error */
- #if !defined(DymosimRealTimePriority_)
- if (dt->num<=10 && 2*diff<dt->minim && diff>0) {
- dt->num=0;
- dt->numAccum=0;
- dt->total=0;
- dt->totalAccum=0;
- }
- #endif
- if (dt->num==0) {
- dt->minim=diff;
- dt->maxim=diff;
- dt->name=name;
- } else {
- #if !defined(DymosimRealTimePriority_)
- if (diff>2000*dt->maxim && dt->maxim>0)
- return;
- #endif
- if (diff>dt->maxim) dt->maxim=diff;
- if (diff<dt->minim) dt->minim=diff;
- }
- dt->total+=diff;
- #if defined(DymosimRealTimePriority_)
- if (num==0 && tz && tz_cnt && 2*(*tz_cnt)+2<tz_len) {
- if (tz[2*(*tz_cnt)]!=Time)
- ++(*tz_cnt);
- tz[2*(*tz_cnt)]=Time;
- tz[2*(*tz_cnt)+1]+=diff;
- }
- #endif
- dt->num++;
- if (dt->numAccum==0) {
- dt->maximAccum=diff;
- dt->minimAccum=diff;
- } else {
- if (diff>dt->maximAccum) dt->maximAccum=diff;
- if (diff<dt->minimAccum) dt->minimAccum=diff;
- }
- dt->totalAccum+=diff;
- dt->numAccum++;
- }
- #define DymolaStartTimer(x) ((x==0)&&((DymolaNoInterrupt()))),DymolaTimerCounterSince(did_->DymolaStartTimers_vec+x,0)
- #define DymolaEndTimer(x) dymolaEndTimer(x,did_->DymolaTimerStructs_vec+x,DYNTime,did_->DymolaStartTimers_vec+x,"",did_->DymolaTimeZero_vec,did_->DymolaTimeZeroLength_var,&(did_->DymolaTimecounter_var))
- #define DymolaEndTimerName(x,n) dymolaEndTimer(x,did_->DymolaTimerStructs_vec+x,DYNTime,did_->DymolaStartTimers_vec+x,n,did_->DymolaTimeZero_vec,did_->DymolaTimeZeroLength_var,&(did_->DymolaTimecounter_var))
- #define DymolaResetTimer(x) did_->DymolaTimerStructs_vec[x].totalAccum=did_->DymolaTimerStructs_vec[x].num=did_->DymolaTimerStructs_vec[x].numAccum=0
- #endif
- #endif
- #endif
|