code 33170;
  procedure RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,K,
  RATECONV,DOMEIGVAL,OUT); value LJ,UJ,LL,UL,A,B;
  integer N,K,LJ,UJ,LL,UL; real A,B,RATECONV,DOMEIGVAL; boolean 
  INAP; array U,DISCR; procedure RESIDUAL,OUT;
  begin integer J,L; real X,Y,Z,Y0,C,D,ALFA,OMEGA,OMEGA0,
    EIGMAX,EIGEUCL,EUCLRES,MAXRES,RCMAX,RCEUCL,MAXRES0,EUCLRES0;
    array V,RES[LJ:UJ,LL:UL];
    procedure CALPAR;
    comment CALPAR CALCULATES THE PARAMETERS ALFA AND OMEGA FOR
    EACH ITERATION;
    begin ALFA:= Z/(Z - ALFA);
      OMEGA:= 1/(X - OMEGA * Y)
    end CALPAR;
    procedure ITERATION;
    comment FIRST THE ITERATION FORMULA IS CONSTRUCTED;
    begin real AUXV,AUXU,AUXRES,EUCLUV,MAXUV;
      EUCLUV:= EUCLRES:= MAXUV:= MAXRES:= 0;
      for J:= LJ step 1 until UJ do 
      for L:= LL step 1 until UL do RES[J,L]:= V[J,L];
      RESIDUAL(RES);
      for J:= LJ step 1 until UJ do 
      for L:= LL step 1 until UL do 
      begin AUXV:= U[J,L]; AUXU:= V[J,L]; AUXRES:= RES[J,L];
        AUXV:= ALFA * AUXU - OMEGA * AUXRES + (1 - ALFA) * AUXV;
        V[J,L]:= AUXV; U[J,L]:= AUXU;
        comment THE NORMS OF THE K-TH RESIDUAL AND THE DIFFERENCE
        BETWEEN THE (K+1)-TH AND K-TH ITERAND ARE CALCULATED;
        AUXU:= ABS(AUXU - AUXV); AUXRES:= ABS(AUXRES);
        MAXUV:= if MAXUV < AUXU then AUXU else MAXUV;
        MAXRES:= if MAXRES < AUXRES then AUXRES else MAXRES;
        EUCLUV:= EUCLUV + AUXU * AUXU;
        EUCLRES:= EUCLRES + AUXRES * AUXRES;
      end;
      EUCLUV:= SQRT(EUCLUV); EUCLRES:= SQRT(EUCLRES);
      DISCR[1]:= EUCLRES; DISCR[2]:= MAXRES;
      comment DOMEIGVAL IS EVALUATED;
      MAXUV:= MAXRES/MAXUV; EUCLUV:= EUCLRES/EUCLUV;
      EIGMAX:= MAXUV * (C - MAXUV)/(.25 * D - MAXUV);
      EIGEUCL:= EUCLUV * (C - EUCLUV)/(.25 * D - EUCLUV);
      DOMEIGVAL:= .5 * (EIGMAX + EIGEUCL);
      comment FINALLY THE RATE OF CONVERGENCE IS CALCULATED;
      RCEUCL:= -LN(EUCLRES/EUCLRES0)/K;
      RCMAX:= -LN(MAXRES/MAXRES0)/K;
      RATECONV:= .5 * (RCEUCL + RCMAX)
    comment THE CONSTANTS FOR STARTING CALPAR ARE CALCULATED;
    ALFA:= 2; OMEGA:= 4/(B + A); Y0:= (B + A)/(B - A);
    X:= .5 * (B + A); Y:= (B - A) * (B - A)/16; Z:= 4 * Y0 * Y0;
    comment THE CONSTANTS NEEDED FOR DOMEIGVAL ARE CALCULATED;
    C:= A * B; C:= SQRT(C); D:= SQRT(A) + SQRT(B); D:= D * D;
    comment THE INITIAL APPROXIMATION IS PUT INTO ARRAY U;
    if ^INAP then 
    begin for J:= LJ step 1 until UJ do 
      for L:= LL step 1 until UL do U[J,L]:= 1
    end;
    comment THE ZEROTH ITERATION IS NOW PERFORMED;
    K:= 0;
    for J:= LJ step 1 until UJ do 
    for L:= LL step 1 until UL do RES[J,L]:= U[J,L];
    RESIDUAL(RES);
    OMEGA0:= 2/(B+A);
    begin real AUXRES0;
      MAXRES0:= EUCLRES0:= 0;
      for J:= LJ step 1 until UJ do 
      for L:= LL step 1 until UL do 
      begin AUXRES0:= RES[J,L];
        V[J,L]:= U[J,L] - OMEGA0 * AUXRES0;
        AUXRES0:= ABS(AUXRES0);
        MAXRES0:= if MAXRES0 < AUXRES0 then AUXRES0 else MAXRES0;
        EUCLRES0:= EUCLRES0 + AUXRES0 * AUXRES0
      end;
    EUCLRES0:= SQRT(EUCLRES0)
    end;
    DISCR[1]:= EUCLRES0; DISCR[2]:= MAXRES0;
    OUT(K);
    if K >= N then goto FINALLY;
  NEXT STEP:
    K:= K + 1; CALPAR; ITERATION; OUT(K);
    if K < N then goto NEXT STEP;
  FINALLY:
  end RICHARDSON

        eop