begin;
    comment  These examples were described in [D.E.Knuth, J.N.Merner.
             Algol 60 Confidential, Comm.A.C.M., 4, 6 (1961), pp.268-272].
             The authors show possibilities of recursive procedures and
             call-by-name parameters that based on examples of usage of
             the procedure GPS (General Problem Solver) in some assignment
             statements. One of those statements performing matrix multi-
             plication was offered as a challenge to Algol 60 translator
             developers
    ;
    real procedure GPS(I, N, Z, V);
      real I, N, Z, V;
    begin;
        for I := 1 step 1 until N do Z := V;
        GPS := 1;
    end ;
    comment  The first example demonstrates usage of the procedure
             GPS to create a matrix A[i,j] = i+j
    ;
    FIRSTEXAMPLE: begin;
        real I, J;
        array A[1 : 2, 1 : 3];
        OUTSTRING(1, "First example\n");
        I := GPS(J, 3.0, I, GPS(I, 2.0, A[I, J], I + J));
        OUTSTRING(1, "Matrix A:\n");
        OUTREAL(1, A[1, 1]);
        OUTREAL(1, A[1, 2]);
        OUTREAL(1, A[1, 3]);
        OUTSTRING(1, "\n");
        OUTREAL(1, A[2, 1]);
        OUTREAL(1, A[2, 2]);
        OUTREAL(1, A[2, 3]);
        OUTSTRING(1, "\n");
    end OFFIRSTEXAMPLE;
    comment  The second example demonstrates usage of the procedure
             GPS to perform matrix multiplication
    ;
    SECONDEXAMPLE: begin;
        procedure TESTMATR(N, A); 
          value N;
          integer N;
          array A;
        comment  Create test matrix (CACM, Algorithm 52);
        begin;
            real C, D;
            integer I, J;
            C := N  (N + 1)  (2  N - 5) ÷ 6;
            D := 1 ÷ C;
            A[N, N] := -D;
            for I := 1 step 1 until N - 1 do
              begin;
                A[I, N] := A[N, I] := D  I;
                A[I, I] := D  (C - I  2);
                for J := 1 step 1 until I - 1 do
                  A[I, J] := A[J, I] := -D  I  J;
            end ;
        end TESTMATR;
        procedure INVERT140(N, EPS, OUT, A); 
          value N, EPS;
          real EPS;
          integer N;
          array A;
          label OUT;
        comment  Invert matrix (CACM, Algorithm 140);
        begin;
            real Q;
            integer I, J, K;
            for I := 1 step 1 until N do
              begin;
                if ABS(A[I, I])  EPS then goto OUT;
                Q := 1 ÷ A[I, I];
                A[I, I] := 1;
                for K := 1 step 1 until N do
                  A[I, K] := A[I, K]  Q;
                for J := 1 step 1 until N do
                  if I  J then begin;
                    Q := A[J, I];
                    A[J, I] := 0;
                    for K := 1 step 1 until N do
                      A[J, K] := A[J, K] - Q  A[I, K];
                end J;
            end I;
        end INVERT140;
        procedure PRINTMATRIX(NAME, N, A); 
          value N;
          string NAME;
          integer N;
          array A;
        begin;
            integer I, J;
            OUTSTRING(1, "Matrix ");
            OUTSTRING(1, NAME);
            OUTSTRING(1, ":\n");
            for I := 1 step 1 until N do
              begin;
                for J := 1 step 1 until N do
                  OUTREAL(1, if ABS(A[I, J]) < 10-12 then 0 else A[I, J]);
                OUTCHAR(1, "\n", 1);
            end I;
        end PRINTMATRIX;
        comment  N IS ORDER OF MATRICES
        ;
        integer N;
        N := 5;
        OUTSTRING(1, "Second example\n");
        begin;
            array A, B, C[1 : N, 1 : N];
            integer I, J, K;
            comment  Create test matrix A
            ;
            TESTMATR(N, A);
            PRINTMATRIX("A", N, A);
            comment  b := INV(a)
            ;
            TESTMATR(N, B);
            INVERT140(N, EPSILON, SING, B);
            goto SKIP;
            SING: FAULT("Matrix is singular", N);
            SKIP: PRINTMATRIX("B = INV(A)", N, B);
            comment  C := A * B using GPS
            ;
            I := GPS(I, 1.0, C[1, 1], 0.0)  GPS(I, (N - 1)  GPS(J, (N - 1)  GPS(K, N, C[I, J], C[I, J] + A[I, K]  B[K, J]), C[I, J + 1], 0.0), C[I + 1, 1], 0.0);
            PRINTMATRIX("C = A * B", N, C);
        end ;
    end OFSECONDEXAMPLE;
end ;