code 33132;
procedure LINIGER1VS(X,XE,M,Y,SIGMA,DERIVATIVE,J,
   JACOBIAN,ITMAX,HMIN,HMAX,AETA,RETA,INFO,OUTPUT);
value M;
integer M,ITMAX;
real X,XE,SIGMA,AETA,RETA,HMIN,HMAX;
array Y,J,INFO;
procedure DERIVATIVE,JACOBIAN,OUTPUT;

begin integer I,ST,LASTJAC;
    real H,HNEW,MU,MU1,BETA,P,E,E1,ETA,ETA1,DISCR;
    boolean LAST,FIRST,EVALJAC,EVALCOEF;
    integer array PI[1:M];
    real array DY,YL,YR,F[1:M],A[1:M,1:M],AUX[1:3];

    real procedure NORM(A); array A;
    NORM:=SQRT(VECVEC(1,M,0,A,A));

    procedure COEFFICIENT;
    begin real B,E; B:=ABS(H*SIGMA);
        if B>40 then 
        begin MU:=1/B; BETA:=1; P:=2+2/(B-2)
        end else 
        if B<.04 then 
        begin E:=B*B/30; P:=3-E;
            MU:=.5-B/12*(1-E/2);
            BETA:=.5+B/6*(1-E)
        end else 
        begin E:=EXP(B)-1;
            MU:=1/B-1/E;
            BETA:=(1-B/E)*(1+1/E);
            P:=(BETA-MU)/(.5-MU)
        end;
        MU1:=H*(1-MU);
        LUDECOMP
    procedure LUDECOMP;
    begin integer I;
        for I:=1 step 1 until M do 
        begin MULROW(1,M,I,I,A,J,-MU1);
            A[I,I]:=A[I,I]+1
        end;
        AUX[2]:=0; DEC(A,M,AUX,PI)
    end LUDECOMP;

    procedure STEPSIZE;
    if ETA<0 then 
    begin real HL; HL:=H;
        H:=HNEW:=HMAX; INFO[5]:=INFO[5]+1;
        if 1.1*HNEW>XE-X then 
        begin LAST:=true; H:=HNEW:=XE-X
        end;
        EVALCOEF:=H^=HL;
    end else 
    if FIRST then 
    begin H:=HNEW:=HMIN; FIRST:=false; INFO[4]:=INFO[4]+1
    end else 
    begin real B,HL;
        B:=DISCR/ETA; HL:=H; if B<.01 then B:=.01;
        HNEW:= if B>0 then H*B**(-1/P) else HMAX;
        if HNEW<HMIN then 
        begin HNEW:=HMIN; INFO[4]:=INFO[4]+1
        end else 
        if HNEW>HMAX then 
        begin HNEW:=HMAX; INFO[5]:=INFO[5]+1 end;
        if 1.1*HNEW>=XE-X then 
        begin LAST:=true; H:=HNEW:=XE-X
        end else 
        if ABS(H/HNEW-1)>.1 then H:=HNEW;
        EVALCOEF:=H^=HL
    end STEPSIZE;

    procedure LINEARITY(ERROR); real ERROR;
    begin integer K;
        for K:=1 step 1 until M do 
        DY[K]:=Y[K]-MU1*F[K];
        SOL(A,M,PI,DY);
        ELMVEC(1,M,0,DY,Y,-1);
        ERROR:=NORM(DY)
    procedure ITERATION(I); integer I;
    if RETA<0 then 
    begin integer K;
        if I=1 then 
        begin MULVEC(1,M,0,DY,F,H);
            for K:=1 step 1 until M do YL[K]:=Y[K]+MU*DY[K];
            SOL(A,M,PI,DY); E:=1;
        end else 
        begin for K:=1 step 1 until M do 
            DY[K]:=YL[K]-Y[K]+MU1*F[K];
            if E*NORM(Y)>E1*E1 then 
            begin EVALJAC:=I>=3;
                if I>3 then 
                begin INFO[3]:=INFO[3]+1; JACOBIAN(J,Y);
                    LUDECOMP
                end;
            end;
            SOL(A,M,PI,DY)
        end;
        E1:=E; E:=NORM(DY);
        ELMVEC(1,M,0,Y,DY,1);
        ETA:=NORM(Y)*RETA+AETA;
        DISCR:=0;
        DUPVEC(1,M,0,F,Y);
        DERIVATIVE(F);
        INFO[2]:=INFO[2]+1;
    end else 
    begin integer K;
        if I=1 then 
        begin LINEARITY(E);
            if E*(ST-LASTJAC)>ETA then 
            begin JACOBIAN(J,Y); LASTJAC:=ST;
                INFO[3]:=INFO[3]+1;
                H:=HNEW; COEFFICIENT;
                LINEARITY(E)
            end;
            EVALJAC:= E*(ST+1-LASTJAC)>ETA;
            MULVEC(1,M,0,DY,F,H);
            for K:=1 step 1 until M do YL[K]:=Y[K]+MU*DY[K];
            SOL(A,M,PI,DY);
            for K:=1 step 1 until M do 
            YR[K]:=H*BETA*MATVEC(1,M,K,J,DY);
            SOL(A,M,PI,YR);
            ELMVEC(1,M,0,YR,DY,1);
        end else 
        begin for K:=1 step 1 until M do 
            DY[K]:=YL[K]-Y[K]+MU1*F[K];
            if E>ETA1 and DISCR>ETA1 then 
            begin INFO[3]:=INFO[3]+1; JACOBIAN(J,Y);
                LUDECOMP
            end;
            SOL(A,M,PI,DY);
            E:=NORM(DY)
        end;
        ELMVEC(1,M,0,Y,DY,1);
        ETA:=NORM(Y)*RETA+AETA;
        ETA1:=ETA/SQRT(RETA);
        DUPVEC(1,M,0,F,Y);
        DERIVATIVE(F);
        INFO[2]:=INFO[2]+1;
        for K:=1 step 1 until M do DY[K]:=YR[K]-H*F[K];
        DISCR:=NORM(DY)/2
    end ITERATION;

    FIRST:=EVALJAC:=true; LAST:=EVALCOEF:=false;
    INIVEC(1,9,INFO,0);
    ETA:=RETA*NORM(Y)+AETA;
    ETA1:=ETA/SQRT(ABS(RETA));
    DUPVEC(1,M,0,F,Y);
    DERIVATIVE(F);
    INFO[2]:=1;
    for ST:=1,ST+1 while ^LAST do 
    begin STEPSIZE; INFO[1]:=INFO[1]+1;
        if EVALJAC then 
        begin JACOBIAN(J,Y);
            INFO[3]:=INFO[3]+1;
            H:=HNEW;
            COEFFICIENT;
            EVALJAC:=false; LASTJAC:=ST
        end else 
        if EVALCOEF then COEFFICIENT;
        for I:=1,I+1 while E>ABS(ETA) and DISCR>1.3*ETA
        and I<=ITMAX do 
        begin ITERATION(I); if I>INFO[6] then INFO[6]:=I;
        end; INFO[7]:=ETA; INFO[8]:=DISCR;
        X:=X+H;
        if DISCR>INFO[9] then INFO[9]:=DISCR;
        OUTPUT;
    end;
end LINIGER1VS;
        eop