"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"