module mStiffbgh;
use zone_environment_types, only: zone_environment_set_type;
use cp_error_handling, only: cp_error_type;
use mHcdfnrad, only: difmat;
implicit real*8 (A-H,O-Z);
-- implicit none;
private;
character(len=*), parameter, private :: mdl_name = 'mStiffbgh';
public vtstif;
contains;
-- GEAR integrator for ODE with sparse Jacobian and option METH=3
-- for BRAYTON, GUSTAVSON & HACHTEL method
_Define
@DNAMEDVAR DNDVAR -- work array of derivatives
@STRING STR020 -- for debugger
@LINE LIN130 -- for debugger
@THR 1.D-12 -- Zlatev threshold
-- @THR 1.D-15 -- Zlatev threshold
--_TRACE "WRITE(@Wres,*)' ifail =',ifail,' NumSing=',NUMSING,"
--_TRACE "WRITE(@Wres,*)' Y11 =',Y(1,1),' Ytp11 =',Y(301,1),"
--_TRACE "IF(NSTEP>8000) write(*,'(a,L3,1p,2e12.3,a)')'EVALJA H Hused=',EVALJA,H,Hused,"
--_TRACE "@wterm' lmax=',lmax,' maxord=', maxord,"
--_TRACE "@wterm' EVALJA=',EVALJA,' maxer=', maxer,"
--_TRACE "write(@term,'(a,4L3,i4,a)')' EVALJA needbr corrco badste ifail=',
-- EVALJA, needbr, corrco, badste, ifail,"
--_TRACE "@wterm' EVALJA=',EVALJA,' ifail =',ifail,' H=',H,"
--_TRACE "WRITE(*,*)' N=',N,' Nzmod=',Nzmod,"
subroutine vtstif(envs, err);
IMPLICIT REAL*8(A-H,O-Z);
/* VTSTIF performs one step of the integration of an initial
value problem for a system of ordinary differential equations.
Communication with VTSTIF is done with the following variables:
Y an NYDIM by LMAX array containing the dependent variables
and their scaled derivatives. LMAX is currently 13 for
the Adams methods and 6 for the Gear methods (though
LMAX=5 is recommended due to stability considerations).
LMAX-1 is MAXORD, the maximum order used. See subroutine
COSET.
Y(I,J+1) contains the J-th derivative of Y(I), scaled
by H**J/Factorial(J). Only Y(I), 1 <= I <= N, need be set
by the calling program on the first entry.
If it is desired to interpolate to non-mesh points,
the Y array can be used. If the current step size is H
and the value at t+E is needed, form S = E/H, and then
compute
NQ
Y(I)(t+E) = SUM Y(I,J+1)*S**J .
J=0
The Y array should not be altered by the calling program.
When referencing Y as a 2-dimensional array, use a column
length of NYDIM, as this is the value used in STIFF.
N The number of first order differential equations. N may be
decreased on later calls if the number of active equations
reduces, but it must not be increased without calling with
JSTART=0.
NYDIM A constant integer >= N, used for dimensioning purposes.
NYDIM must not be changed without setting JSTART=0.
t The independent variable. t is updated on each step taken.
H The step size to be attempted on the next step. H may be
adjusted up or down by the routine in order to achieve
an economical integration. However, if the H provided by
the user does not cause a larger error than requested, it
will be used. To save computer time, the user is advised
to use a fairly small step for the first call. It will be
automatically increased later. H can be either positive
or negative, but its sign must remain constant throughout
the problem.
HMIN The minimum absolute value of the step size that will be
used for the problem. On starting this must be much smaller
than the average Abs(H) expected, since a first order
method is used initially.
HMAX The maximum absolute value of the step size that will be
used for the problem.
EPSJ The relative error test constant. Single step error estimates
divided by YMAX(I) must be less than this in the Euclidean
norm. The step and/or order is adjusted to achieve this.
METH The method flag.
METH=1 means the implicit ADAMS methods
METH=2 means the GEAR method for stiff problems
METH=3 means the GEAR method, corrected by
R.K.Brayton, F.G.Gustavson & G.D.Hachtel in the
Proceedings of the IEEE v.60, No.1, January 1972, p.98-108
YMAX An array of N locations which contains the maximum absolute
value of each Y seen so far.
ERROR An array of N elements proportional to the estimated one step
error in each component.
KFLAG A completion code with the following meanings:
0 the step was succesful
-1 the requested error could not be achieved with
abs(H)=HMIN
-2 corrector convergence could not be achieved for
abs(H)>HMIN.
On a return with KFLAG negative, the values of t and the Y
array are as of the beginning of the last step, and H is
the last step size attempted.
JSTART An integer used on input and output.
On input it has the following values and meanings:
0 perform the first step. This value enables the
subroutine to initialize itself.
>0 take a new step continuing from the last.
Assumes the last step was succesful and user has not
changed any parameters.
<0 repeat the last step with a new value of H and/or
EPS and/or METH. This may be either in redoing a step
that failed, or in continuing from a succesful step.
On exit, JSTART is set to NQ, the current order of the
method. This is also the order of the maximum derivative
available in the Y array. After a succesful step, JSTART
need not be reset for the next call.
AJAC A block used for partial derivatives. It keeps the
matrix AJAC=1-EL(1)*H*Jacobian & is computed by DIFMAT.
HL = H*EL(1) is communicated to DIFMAT to compute AJAC.
FSAVE A block of at least 2*NYDIM locations for temporary storage.
The parameters which must be input by the user are:
N, NYDIM, t, Y, H, HMIN, HMAX, EPS, METH, JSTART, MAXORD.
Additional Subroutines required are:
DIFMAT computes DY/Dt given t and Y, and if EVALJA computes AJAC,
stored in the form appropriate for M28Y12,
COSET,
RESCAL,
M28Y12, M30Y12, M28BYS - sparse matrix solvers.
The calling program must contain the common declarations given in
node %C:
*/
type(zone_environment_set_type), pointer :: envs;
type(cp_error_type), intent(inout), optional :: err;
character(len=*), parameter :: subrtn_name = 'vtstif'
, fullPathSubrtn = mdl_name//'.'//subrtn_name;
_INCLUDE snrad;
<*C: COMMONS & VAR STIFF *>;
<*B: DATA FOR STIFF *>;
<*D: DATA FOR SPARSE MATRIX SubroutineS M28Y12 *>;
<*F: PREPARATIONS DEPENDING ON JSTART *>;
HNT(1)=H;
CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);-- needed for continuation
--runs
RC=RC*EL(1)/OLDL0;OLDL0=EL(1);
-- RC is the ratio of new to old values
-- of the coefficient H*EL(1).
-- When RC differs from 1 by more than RCTEST,
-- EVALJA is set to .TRUE. to force
-- jacobian to be calculated in DIFMAT.
<*A: the constants E,EUP,EDN,BND *>;
IF(MOD(NSTEP,MBATCH)==0) NEEDBR=.TRUE.;
RETCOD=0;
<*G: GENERAL STEP OF STIFF WHILE RETCOD==0 *>;
<*X: ALL RETURNS ARE MADE THROUGH THIS SECTION DEPENDING ON RETCOD.
H IS USED IN HOLD TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. *>;
RETURN;
end subroutine vtstif;
end module mStiffbgh;
%C:O
_include nstep;
-- COMMON/NSTEP/NSTEP,NDebug,MAXER,IOUT,NOUT;
-- COMMON/CREAD/TAUOLD,NSTMAX,MBATCH,MAXORD;
_include stsave;
-- COMMON/STSAVE/RMAX,TREND,OLDL0,RC,HOLD,EDN,E,EUP,BND,EPSOLD,TOLD,
-- MEO,NOLD,NQ,LNQ,IDOUB;
REAL*8 ERROR(NYDIM),EL(13),TQ(4);
REAL*8 XMOD,ERREST,ERROLD;
PARAMETER (NCHANL=6,
ETA=1.D-04); -- PIVOT THRESHOLD IN M28BYS
PARAMETER (DeltaDeb=1.d-10); -- FOR DEBUGGER
CHARACTER @STRING*20,@LINE*130; -- FOR DEBUGGER
CHARACTER*20 CharVarName(5); -- FOR DEBUGGER
Logical FOUNDI; -- FOR DEBUGGER
REAL*8 @DNAMEDVAR(NYDIM),DSAVE(NYDIM);-- FOR DEBUGGER
Integer RETCOD; -- RETURN CODE
Real Tm1,Tm2; -- Timer for CMS -- CRAY REAL
LOGICAL CORRCO; -- CORRECTOR ITERATIONS HAVE CONVERGED
LOGICAL NUMSING; -- Numerical singularity in decomposition
<*L: COMMON & VARIABLES FOR LINEAR EQS SOLVER F01BRF or M28Y12 *>;
%CL:O
Parameter(MAXIT=15);
LOGICAL LBLOCK,GROW,ABORT(4);
COMMON/JSAVE/ASAVE(NZ),IRS(NZ),ICS(NZ);
Dimension DB(NYDIM), XSAVE(NYDIM);
%B:O
DATA ANOISE/1.D-12/;--should be set to the noise level of the machine
DATA DFLTZR/1.D-75/, MAXITE/3/, MAXFAI/3/;
DATA RMXINI/1.D+4/, RMXNOR /10.D0/, RMXFAI/2.0D0/, IDELAY /10/;
DATA RHCORR/.25D0/, RHERR3/.1D0/, RCTEST/.3D0/;
DATA BIAS1/1.3D0/, BIAS2/1.2D0/, BIAS3/1.4D0/;
DATA CharVarName/'radius','velocity','temperature','FJ','FH'/;
%D:O
DATA LBLOCK /.True./,GROW/.TRUE./,
ABORT/.FALSE.,.TRUE.,.FALSE.,.TRUE./,
PIVOT/1.D-01/;
%F:
KFLAG=0;tOLD=t; NUMSING=.FALSE.;
-- WRITE(@term,*)'Entering Stiff NSTEP=', NSTEP;
/* On the first call, the order is set to 1 and the initial
derivatives are calculated. RMAX is the maximum ratio
by which H can be increased in a single step. It is
initially RMXINI to compensate for the small initial H,
but then is normally equal to RMXNOR. If a failure occurs
(in corrector convergence or error test), RMAX is set
at RMXFAI for the next increase. */
_SELECT
_ JSTART==0 [<*B: INITIAL *> ]
_ JSTART<0 [
<*C: reset Y & other variables for changed
METH, EPS, N, H *> ]
_END;
%FB:
HUSED=0.D0; NQUSED=0; NOLD=N; EVALJA=.FALSE.;
NQ=1;
_DO I=1,N; FSAVE(I)=Y(I,1)_OD;
HL=H*EL(1);
CALL DIFMAT(envs, err);
_DO I=1,N; Y(I,2)=FSAVE(NYDIM+I)*H _OD;
LNQ=2; IDOUB=LNQ+1; RMAX=RMXINI; EPSJ=EPS*SQRT(Real(N));
TREND=1.D0; OLDL0=1.D0; RC=0.D0;
-- RC is the ratio of new to old values
-- of the coefficient H*EL(1).
-- When RC differs from 1 by more than RCTEST,
-- EVALJA is set to .TRUE. to force
-- Jacobian to be calculated in DIFMAT.
HOLD=H; MEO=METH;
EVALJA=.TRUE.; OLDJAC=.FALSE.;
LMAX=MAXORD+1;
EPSOLD=EPSJ;
%A:
EDN=(TQ(1)*EPSJ)**2;-- is to test for decreasing the order,
E=(TQ(2)*EPSJ)**2; -- comparison for errors of current order NQ,
EUP=(TQ(3)*EPSJ)**2;-- is to test for increasing the order,
BND=(TQ(4)*EPSJ)**2;-- is used to test for convergence of the
-- correction iterates
%FC:
-- If the caller change EPSJ or METH , the constants E,EUP,EDN,BND
-- must be reseted.
EPSJ=EPS*SQRT(Real(N));
IF(METH==MEO) THEN;
IF(EPSJ^=EPSOLD)THEN;
EPSOLD=EPSJ;
ENDIF;
ELSE;
IDOUB=LNQ+1;
ENDIF;
<*C: TEST N==NOLD, H==HOLD ? *>;EPSOLD=EPSJ;
%FCC:
IF(N^=NOLD)THEN;
IDOUB=LNQ+1;EVALJA=.TRUE.;OLDJAC=.FALSE.;NOLD=N;
ENDIF;
IF(H^=HOLD) THEN;
RH=H/HOLD;H=HOLD;
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
ENDIF;
%G:
_WHILE RETCOD==0 _DO
IF(ABS(RC-1.D0)>RCTEST .or. MOD(NSTEP,MBATCH)==0) EVALJA=.TRUE.;
If(ifail==13)Then;
Needbr=.true.;
Evalja=.true.;
Endif;
t=t+h;
<*E: Predictor.
This section computes the predicted values by effectively
multiplying the Y array by the Pascal triangle matrix *>
_Repeat -- loop until .not.EVALJA
ITER=0;
_DO I=1,N; ERROR(I)=0.D0 _OD;
_DO I=1,N; FSAVE(I)=Y(I,1) _OD;
CORRCO=.FALSE.; -- Corrector convergence. It is tested
-- by requiring changes relative to YMAX(I)
-- to be less, in Euclidian norm, than BND,
-- which is dependent on EPSJ.
_WHILE ^CORRCO & ITER<MAXITE _DO
<*Cor: Up to MAXITE corrector iterations are taken.
The sum of the corrections is accumulated in the
vector ERROR(i). It is approximately equal to the
Lnq-th derivative of Y multiplied by
H**Lnq/(factorial(Lnq-1)*EL(Lnq)),
and is thus proportional to the actual errors
to the lowest power of H present (H**Lnq).
The Y array is not altered in the correction
loop. The updated Y vector is stored temporarily
in FSAVE. The norm of the iterate difference
is stored in D. *>
_OD; -- CORRCO ! ITER==MAXITE
_IF CORRCO _THEN
<*X: The corrector has converged. OLDJAC is set to .true.
if partial derivatives were used, to signal that they
may need updating on subsequent steps. The error test
is made and control passes to node %GXA if it
succeeds. *>;
_ELSE -- ITER==MAXITE
-- The corrector iterations failed to converge in MAXITE
-- tries.
NITER=NITER+1;
IF(OLDJAC)THEN;
-- If partials are not up to date , they are
-- reevaluated for the next try.
EVALJA=.TRUE.; -- Refresh JACOBIAN
ELSE;
<*Y: Otherwise, the Y array is retracted to its values
before prediction, and H is reduced, if possible.
If not, a no convergence exit is taken. *>;
ENDIF;
_FI;
_Until ^EVALJA .or. RETCOD^=0;
_OD; -- while RETCOD==0
%GE:
do J1=1,NQ;
do J2=J1,NQ;
J=NQ-J2+J1;
do I=1,N;
Y(I,J)=Y(I,J)+Y(I,J+1);
enddo;
enddo;
enddo; --_OD_OD_OD;
%G_Cor:
BADSTE=.FALSE.; -- Logical variable BADSTE may be changed by DIFMAT
HL=H*EL(1);
-- if(Nstep>1750 & EVALJA)then;
-- if(Nstep>1586 & EVALJA)then;
if(Nstep>=NDebug & EVALJA)then;
<*DEB: full debugger of Jacobian *>;
end if;
if(ifail==13)then;
@wterm ' G_Cor: ifail 13';
endif;
CALL DIFMAT(envs, err); -- if necessary (i.e. EVALJA=.true.) the partials
-- are reevaluated by DIFMAT
<*MatrOut : output Jacobian matrix *>;
-- stop ' Jacobian is output';
_IF BADSTE _THEN
-- Something is wrong in DIFMAT (i.e. Y(I) is out of physical
-- range for some I)
WRITE(*,*)' G_Cor: BADSTE=',BADSTE;
ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.;
Ifail=13;
-- Variables are given the values such as to guarantee that
-- control passes to the node %GY - to try a new smaller step
_ELSE
_DO Isw=1,N;
XSAVE(Isw)=0.D0;
_OD;
IF(EVALJA) THEN;
<*L: L-U Decomposition of matrix AJAC=1-El*H*Jacobian *>;
EVALJA=.FALSE.;OLDJAC=.FALSE.;RC=1.D0;
NJAC=NJAC+1;
ENDIF;
<*U: Find in FSAVE(1:N) the solution X of equation
AJAC!X>=H*F(Y)-Y(I,2)-ERROR(I),
where AJAC is the matrix AJAC=1-El*H*Jacobian,
then obtain in ERROR the full correction,
in D - the norm of the iterate difference &
in FSAVE a new approximation to Y *>
If(Ifail^=13)Then;
IF(ITER^=0) TREND=MAX(.9D0*TREND,D/D1);
-- If ITER > 0 an estimate of the convergence rate constant
-- is stored in TREND, and this is used in the convergence test
CORRCO=D*MIN(1.D0,2.D0*TREND)<=BND;
D1=D;
ITER=ITER+1;
Else;
--WRITE(@Wres,*)'Ifail=',Ifail;
ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.;
-- Variables are given the values such as to guarantee that
-- control passes to the node %GY - to try a new smaller step
Endif;
_FI; -- FOR BADSTE
%G_Cor.L:
/*_Do I=1,NZMOD;
If(IRS(I)==201 & Nstep>=160 & Nstep<=180) then;
WRITE(NCHANL,'(2(A,I5),1P,A,E15.8,A,L3)')
' ROW=',IRS(I),' COL=',ICS(I),' AJAC=',AJAC(I),
' Needbr=',Needbr;
endif
_od; */
write(Nchanl,*)' Nzmod=',Nzmod;
_Do Isw=1,NZMOD;
ASAVE(Isw)=AJAC(Isw);
_od;
--_repeat
If(Needbr)then;
<*SOLY12: *>
IRW=Idisp(11);
IF (IFAIL==2)then;
NUMSING=.true.;
NEEDBR=.true.;
ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.;
Ifail=13;
-- Variables are given the values such as to guarantee that
-- control passes to the node %GY - to try a new smaller step
else;
NUMSING=.false.;
endif;
else;
<*SOLBSF: *>
IRW=IW(1);
IF (IFAIL==6)then;
NUMSING=.true.;
NEEDBR=.true.;
ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.;
Ifail=13;
-- Variables are given the values such as to guarantee that
-- control passes to the node %GY - to try a new smaller step
else;
NUMSING=.false.;
NEEDBR=.false.;
endif;
endif;
IF (NUMSING) THEN;
<*NUMSING: process the numerical singularity in row IR *>;
-- stop 26;
ENDIF;
--_until ^NUMSING ;
%G_Cor.L_SOLY12:
-- IFAIL=110; -- hard failure
IFAIL=111; -- soft failure
/* if(NUMSING)then;
NZMOD=N;
do isw=1,N;
IRN(isw)=isw;
ICN(isw)=isw;
AJAC(isw)=1.d0;
enddo;
else;
*/
_DO Isw=1,NZMOD;
IRN(Isw)=IRS(Isw);
ICN(Isw)=ICS(Isw);
_OD;
/* endif; */
-- Call VTIME(Tm1);
-- CALL F01BRF(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT,
-- IKEEP, IW, WJAC, LBLOCK, GROW, ABORT, IDISP, IFAIL);
-- write(*,*) ' In stiffbgh.trf NZMOD=',NZMOD;
-- pause;
CALL M28Y12(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT,
IKEEP, IW, WJAC, LBLOCK, GROW, ABORT,
IDISP, IFAIL, @THR);
NEEDBR=.FALSE.; -- F01BRF or M28Y12 NOT NEEDED
-- Call VTIME(Tm2);
IF (GROW) WRITE (NCHANL,89996) WJAC(1);
89996:FORMAT (' ON EXIT FROM M28Y12: W(1) = ',1P,E12.4);
WRITE (NCHANL,99996)IDISP(2), IDISP(6),
IDISP(7), IDISP(3),IDISP(4);
IF (LBLOCK) WRITE (NCHANL,99994) (IDISP(I),I=8,10);
99999:FORMAT (6A4, A3);
99998:FORMAT (4(1X/), 1X , 5A4, A3, 'RESULTS'/1X);
99997:FORMAT (5(2X, G8.0, 2X, I1, 2X, I1));
99996:FORMAT (' NUMBER OF NON-ZEROS IN DECOMPOSITION', 9X, '=', I8
/' MINIMUM SIZE OF ARRAY IRN', 20X, '=', I8
/' MINIMUM SIZE OF ARRAYS A AND ICN', 13X, '=', I8
/' NUMBER OF COMPRESSES ON IRN (IDISP(3))', 7X, '=', I8
/' NUMBER OF COMPRESSES ON A AND ICN (IDISP(4)) =', I8);
99994:FORMAT (' STRUCTURAL RANK', 16X, '=', I7
/' NUMBER OF DIAGONAL BLOCKS', 6X, '=', I7
/' SIZE OF LARGEST DIAGONAL BLOCK =', I7);
-- WRITE(NCHANL,
-- '('' M28Y12-1 done at step :'',I6,'' Time(mS):'',I6)')
-- CRAY '('' M28Y12 done at step :'',I6,'' Time(S):'',1P,G11.5)')
-- NSTEP,Tm2-Tm1;
%G_Cor.L_SOLBSF:
IFAIL=11;
-- Call VTIME(Tm1);
CALL M28BYS(N , NZMOD, AJAC, LICN, IRS, ICS, ICN,
IKEEP, IW, WJAC, GROW, ETA, RPMIN, ABORT(4), IDISP, IFAIL);
-- Call VTIME(Tm2);
-- WRITE(NCHANL,'('' M28BYS done at step :'',I6,'' Time(ms):'',I6)')
-- NSTEP,Tm2-Tm1;
/* IF (GROW) WRITE (NCHANL,89995) WJAC(1);
WRITE (NCHANL,89994) RPMIN;
89995:FORMAT (' ON EXIT FROM M28BYS: W(1) = ',1P,E12.4);
89994:FORMAT (' VALUE OF RPMIN = ', G12.4); */
IF(IFAIL^=0 .OR. WJAC(1)>1.D50)THEN;
IF(IFAIL==6)THEN;
<*NUMSING: process the numerical singularity in row IR *>;
ENDIF;
-- IFAIL=110; -- hard failure
IFAIL=111; -- soft failure
_Do Isw=1,NZMOD;
IRN(Isw)=IRS(Isw);
ICN(Isw)=ICS(Isw);
AJAC(Isw)=ASAVE(Isw);
_od;
-- Call VTIME(Tm1);
CALL M28Y12(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT,
IKEEP, IW, WJAC, LBLOCK, GROW, ABORT,
IDISP, IFAIL, @THR);
NEEDBR=.FALSE.; -- F01BRF or M28Y12 NOT NEEDED
-- Call VTIME(Tm2);
/* IF (GROW) WRITE (NCHANL,89996) WJAC(1);
WRITE (NCHANL,99996)IDISP(2), IDISP(6),
IDISP(7), IDISP(3),IDISP(4);
IF (LBLOCK) WRITE (NCHANL,99994) (IDISP(I),I=8,10); */
-- WRITE(NCHANL,
-- '('' M28Y12 done again at step :'',I6,'' Time(ms):'',I6)')
-- CRAY '('' M28Y12 done at step :'',I6,'' Time(S):'',1P,G11.5)')
-- NSTEP,Tm2-Tm1;
ENDIF;
%G_Cor.L_NUMSING:
WRITE(@term,'(3(A,1P,E12.4))')
'Y(IR)=',Y(IRW,1),' Fsave(IR)=',Fsave(IRW),
' Ydot(IR)=',FSAVE(NYDIM+IRW);
_DO I=1,NZMOD;
IF(IRS(I)==IRW)THEN;
WRITE(@term,'(2(A,I5),A,1P,E12.4)') 'IR=',IRS(I),
' IC=',ICS(I),' AJAC=',ASAVE(I);
ICW=ICS(I);
If(ICW<=NZON*NVARS)then;
Izon=MOD(ICW-1,NZON)+1;
Ivar=(ICW-1)/NZON +1;
WRITE(@term,'(2(A,I5))') 'Ivar=',Ivar,' Izon=',Izon;
else;
ICW=ICW-NZON*NVARS;
If(ICW<=KRAD)Then;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@term,'(2(A,I5))') 'FJ in L=',L,' Izon=',Izon;
else;
ICW=ICW-KRAD;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@term,'(2(A,I5))') 'FH in L=',L,' Izon=',Izon;
endif;
endif;
ENDIF;
_OD;
WRITE(@term,'(A)') 'For Ydot of: ';
If(IRW<=NZON*NVARS)then;
Izon=MOD(IRW-1,NZON)+1;
Ivar=(IRW-1)/NZON +1;
WRITE(@term,'(2(A,I5))') 'Ivar=',Ivar,' Izon=',Izon;
else;
IRW=IRW-NZON*NVARS;
If(IRW<=KRAD)Then;
L=(IRW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(IRW-1,(NZON-NCND))+1;
WRITE(@term,'(2(A,I5))') 'FJ in L=',L,' Izon=',Izon;
else;
IRW=IRW-KRAD;
L=(IRW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(IRW-1,(NZON-NCND))+1;
WRITE(@term,'(2(A,I5))') 'FH in L=',L,' Izon=',Izon;
endif;
endif;
%G_Cor.L_SOLBSF_NUMSING:=G_Cor.L_NUMSING:
%G_Cor.U:
<*K: Find solution by M28CYN *>;
If (Ifail^=13) Then;
<*V: Find ERROR, D and put new Y in FSAVE *>;
Endif;
%G_Cor.UK:
_DO I=1,N;
FSAVE(I+NYDIM)=FSAVE(I+NYDIM)*H-Y(I,2)-ERROR(I);
-- R-H SIDE IN FSAVE(NYDIM+1:NYDIM+N)
DB(I)=FSAVE(I+NYDIM); -- TO SAVE R-H SIDE
_OD;
-- ITERATIVE REFINEMENT:
ITQ = 0;
If(Ifail^=13)Then;
_Repeat
CALL M28CYN(N,AJAC,LICN,ICN,IKEEP,FSAVE(NYDIM+1),
WJAC,IW,1,IDISP,RESID,Ifail);
-- SOLUTION IN FSAVE(NYDIM+1:NYDIM+N)
If(Ifail^=13)Then;
ERROLD=ERREST; ERREST=0.D0; XMOD=0.D0;
_DO ID=1,N;
XSAVE(ID)=XSAVE(ID)+FSAVE(ID+NYDIM);
XMOD=XMOD+ABS(XSAVE(ID)); -- NORMA
ERREST=ERREST+ABS(FSAVE(ID+NYDIM));
FSAVE(ID)=0.D0;
_OD;
_DO K=1,NZMOD;
FSAVE(IRS(K))=FSAVE(IRS(K)) -- Polish strategy
+XSAVE(ICS(K))*ASAVE(K);
_OD;
_DO IL=1,N;
FSAVE(IL+NYDIM)=DB(IL)-FSAVE(IL);
_OD;
-- WRITE(NCHANL,'(2(A,1P,G12.3))')' ERROLD=',ERROLD,' ERREST=',ERREST;
ITQ = ITQ+1;
Else;
_DO Isw=1,NZMOD;
IRN(Isw)=IRS(Isw);
ICN(Isw)=ICS(Isw);
AJAC(Isw)=ASAVE(Isw);
_OD;
NEEDBR=.TRUE.; EVALJA=.TRUE.;
@wterm ' G_Cor.UK: ifail 13';
Endif;
-- _Until ERREST<=ANOISE*XMOD ! ITQ > MAXIT;
_Until ERREST<=ANOISE*XMOD .or. Ifail==13
.or. (ITQ>2 & ERREST>ERROLD) .or. ITQ > MAXIT;
IF(ERREST > XMOD*.1d0*EPS) NEEDBR=.TRUE.; -- Esaulov correction
-- IF(XMOD>EPS**3)WRITE(NCHANL,*)' RELEST=',ERREST/XMOD;
Endif;
%G_Cor.UV:
D=0.D0;
_DO I=1,N;
ERROR(I)=ERROR(I)+XSAVE(I);
D=D+(XSAVE(I)/YMAX(I))**2;
FSAVE(I)=Y(I,1)+EL(1)*ERROR(I)
_OD;
%GX:
D=0.D0;
ERMAX=0.D0; MAXER=0;
do I=1,N;
IF((ERROR(I)/YMAX(I))**2>ERMAX)THEN;
ERMAX=(ERROR(I)/YMAX(I))**2; MAXER=I;
ENDIF;
D=D+(ERROR(I)/YMAX(I))**2;
enddo;
OLDJAC=.TRUE.;
if( D<=E) then;
-- After a succesful step, update the Y array and YMAX.
-- Consider changing H if IDOUB==1. Otherwise decrease
-- IDOUB by 1.
-- If a change in H is considered, an increase
-- or decrease in order by one is considered also.
HUSED=H;NQUSED=NQ;KFLAG=0;
IF(METH==3)THEN; -- Save time intervals in HNT
do J=7,3,-1;
HNT(J)=HNT(J-1) + H;
-- HNT(LNQ+1) saves old H for possible order increase
-- Initially HNT(2:6) is undefined here!
enddo;
HNT(2)=H;
ENDIF;
do J=1,LNQ;
do I=1,N;
Y(I,J)=Y(I,J)+EL(J)*ERROR(I);
enddo; enddo;
if ( IDOUB==1 ) then;
<*A: Factors PR1, PR2 and PR3 are computed, by which H
could be multiplied at order NQ-1, order NQ,
or order NQ+1, respectively. The largest of these
is determined and the new order is chosen accordingly.
If the order is to be increased, we compute one
additional scaled derivative. *>
RETCOD=3; -- EXIT
else;
IDOUB=IDOUB-1;
IF((IDOUB==1) & (NQ^=MAXORD)) THEN;
-- Error is saved for use in a possible
-- order increase on the next step.
-- _DO I=1,N; Y(I,LMAX)=ERROR(I) _OD;
Y(1:N,LMAX)=ERROR(1:N);
ENDIF;
IF(METH==3)THEN;
<*TWO: RESET PR2 *>;
RH=PR2;
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
-- NO IDOUB HERE!!!
ENDIF;
RETCOD=4; -- EXIT
endif;
else;
-- The error test failed. KFLAG keeps track
-- of multiple failures. Restore t and the Y array to their
-- previous values, and prepare to try the step again.
-- Compute the optimum step size for this or one lower order.
NFAIL=NFAIL+1;KFLAG=KFLAG-1;t=tOLD;
do J1=1,NQ;
do J2=J1,NQ;
J=NQ-J2+J1;
do I=1,N;
Y(I,J)=Y(I,J)-Y(I,J+1);
enddo;enddo;enddo;
RMAX=RMXFAI;
if ( ABS(H)<=(HMIN*1.00001D0) ) then;
RETCOD=1;
else;
if( KFLAG<=-MAXFAI) then;
<*Z: Control reaches this section if MAXFAI or more
failures have occured. It is assumed that the
derivatives that have accumulated in the Y array
have errors of the wrong order.
Hence the first derivative is recomputed, and the order
is set to one. Then H is reduced by factor of RHERR3,
and the step is retried. *>
else;
<*2: RESET PR2 *>;
if(NQ==1) then;
RH=PR2;
CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
else;
<*1: RESET PR1 *>;
if(PR1<=PR2) then;
RH=PR2;
if(meth==3)then;
HNT(1)=H*RH;
CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);
RC=RC*EL(1)/OLDL0;OLDL0=EL(1);
<*E: the constants E,EUP,EDN,BND *>;
CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ);
else;
CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
endif;
else;
<*R: CHANGE NQ,RESET CONST. *>;
endif;
endif;
endif;
endif;
endif;
%GXE=A:
%GX1:
SUM=0.D0;
_DO I=1,N;SUM=SUM+(Y(I,LNQ)/YMAX(I))**2 _OD;
PR1=1.D0/(((SUM/EDN)**(.5D0/DBLE(NQ)))*BIAS1+DFLTZR);
%GX2:
PR2=1.D0/(((D/E)**(.5D0/DBLE(LNQ)))*BIAS2+DFLTZR);
%GX_TWO=GX2:
%GXR:
LNQ=NQ;NQ=NQ-1;RH=PR1;
HNT(1)=H*RH;
CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);
RC=RC*EL(1)/OLDL0;OLDL0=EL(1);
<*A: EDN,E,EUP,BND *>;
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
%GXRA=A:
%GXZ:
RH=RHERR3;RH=MAX(HMIN/ABS(H),RH);H=H*RH;
HL=H*EL(1);
_DO I=1,N; FSAVE(I)=Y(I,1) _OD;
CALL DIFMAT(envs, err);
_DO I=1,N; Y(I,2)=FSAVE(NYDIM+I)*H _OD;
-- EVALJA=.TRUE.;
EVALJA=.False.; -- to get out the _repeat loop
RC=0.d0; -- to force a new Jacobian
IDOUB=IDELAY;
IF(NQ^=1) THEN;
NQ=1;LNQ=2;
CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);
OLDL0=EL(1);
<*A: EDN,E,EUP,BND *>;
ENDIF;
%GXZA=A:
%GXA:
PR3=DFLTZR;
IF(NQ^=MAXORD) THEN;
<*3: RESET PR3 *>;
ENDIF;
<*2: RESET PR2 *>;
PR1=DFLTZR;
IF(NQ^=1) THEN;<*1: RESET PR1 *>;ENDIF;
IF((PR2>PR1) & (PR2>PR3))THEN;
RH=PR2;
ELSE;
/*
IF(PR1>PR3) THEN;
NEWQ=NQ-1;RH=PR1;
ELSE;
NEWQ=LNQ;RH=PR3; -- PR3 IS GREATER
_DO I=1,N;Y(I,NEWQ+1)=ERROR(I)*EL(LNQ)/DBLE(LNQ)_OD;
ENDIF;
*/
-- bug (for METH==3) in above lines corrected:
IF(PR1>PR3) THEN;
NEWQ=NQ-1;RH=PR1;
ELSEIF(NQ^=MAXORD)THEN;
NEWQ=LNQ;RH=PR3; -- PR3 IS GREATER
_DO I=1,N;Y(I,NEWQ+1)=ERROR(I)*EL(LNQ)/DBLE(LNQ) _OD;
ELSE; -- NQ==MAXORD
NEWQ=NQ;RH=PR3; -- PR3 IS GREATER
ENDIF;
-- If there is a change of order, reset NQ and LNQ.
-- In any case H is reset according to
-- RH and the Y array and the coefficients EL()
-- are rescaled. Then exit.
NQ=NEWQ;LNQ=NQ+1;
ENDIF;
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
%GXA3:
SUM=0.D0;
IF(METH^=3)THEN;
_DO I=1,N;SUM=SUM+((ERROR(I)-Y(I,LMAX))/YMAX(I))**2 _OD;
PR3=1.D0/(((SUM/EUP)**(.5D0/DBLE(LNQ+1)))*BIAS3+DFLTZR);
ELSE;
_DO I=1,N;
SUM=SUM+
((ERROR(I)*EL(1)/DBLE(NQ+2)-Y(I,LMAX)/TQ(3))/YMAX(I))**2
_OD;
PR3=1.D0/(((SUM/EPSJ**2)**(.5D0/DBLE(LNQ+1)))*BIAS3+DFLTZR);
ENDIF;
%GXA2=GX2:
%GXA1=GX1:
%GY:
t=tOLD;RMAX=RMXFAI;
_DO J1=1,NQ;
_DO J2=J1,NQ;
J=NQ-J2+J1;
_DO I=1,N;
Y(I,J)=Y(I,J)-Y(I,J+1)
_OD_OD_OD;
IF(ABS(H)<=(HMIN*1.00001D0)) THEN;
RETCOD=2;
ELSE;
write(*,'(a,i7,i3,1p,2e12.3)') ' GY: maxer, kflag, H, Hu:',
maxer, kflag, H, Hu;
RH=RHCORR;
IF(METH==3)THEN;
HNT(1)=H*RH;
CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);
<*A: the constants E,EUP,EDN,BND *>;
RC=RC*EL(1)/OLDL0;OLDL0=EL(1);
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
ELSE;
CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ);
IDOUB=LNQ+1;
ENDIF;
ENDIF;
if(ifail==13)then;
@wterm ' GY: ifail 13';
endif;
%GYA=A:
%G_COR_MatrOut:
_DO JTST=1,N; --@VAR
<*ROW: put in IRN(NZ+1:ILAST) all indeces I such that
ICN(I)==JTST *>:
_Do KTST=1,N; --@NAMELEFT,@NAMERIGHT;
<*COMPR: IF in IRN(NZ+1:ILAST) there is I such that
IRN(I)==KTST THEN
output relevant vars *>
_od;
_OD
%G_COR_MatrOut_ROW:
ILAST=0;
do I=1,NZMOD;
IF(ICS(I)==JTST)THEN;
ILAST=ILAST+1;
IRN(NZ+ILAST)=I;
ENDIF;
enddo;
-- WRITE(5,*)'ILAST=',ILAST;
%G_COR_MatrOut_COMPR:
FOUNDI=.FALSE.; IL=1;
_WHILE IL<=ILAST & ^FOUNDI _DO
IF(IRS(IRN(NZ+IL))==KTST)THEN;
FOUNDI=.TRUE.;
ELSE;
IL=IL+1;
ENDIF;
_OD;
IF(FOUNDI)THEN;
If(KTST^=JTST)then;
if(AJAC(IRN(NZ+IL))>1.d-6) WRITE(@Wres,'(A,2I7,1P,3E12.3)')
' K,J,AJAC : ',
KTST,JTST,AJAC(IRN(NZ+IL));
else;
WRITE(@Wres,'(A,2I7,1P,3E12.3)')
' K,J,AJAC,1-AJAC : ',
KTST,JTST,AJAC(IRN(NZ+IL)),1.d0-AJAC(IRN(NZ+IL));
endif;
ENDIF;
<*PhysInfo : give info on variable number KTST *>;
%G_COR_MatrOut_COMPR_PhysInfo:
do I=1,NZMOD;
if(IRS(I)==KTST .and. ICS(I)==JTST)then;
WRITE(@Wres,'(3(A,1P,E12.4))')
' Y(K)=',Y(KTST,1),' Fsave(K)=',Fsave(KTST),
' Ydot(K)=',FSAVE(NYDIM+KTST);
WRITE(@Wres,'(2(A,I5),A,1P,E12.4)') ' Control IR=',KTST,
' IC=',JTST,' AJAC=',AJAC(I);
endif;
enddo; -- this is for control of found AJAC;
ICW=JTST;
If(ICW<=NZON*NVARS)then;
Izon=MOD(ICW-1,NZON)+1;
Ivar=(ICW-1)/NZON +1;
WRITE(@Wres,'(3A,I3)')
' Derivative over ',CharVarName(Ivar)(1:length(CharVarName(Ivar))),
' in zon=',Izon;
else;
ICW=ICW-NZON*NVARS;
If(ICW<=KRAD)Then;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I3))') ' over FJ in L=',L,' Izon=',Izon;
else;
ICW=ICW-KRAD;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I3))') ' over FH in L=',L,' Izon=',Izon;
endif;
endif;
If(KTST<=NZON*NVARS)then;
Izon=MOD(KTST-1,NZON)+1;
Ivar=(KTST-1)/NZON +1;
WRITE(@Wres,'(3A,I3)')
' for Ydot of ',CharVarName(Ivar)(1:length(CharVarName(Ivar))),
' in zon=',Izon;
else;
KTST1=KTST-NZON*NVARS;
If(KTST1<=KRAD)Then;
L=(KTST1-1)/(NZON-NCND)+1;
Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FJ) in L=',L,' Izon=',Izon;
else;
KTST1=KTST1-KRAD;
L=(KTST1-1)/(NZON-NCND)+1;
Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FH) in L=',L,' Izon=',Izon;
endif;
endif;
%G_COR_DEB:
_DEFINE
@NOISE 1.D-4 -- for IBM Double Precision 3.D-8
-- 1.D-7 -- for Cray
-- input variables:
-- @NAME - zero starting position in FSAVE-array for testing
-- @VAR - zero starting position in FSAVE-array for indep.var
-- @NOISE- ignore derivatives less than @NOISE
If(^BADSTE)Then;
t=tOLD;
_DO J1=1,NQ;
_DO J2=J1,NQ;
J=NQ-J2+J1;
_DO I=1,N;
Y(I,J)=Y(I,J)-Y(I,J+1)
_OD_OD_OD;
_Do Ich=1,Mzon;
Chem0(Ich)=0.1; -- check if we need this
_od;
_Do I=1,N; FSAVE(I)=Y(I,1)_od;
--************** ************** ************** **************
-- EVALJA=.FALSE.; -- never put it here: some errors may be lost
--************** ************** ************** **************
Nperturb=0; -- index of perturbation, now NO perturbation
CALL DIFMAT(envs, err);
_Do KTST=1,N; -- save old YDOT
DSAVE(KTST)=FSAVE(NYDIM+KTST);
_od;
_DO JTST=1,N; --@VAR
-- _DO JTST=Ncnd+1,Ncnd+2; -- only d over dr is tested
@wterm' JTST=',JTST;
EVALJA=.FALSE.;
OLD=FSAVE(JTST);
If(JTST<Nzon)then; -- radius
If(JTST==1)then;
DeltY=DeltaDeb*min(FSAVE(JTST+1)-OLD,OLD-Rce);
else;
DeltY=DeltaDeb*min(FSAVE(JTST+1)-OLD,OLD-FSAVE(JTST-1));
endif;
elseIf(Nzon<JTST & JTST<2*Nzon)then; -- velocity
If(JTST==Nzon+1)then;
DeltY=DeltaDeb*min(abs(FSAVE(JTST+1)-OLD),abs(OLD));
else;
DeltY=DeltaDeb*min(abs(FSAVE(JTST+1)-OLD),
abs(OLD-FSAVE(JTST-1)));
endif;
else;
DeltY=OLD*DeltaDeb;
endif;
FSAVE(JTST)=OLD+DeltY;
_Do Ich=1,Mzon;
Chem0(Ich)=0.1
_od;
Nperturb=1; -- index of perturbation, now IS perturbation
CALL DIFMAT(envs, err);
FSAVE(JTST)=OLD;
<*ROW: put in IRN(NZ+1:ILAST) all indeces I such that
ICN(I)==JTST *>;
_Do KTST=1,N; --@NAMELEFT,@NAMERIGHT;
@DNAMEDVAR(KTST) = (FSAVE(NYDIM+KTST)-DSAVE(KTST));
<*COMPR: IF in IRN(NZ+1:ILAST) there is I such that
IRN(I)==KTST THEN
IF AJAC(I) differs strongly
from -HL*@DNAMEDVAR for KTST^=JTST or
from 1.-HL*@DNAMEDVAR for KTST==JTST
THEN output relevant vars
ELSE
IF @DNAMEDVAR is HIGHER than @NOISE
THEN output relevant vars *>
_od;
_OD;
STOP;
else;
WRITE(@Wres,*)' Testing not possible: BADSTE=',BADSTE;
endif;
%G_COR_DEB_ROW:
ILAST=0;
_DO I=1,NZMOD;
IF(ICS(I)==JTST)THEN;
ILAST=ILAST+1; IRN(NZ+ILAST)=I;
ENDIF;
_OD;
-- WRITE(5,*)'ILAST=',ILAST;
%G_COR_DEB_COMPR:
FOUNDI=.FALSE.; IL=1;
_WHILE IL<=ILAST & ^FOUNDI _DO
IF(IRS(IRN(NZ+IL))==KTST)THEN;
FOUNDI=.TRUE.;
ELSE;
IL=IL+1;
ENDIF;
_OD;
IF(FOUNDI)THEN;
If(KTST^=JTST)THEN;
If( ABS(AJAC(IRN(NZ+IL))*DeltY
+HL*@DNAMEDVAR(KTST)) >
1.D-1*ABS(AJAC(IRN(NZ+IL))*DeltY)
& ABS(AJAC(IRN(NZ+IL))/HL) > @NOISE ) Then;
If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then;
WRITE(@Wres,'(/A/2I7,1P,3G12.3)')
' K,J,AJAC,Stella df/dY,Num DN(K)/DV(J): ',
KTST,JTST,AJAC(IRN(NZ+IL)),
-AJAC(IRN(NZ+IL))/HL,@DNAMEDVAR(KTST)/DeltY;
else;
WRITE(@Wres,'(/A/2I7,1P,3G12.3)')
' K,J,AJAC,(df/dY,DN/DV)*DeltaDeb: ',
KTST,JTST,AJAC(IRN(NZ+IL)),
-AJAC(IRN(NZ+IL))*DeltY/HL,@DNAMEDVAR(KTST);
endif;
<*PhysInfo: give info on variable number KTST *>;
endif;
else;
If( ABS((AJAC(IRN(NZ+IL))-1.D0)*DeltY+HL*@DNAMEDVAR(KTST)) >
1.D-1*ABS((AJAC(IRN(NZ+IL))-1.D0)*DeltY)
& ABS((AJAC(IRN(NZ+IL))-1.D0)/HL) > @NOISE ) Then;
If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then;
WRITE(@Wres,'(/A/2I7,1P,3G12.3)')
' K,J,AJAC,Stella df/dY,Num DN(K)/DV(J): ',
KTST,JTST,AJAC(IRN(NZ+IL)),(1.D0 -AJAC(IRN(NZ+IL)) )/HL,
@DNAMEDVAR(KTST)/DeltY ;
else;
WRITE(@Wres,'(/A/2I7,1P,3G12.3)')
' K,J,AJAC,(df/dY,DN/DV)*DeltaDeb: ',
KTST,JTST,AJAC(IRN(NZ+IL)),(1.D0-AJAC(IRN(NZ+IL)))*DeltY/HL,
@DNAMEDVAR(KTST);
endif;
<*PhysInfo2: give info on variable number KTST *>;
endif;
endif;
ELSE;
IF(ABS(@DNAMEDVAR(KTST)) > @NOISE) then;
If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then;
WRITE(@Wres,'(/A,2I7,1P,2G12.3)')' NO JAC: K,J,DN(K)/DV(J) ',
KTST,JTST, @DNAMEDVAR(KTST)/DeltY;
else;
WRITE(@Wres,'(/A,2I7,1P,2G12.3)')' NO JAC: K,J,DN(K)/DV(J)*DeltaDeb ',
KTST,JTST, @DNAMEDVAR(KTST);
endif;
endif;
ENDIF;
%G_COR_DEB_COMPR_PhysInfo:
do I=1,NZMOD;
if(IRS(I)==KTST .and. ICS(I)==JTST)then;
WRITE(@Wres,'(3(A,1P,E12.4))')
' Y(K)=',Y(KTST,1),' Fsave(K)=',Fsave(KTST),
' Ydot(K)=',FSAVE(NYDIM+KTST);
WRITE(@Wres,'(2(A,I5),A,1P,E12.4)') ' Control IR=',KTST,
' IC=',JTST,' AJAC=',AJAC(I);
endif;
enddo; -- this is for control of found AJAC;
ICW=JTST;
If(ICW<=NZON*NVARS)then;
Izon=MOD(ICW-1,NZON)+1;
Ivar=(ICW-1)/NZON +1;
WRITE(@Wres,'(3A,I3)')
' Derivative over ',CharVarName(Ivar)(1:length(CharVarName(Ivar))),
' in zon=',Izon;
else;
ICW=ICW-NZON*NVARS;
If(ICW<=KRAD)Then;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I3))') ' over FJ in L=',L,' Izon=',Izon;
else;
ICW=ICW-KRAD;
L=(ICW-1)/(NZON-NCND)+1;
Izon=NCND+MOD(ICW-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I3))') ' over FH in L=',L,' Izon=',Izon;
endif;
endif;
If(KTST<=NZON*NVARS)then;
Izon=MOD(KTST-1,NZON)+1;
Ivar=(KTST-1)/NZON +1;
WRITE(@Wres,'(3A,I3)')
' for Ydot of ',CharVarName(Ivar)(1:length(CharVarName(Ivar))),
' in zon=',Izon;
else;
KTST1=KTST-NZON*NVARS;
If(KTST1<=KRAD)Then;
L=(KTST1-1)/(NZON-NCND)+1;
Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FJ) in L=',L,' Izon=',Izon;
else;
KTST1=KTST1-KRAD;
L=(KTST1-1)/(NZON-NCND)+1;
Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1;
WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FH) in L=',L,' Izon=',Izon;
endif;
endif;
%G_COR_DEB_COMPR_PhysInfo2=G_COR_DEB_COMPR_PhysInfo:
%X:
-- _Case RETCOD _Of
-- -- Control passes to _Esac for RETCOD==4
-- _1 KFLAG=-1
-- _2 KFLAG=-2
-- _3 RMAX=RMXNOR
-- _Esac;
select case(RETCOD);
case(1);
KFLAG=-1;
case(2);
KFLAG=-2;
case(3);
RMAX=RMXNOR;
end select;
HOLD=H;
JSTART=NQ;
English
Русский