1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR: M. BAKKER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
 
 
 RECEIVED: 751231/ REVISED 791231.
 
 
 BRIEF DESCRIPTION:
 
       THIS SECTION CONTAINS THREE PROCEDURES FOR THE SOLUTION
       OF SECOND ORDER SELF-ADJOINT LINEAR TWO POINT
       BOUNDARY VALUE PROBLEMS;
 
       (1) FEM LAG SYM;
 
       THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
 
              - (P(X)*Y')' + R(X)*Y = F(X), A < X < B,
 
         WITH BOUNDARY CONDITIONS
 
             E[1]*Y(A) + E[2]*Y'(A) = E[3],
 
             E[4]*Y(B) + E[5]*Y'(B) = E[6].
 
       (2) FEM LAG;
 
       THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
 
                     - Y'' + R(X)*Y = F(X), A < X < B,
 
         WITH BOUNDARY CONDITIONS
 
             E[1]*Y(A) + E[2]*Y'(A) = E[3],
 
             E[4]*Y(B) + E[5]*Y'(B) = E[6].
 
       (3) FEM LAG SPHER:
 
       THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
 
         WITH SPHERICAL COORDINATES
 
              - (X**NC*Y')'/X**NC + R(X)*Y = F(X),  A < X < B,
 
         WITH BOUNDARY CONDITIONS
 
             E[1]*Y(A) + E[2]*Y'(A) = E[3],
 
             E[4]*Y(B) + E[5]*Y'(B) = E[6].
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 2
 
 
 
 KEY WORDS AND FRASES:
 
       SECOND ORDER DIFFERENTIAL EQUATIONS,
       TWO POINT BOUNDARY VALUE PROBLEMS,
       SELF-ADJOINT BOUNDARY VALUE PROBLEMS,
       RITZ-GALERKIN METHOD,
       SPHERICAL COORDINATES,
       GLOBAL METHODS.
 
 
 LANGUAGE: ALGOL 60.
 
 
 REFERENCES:
 
       [1]  STRANG, G. AND G.J. FIX,
            AN ANALYSIS OF THE FINITE ELEMENT METHOD,
            PRENTICE-HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1973.
 
       [2]  BAKKER, M., EDITOR,
            COLLOQUIUM ON DISCRETIZATION METHODS, CHAPTER 3 (DUTCH),
            MATHEMATISCH CENTRUM, MC-SYLLABUS, TO APPEAR.
 
       [3]  HEMKER, P.W.,
            GALERKIN'S METHOD AND LOBATTO POINTS,
            MATHEMATISCH CENTRUM, REPORT 24/75 (1975).
 
       [4]  BABUSKA, I.,
            NUMERICAL STABILITY IN PROBLEMS OF LINEAR ALGEBRA,
            S.I.A.M.  J. NUM. ANAL., VOL.9, P. 53-77 (1972).
 
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 3
 
 
 
 SUBSECTION: FEM LAG SYM.
 
 
 CALLING SEQUENCE:
 
       THE HEADING OF THE PROCEDURE READS:
 
       "PROCEDURE" FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
       "VALUE" N, ORDER; "INTEGER" N, ORDER;
       "ARRAY" X, Y, E;
       "REAL" "PROCEDURE" P, R, F;
       "CODE" 33300;
 
       THE MEANING OF THE FORMAL PARAMETERS IS:
 
       N:  <ARITHMETIC EXPRESSION>;
           THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
 
       X:  <ARRAY IDENTIFIER>;
           "ARRAY" X[0:N];
           ENTRY: A = X[0] < X[1] < ... < X[N] = B
           IS A PARTITION OF THE INTERVAL [A,B];
 
       Y:  <ARRAY IDENTIFIER>;
           "ARRAY" Y[0:N];
           EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
           SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
 
           (1)      - (P(X)*Y')' + R(X)*Y = F(X), A < X < B,
 
           WITH BOUNDARY CONDITIONS
 
                  E[1]*Y(A) + E[2]*Y'(A) = E[3],
           (2)
                  E[4]*Y(B) + E[5]*Y'(B) = E[6];
 
 
       P:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF P READS:
           "REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X;
           P(X) IS THE COEFFICIENT OF Y' IN (1);
 
       R:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF R READS:
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R(X) IS THE COEFFICIENT OF Y IN (1);
 
       F:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF F READS:
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F(X) IS THE RIGHT HAND SIDE OF (1);
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 4
 
 
 
       ORDER: <ARITHMETIC EXPRESSION>;
           ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
           APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
           THEN  ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
           ORDER CAN BE CHOSEN EQUAL TO 2, 4 OR 6 ONLY;
 
       E:  <ARRAY IDENTIFIER>;
           "ARRAY" E[1:6];
           E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (2);
           E[1] AND E[4] ARE NOT ALLOWWED TO VANISH BOTH.
 
 
 PROCEDURES USED: NONE.
 
 
 REQUIRED CENTRAL MEMORY:
       FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
 
 
 RUNNING TIME:
 
       LET K = ORDER/2; THEN
       (A)   K*N + 1 EVALUATIONS OF P(X), R(X) AND F(X) ARE NEEDED;
       (B)   ABOUT 17*2**(K-1)*N MULTIPLICATIONS/DIVISIONS ARE NEEDED.
 
 
  DATA AND RESULTS:
 
       THE PROCEDURE FEM LAG SYM HAS SOME RESTRICTIONS IN ITS USE:
       (I) P(X) SHOULD BE POSITIVE ON THE CLOSED INTERVAL <X[J-1],X[J]>;
       (II) P(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY SMOOTH
         ON <X[0],X[N]> EXCEPT AT THE GRID POINTS WHERE P(X) SHOULD BE
         AT LEAST CONTINUOUS;
         IN THAT CASE THE ORDER OF ACCURACY (2, 4, OR 6) IS PRESERVED;
       (III) R(X) SHOULD BE NONNEGATIVE ON <X[0],X[N]>;
         IF, HOWEVER, THE PROBLEM HAS PURE DIRICHLET BOUNDARY CONDITIONS
         (I.E. E[2] = E[5] = 0) THIS CONDITION CAN BE WEAKENED TO THE
         REQUIREMENT THAT
 
               R(X) > - P0*(PI/(X[N] - X[0]))**2,
 
         WHERE P0 IS THE MINIMUM OF P(X) ON <X[0],X[N]> AND PI HAS
         THE VALUE 3.14159...; HOWEVER, ONE SHOULD NOTE THAT THE
         PROBLEM MAY BE ILL-CONDITIONED WHEN R(X) IS QUITE NEAR THAT
         LOWER BOUND; FOR OTHER NEGATIVE VALUES OF R(X) THE EXISTENCE
         OF A SOLUTION REMAINS AN OPEN QUESTION;
       (IV) THE USER SHOULD NOT EXPECT GREATER ACCURACY THAN 12 DECIMAL
         DUE TO THE LOSS OF DIGITS DURING THE EVALUATION OF THE MATRIX
         AND THE VECTOR OF THE LINEAR SYSTEM TO BE SOLVED AND DURING ITS
         REDUCTION TO A TRIDIAGONAL SYSTEM; WHEN THE SOLUTION OF THE
         PROBLEM IS NOT TOO WILD, THIS 12-DIGIT ACCURACY CAN ALREADY BE
         OBTAINED WITH A MODERATE MESH SIZE (E.G. < 0.1), PROVIDED THAT
         A SIXTH ORDER METHOD IS USED.
 
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 5
 
 
 
  METHOD AND PERFORMANCE:
 
       PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
       CONTINUOUS PIECEWISE POLYNOMIALS (SEE [1], [2]);
       THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS ON
       THE CLOSED INTERVAL <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS
       THAN OR EQUAL TO K (K = ORDER//2) ON EACH SEGMENT <X[J-1],X[J]>
       (J = 1, ..., N); THIS  PIECEWISE POLYNOMIAL IS ENTIRELY
       DETERMINED BY THE VALUES IT HAS
       AT THE KNOTS X[J] AND ON (K-1) INTERIOR KNOTS ON EACH SEGMENT
       <X[J-1],X[J]>; THESE VALUES ARE OBTAINED BY THE SOLUTION OF AN
       (ORDER + 1)-DIAGONAL LINEAR SYSTEM WITH A SPECIALLY STRUCTURED
       MATRIX (SEE [2]); THE ENTRIES OF THE MATRIX AND THE VECTOR ARE
       INNER PRODUCTS WHICH ARE APPROXIMATED BY PIECEWISE (K+1)-POINT
       LOBATTO QUADRATURE (SEE [3]); THE EVALUATION OF THE MATRIX AND
       THE VECTOR IS DONE SEGMENT BY SEGMENT: ON EACH SEGMENT
       THE CONTRIBUTIONS TO THE ENTRIES OF THE MATRIX AND THE
       VECTOR ARE COMPUTED AND EMBEDDED IN THE GLOBAL MATRIX AND
       VECTOR; SINCE THE FUNCTION VALUES ON THE INTERIOR
       POINTS OF EACH SEGMENT ARE NOT COUPLED WITH THE FUNCTION
       VALUES OUTSIDE THAT SEGMENT, THE RESULTING LINEAR SYSTEM
       CAN BE REDUCED TO A TRIDIAGONAL SYSTEM BY MEANS OF STATIC
       CONDENSATION (SEE [2]); THE FINAL TRIDIAGONAL SYSTEM,
       SINCE IT IS OF FINITE DIFFERENCE TYPE, IS SOLVED BY
       MEANS OF BABUSKA'S METHOD (SEE [4]).
 
 
 EXAMPLE OF USE:
 
       WE SOLVE THE BOUNDARY VALUE PROBLEM
 
 
       -(Y'*EXP(X))'+Y*COS(X)=EXP(X)*(SIN(X)-COS(X))+SIN(2*X)/2,
       0 < X < PI = 3.14159265358979,
       Y(0) = Y(PI) = 0;
 
       FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
 
           E[1] = E[4] = 1; E[2] = E[3] = E[5] = E[6] = 0;
 
       THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
       THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
       I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4,6
       THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 6
 
 
 
       "BEGIN" "INTEGER" N; "FOR" N:= 10, 20 "DO"
       "BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X, Y[0:N], E[1:6];
 
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R:= COS(X);
 
           "REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X;
           P:= EXP(X);
 
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F:= EXP(X)*(SIN(X)-COS(X)) + SIN(2*X)/2;
 
           "PROCEDURE" FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
           "CODE" 33300;
           E[1]:= E[4]:= 1; E[2]:= E[3]:= E[5]:= E[6]:= 0;
           PI:= 3.14159265358979;
           "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
           OUTPUT(61,"("//,6B"("N=")"D")",N);
           "FOR" ORDER:= 2, 4, 6 "DO"
           "BEGIN" "REAL" RHO, D;
               FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
               RHO:= 0;
               "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
               "BEGIN" D:= ABS(Y[I] - SIN(X[I]));
                   "IF" RHO < D "THEN" RHO:= D
               "END";
               OUTPUT(61,"("/,16B"("ORDER=")"DD,4B"("MAX.ERROR= ")",
               D.DD"+ZD")",ORDER,RHO)
           "END"
       "END"
       "END"
 
       RESULTS:
 
       N=10
                 ORDER=2  MAX. ERROR= 1.36" -2
                 ORDER=4  MAX. ERROR= 7.55" -5
                 ORDER=6  MAX. ERROR= 3.48" -8
 
       N=20
                 ORDER=2  MAX. ERROR= 3.41" -3
                 ORDER=4  MAX. ERROR= 4.79" -6
                 ORDER=6  MAX. ERROR= 5.51"-10
 
       ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
       2**(-ORDER)  WHEN THE MESH SIZE IS HALVED.
 
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 7
 
 
 
 SUBSECTION: FEM LAG.
 
 
 CALLING SEQUENCE:
 
       THE HEADING OF THE PROCEDURE READS:
 
       "PROCEDURE" FEM LAG(X, Y, N, R, F, ORDER, E);
       "VALUE" N, ORDER; "INTEGER" N, ORDER;
       "ARRAY" X, Y, E;
       "REAL" "PROCEDURE" R, F;
       "CODE" 33301;
 
       THE MEANING OF THE FORMAL PARAMETERS IS:
 
       N:  <ARITHMETIC EXPRESSION>;
           THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
 
       X:  <ARRAY IDENTIFIER>;
           "ARRAY" X[0:N];
           ENTRY: A = X[0] < X[1] < ... < X[N] = B IS  A
           PARTITION OF THE SEGMENT [A,B];
 
       Y:  <ARRAY IDENTIFIER>;
           "ARRAY" Y[0:N];
           EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
           SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
 
           (3)  - Y''+ R(X)*Y = F(X), A < X < B,
 
           WITH BOUNDARY CONDITIONS
           (4)    E[1]*Y(A) + E[2]*Y'(A) = E[3],
                  E[4]*Y(B) + E[5]*Y'(B) = E[6];
 
       R:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF R READS:
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R(X) IS THE COEFFICIENT OF Y IN (3);
 
       F:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF F READS:
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F(X) IS THE RIGHT HAND SIDE OF (3);
 
       ORDER: <ARITHMETIC <EXPRESSION>;
           ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
           APPROXIMATE SOLUTION OF (3)-(4); LET H = MAX(X[I] - X[I-1]);
           THEN  ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
           ORDER CAN CAN BE CHOSEN EQUAL TO 2, 4 OR 6 ONLY;
 
       E:  <ARRAY IDENTIFIER>;
           "ARRAY" E[1:6];
           E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (4);
           E[1] AND E[4] ARE NOT ALLOWWED TO VANISH BOTH.
 
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 8
 
 
 
 PROCEDURES USED: NONE.
 
 
 REQUIRED CENTRAL MEMORY:
 
       FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
 
 
 RUNNING TIME:
 
       LET K = ORDER/2; THEN
 
       (A)   K*N + 1 EVALUATIONS OF R(X) AND F(X) ARE NEEDED;
 
       (B)   ABOUT 12*2**(K-1)*N MULTIPLICATIONS/DIVISIONS ARE NEEDED.
 
 
 DATA AND RESULTS: SEE PREVIOUS SUBSECTION.
 
 
 METHOD AND PERFORMANCE: SEE PREVIOUS SUBSECTION.
 
 
 EXAMPLE OF USE:
 
       WE SOLVE THE BOUNDARY VALUE PROBLEM
 
 
           - Y'' + Y*EXP(X) = SIN(X)*(1+EXP(X),
           0 < X < PI = 3.14159265358979,
           Y(0) = Y(PI) = 0;
 
       FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
 
           E[1] = E[4] = 1; E[2] = E[3] = E[5] = E[6] = 0;
 
       THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
       THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
       I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4,6
       THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                     PAGE 9
 
 
 
       "BEGIN" "INTEGER" N; "FOR" N:= 10, 20 "DO"
       "BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X, Y[0:N], E[1:6];
 
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R:= EXP(X);
 
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F:= SIN(X)*(1 + EXP(X));
 
           "PROCEDURE" FEM LAG(X, Y, N, R, F, ORDER, E);
           "CODE" 33301;
           E[1]:= E[4]:= 1; E[2]:= E[3]:= E[5]:= E[6]:= 0;
           PI:= 3.14159265358979;
           "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
           OUTPUT(61,"("//,6B"("N=")"D")",N);
           "FOR" ORDER:= 2, 4, 6 "DO"
           "BEGIN" "REAL" RHO, D;
               FEM LAG(X, Y, N, R, F, ORDER, E);
               RHO:= 0;
               "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
               "BEGIN" D:= ABS(Y[I] - SIN(X[I]));
                   "IF" RHO < D "THEN" RHO:= D
               "END";
               OUTPUT(61,"("/,16B"("ORDER=")"DD,4B"("MAX.ERROR= ")",
               D.DD"+ZD")",ORDER,RHO)
           "END"
       "END"
       "END"
 
       RESULTS:
 
       N=10
                 ORDER=2    MAX. ERROR= 1.60" -3
                 ORDER=4    MAX. ERROR= 1.55" -5
                 ORDER=6    MAX. ERROR= 7.28"-10
 
 
       N=20
                 ORDER=2    MAX. ERROR= 4.01" -4
                 ORDER=4    MAX. ERROR= 9.80" -7
                 ORDER=6    MAX. ERROR= 9.38"-12
 
       NOTICE THAT THE MAXIMUM ERROR DECREASES BY ABOUT
       2**(-ORDER) WHEN THE MESH SIZE IS HALVED.
 
 
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 10
 
 
 
 
 SUBSECTION: FEM LAG SPHER.
 
 
 CALLING SEQUENCE:
 
 
       THE HEADING OF THE PROCEDURE READS:
 
       "PROCEDURE" FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
       "VALUE" N, NC, ORDER; "INTEGER" N, NC, ORDER;
       "ARRAY" X, Y, E;
       "REAL" "PROCEDURE" R, F;
       "CODE" 33308;
 
       THE MEANING OF THE FORMAL PARAMETERS IS:
 
       N:  <ARITHMETIC EXPRESSION>;
           THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
 
      NC:   <EXPRESSION>;
            IF NC = 0, CARTESIAN COORDINATES ARE USED;
            IF NC = 1, POLAR COORDINATES ARE USED;
            IF NC = 2, SPHERICAL COORDINATES ARE USED;
 
       X:  <ARRAY IDENTIFIER>;
           "ARRAY" X[0:N];
           ENTRY: A = X[0] < X[1] < ... < X[N] = B IS  A
           PARTITION OF THE INTERVAL [A,B];
 
       Y:  <ARRAY IDENTIFIER>;
           "ARRAY" Y[0:N];
           EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
           SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
 
           (1)      - (X**NC*Y')'/X**NC + R(X)*Y = F(X), A < X < B,
 
           WITH BOUNDARY CONDITIONS
 
                  E[1]*Y(A) + E[2]*Y'(A) = E[3],
           (2)
                  E[4]*Y(B) + E[5]*Y'(B) = E[6];
 
 
       R:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF R READS:
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R(X) IS THE COEFFICIENT OF Y IN (1);
 
       F:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF F READS:
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F(X) IS THE RIGHT HAND SIDE OF (1);
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 11
 
 
       ORDER: <ARITHMETIC EXPRESSION>;
           ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
           APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
           THEN  ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
           ORDER CAN BE CHOSEN EQUAL TO 2 OR 4 ONLY;
 
       E:  <ARRAY IDENTIFIER>;
           "ARRAY" E[1:6];
           E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (2);
           E[1] AND E[4] ARE NOT ALLOWED TO VANISH BOTH.
 
 
 PROCEDURES USED: NONE.
 
 
 REQUIRED CENTRAL MEMORY:
 
       FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
 
 
 RUNNING TIME:
 
       LET K = ORDER/2; THEN
       (A)   K*N EVALUATIONS OF R(X) AND F(X) ARE NEEDED;
       (B)   IF NC > 0 AND ORDER=4, THEN N SQUARE ROOTS ARE EVALUATED;
 
 
  DATA AND RESULTS:
 
       THE PROCEDURE FEM LAG SPHER HAS SOME RESTRICTIONS IN ITS USE:
       R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY SMOOTH
       ON <X[0],X[N]> EXCEPT AT THE GRID POINTS; FURTHERMORE R(X)
       SHOULD BE NONNEGATIVE.
 
 
  METHOD AND PERFORMANCE:
 
       PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
       CONTINUOUS PIECEWISE POLYNOMIALS (SEE [1], [2]);
       THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS ON
       THE CLOSED INTERVAL <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS
       THAN OR EQUAL TO K (K = ORDER//2) ON EACH SEGMENT <X[J-1],X[J]>
       (J = 1, ..., N); THIS  PIECEWISE POLYNOMIAL IS ENTIRELY
       DETERMINED BY THE VALUES IT HAS
       AT THE KNOTS X[J] AND ON (K-1) INTERIOR KNOTS ON EACH SEGMENT
       <X[J-1],X[J]>; THESE VALUES ARE OBTAINED BY THE SOLUTION OF AN
       (ORDER + 1)-DIAGONAL LINEAR SYSTEM WITH A SPECIALLY STRUCTURED
       MATRIX (SEE [2]); THE ENTRIES OF THE MATRIX AND THE VECTOR ARE
       INNER PRODUCTS WHICH ARE APPROXIMATED BY SOME PIECEWISE K-POINT
       GAUSSIAN QUADRATURE (SEE [4]); THE EVALUATION OF THE MATRIX AND
       THE VECTOR IS DONE SEGMENT BY SEGMENT: ON EACH SEGMENT
       THE CONTRIBUTIONS TO THE ENTRIES OF THE MATRIX AND THE VECTOR
       ARE COMPUTED AND EMBEDDED IN THE GLOBAL MATRIX AND VECTOR;
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 12
 
 
 
       SINCE THE FUNCTION VALUES ON THE INTERIOR
       POINTS OF EACH SEGMENT ARE NOT COUPLED WITH THE FUNCTION
       VALUES OUTSIDE THAT SEGMENT, THE RESULTING LINEAR SYSTEM
       CAN BE REDUCED TO A TRIDIAGONAL SYSTEM BY MEANS OF STATIC
       CONDENSATION (SEE [2]); THE FINAL TRIDIAGONAL SYSTEM,
       SINCE IT IS OF FINITE DIFFERENCE TYPE, IS SOLVED BY
       MEANS OF BABUSKA'S METHOD (SEE [3]).
 
 
 EXAMPLE OF USE:
 
       WE SOLVE THE BOUNDARY VALUE PROBLEM
 
 
       -(Y'*X**NC)'/X**NC + Y = 1 - X**4 + (12 + 4*NC)*X**2,
           0 < X < 1; Y'(0) = Y(1) = 0;
 
       FOR THE BOUNDARY CONDITIONS THIS IMPLIES THAT
 
           E[2] = E[4] = 1; E[1] = E[3] = E[5] = E[6] = 0;
 
       THE ANALYTIC SOLUTION IS Y(X) = 1 - X**4; WE APPROXIMATE
       THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I/N,  I = 0, ..., N;
       I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4
       THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
 
       "BEGIN" "INTEGER" N, NC;
           "FOR" N:= 10, 20 "DO" "FOR" NC:= 0, 1, 2 "DO"
       "BEGIN" "INTEGER" I, ORDER; "ARRAY" X, Y[0:N], E[1:6];
 
           "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
           R:= 1;
 
           "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
           F:= (12 + 4*NC)*X**2 + 1 - X**4;
 
           "PROCEDURE" FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
           "CODE" 33308;
 
           E[2]:= E[4]:= 1; E[1]:= E[3]:= E[5]:= E[6]:= 0;
           "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= I/N;
           OUTPUT(61,"("//,6B"("N=")"ZZD,6B"("NC=")"ZD")",N,NC);
           "FOR" ORDER:= 2, 4 "DO"
           "BEGIN" "REAL" RHO, D;
               FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
               RHO:= 0;
               "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
               "BEGIN" D:= ABS(Y[I] - 1 + X[I]**4);
                   "IF" RHO < D "THEN" RHO:= D
               "END";
               OUTPUT(61,"("/,16B"(" ORDER=")"ZD,4B"("MAX.ERROR= ")",
               D.DD"+ZD")",ORDER,RHO)
           "END"
       "END"
       "END"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 13
 
 
 
       RESULTS:
 
       N= 10      NC= 0
                  ORDER= 2    MAX.ERROR= 4.37" -3
                  ORDER= 4    MAX.ERROR= 2.93" -6
 
       N= 10      NC= 1
                  ORDER= 2    MAX.ERROR= 1.42" -2
                  ORDER= 4    MAX.ERROR= 5.49" -5
 
       N= 10      NC= 2
                  ORDER= 2    MAX.ERROR= 2.46" -2
                  ORDER= 4    MAX.ERROR= 1.27" -4
 
       N= 20      NC= 0
                  ORDER= 2    MAX.ERROR= 1.09" -3
                  ORDER= 4    MAX.ERROR= 1.83" -7
 
       N= 20      NC= 1
                  ORDER= 2    MAX.ERROR= 3.53" -3
                  ORDER= 4    MAX.ERROR= 3.91" -6
 
       N= 20      NC= 2
                  ORDER= 2    MAX.ERROR= 6.10" -3
                  ORDER= 4    MAX.ERROR= 9.26" -6
 
       ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
       2**(-ORDER)  WHEN THE MESH SIZE IS HALVED.
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 14
 
 
 
  SOURCE TEXT(S):
0"CODE" 33300;
  "PROCEDURE" FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
  "INTEGER" N, ORDER;
  "REAL" "PROCEDURE" P, R, F;
  "ARRAY" X, Y, E;
  "BEGIN" "INTEGER" L, L1;
    "REAL" XL1, XL, H, A12, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
        P1, P2, P3, P4, R1, R2, R3, R4, F1, F2, F3, F4,
        E1, E2, E3, E4, E5, E6;
    "ARRAY" T, SUB, CHI, GI[0:N-1];
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 1;
    "BEGIN" "REAL"  H2;
        "IF" L=1 "THEN"
        "BEGIN" P2:= P(XL1); R2:= R(XL1); F2:= F(XL1) "END";
        P1:= P2; P2:= P(XL); R1:= R2; R2:= R(XL); F1:= F2; F2:= F(XL);
        H2:= H/2; B1:= H2*F1; B2:= H2*F2; TAU1:= H2*R1; TAU2:= H2*R2;
        A12:= -0.5*(P1 + P2)/H
    "END" ELAN. M.V. EV.;
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 2;
    "BEGIN" "REAL"  X2, H6, H15, B3, TAU3, C12, C32, A13, A22, A23;
        "IF" L=1 "THEN"
        "BEGIN" P3:= P(XL1); R3:= R(XL1); F3:= F(XL1) "END";
        X2:= (XL1 + XL)/2; H6:= H/6; H15:= H/1.5;
        P1:= P3; P2:= P(X2); P3:= P(XL);
        R1:= R3; R2:= R(X2); R3:= R(XL);
        F1:= F3; F2:= F(X2); F3:= F(XL);
        B1:= H6*F1; B2:= H15*F2; B3:= H6*F3;
        TAU1:= H6*R1; TAU2:= H15*R2; TAU3:= H6*R3;
        A12:= -(2*P1 + P3/1.5)/H; A13:= (0.5*(P1 + P3) - P2/1.5)/H;
        A22:= (P1 + P3)/H/0.375 + TAU2; A23:= -(P1/3 + P3)*2/H;
        "COMMENT" STATIC CONDENSATION;
        C12:= - A12/A22; C32:= - A23/A22; A12:= A13 + C32*A12;
        B1:= B1 + C12*B2; B2:= B3 + C32*B2;
        TAU1:= TAU1 + C12*TAU2; TAU2:= TAU3 + C32*TAU2
    "END" ELEMENT MAT VEC EVALUATION 2;
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 3;
    "BEGIN" "REAL"  X2, X3, H12, H24, DET, C12, C13, C42, C43,
        A13, A14, A22, A23, A24, A33, A34, B3, B4, TAU3, TAU4;
        "IF" L=1 "THEN"
        "BEGIN" P4:= P(XL1); R4:= R(XL1); F4:= F(XL1) "END";
        X2:= XL1 + 0.27639320225*H; X3:= XL - X2 + XL1;
        H12:= H/12; H24:= H/2.4;
        P1:= P4; P2:= P(X2); P3:= P(X3); P4:= P(XL);
        R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
        F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
        B1:= H12*F1; B2:= H24*F2; B3:= H24*F3; B4:= H12*F4;
        TAU1:= H12*R1; TAU2:= H24*R2; TAU3:= H24*R3; TAU4:= H12*R4;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 15
                                                                 ;
 
 
        A12:= -(+ 4.04508497187450*P1
                + 0.57581917135425*P3
                + 0.25751416197911*P4)/H;
        A13:= (+ 1.5450849718747*P1
               - 1.5075141619791*P2
               + 0.6741808286458*P4)/H;
        A14:= ((P2 + P3)/2.4 - (P1 + P4)/2)/H;
        A22:= (5.454237476562*P1 + P3/.48 +.79576252343762*P4)/H + TAU2;
        A23:= - (P1 + P4)/(H*0.48);
        A24:= (+ 0.67418082864575*P1
               - 1.50751416197910*P3
               + 1.54508497187470*P4)/H;
        A33:= (.7957625234376*P1 + P2/.48 + 5.454237476562*P4)/H + TAU3;
        A34:= -(+ 0.25751416197911*P1
                + 0.57581917135418*P2
                + 4.0450849718747*P4)/H;
        "COMMENT" STATIC CONDENSATION;
        DET:= A22*A33 - A23*A23;
        C12:= (A13*A23 - A12*A33)/DET;
        C13:= (A12*A23 - A13*A22)/DET;
        C42:= (A23*A34 - A24*A33)/DET;
        C43:= (A24*A23 - A34*A22)/DET;
        TAU1:= TAU1 + C12*TAU2 + C13*TAU3;
        TAU2:= TAU4 + C42*TAU2 + C43*TAU3;
        A12:= A14 + C42*A12 + C43*A13;
        B1:= B1 + C12*B2 + C13*B3;
        B2:= B4 + C42*B2 + C43*B3
    "END" ELEMENT MAT VEC EVALUATION 3;
 
    "PROCEDURE" BOUNDARY CONDITIONS;
    "IF" L=1 "AND" E2 = 0 "THEN"
    "BEGIN"  TAU1:= 1; B1:= E3/E1;B2:= B2 - A12*B1;
        TAU2:= TAU2 - A12; A12:= 0 "END"
    "ELSE" "IF" L=1 "AND" E2 ^= 0 "THEN"
    "BEGIN" "REAL" AUX; AUX:= P1/E2; TAU1:= TAU1 - AUX*E1 ;
        B1:= B1 - E3*AUX
    "END" "ELSE" "IF" L=N "AND" E5 = 0 "THEN"
    "BEGIN" TAU2:= 1; B2:= E6/E4;
        B1:= B1 - A12*B2; TAU1:= TAU1 - A12; A12:= 0
    "END" "ELSE" "IF" L=N "AND" E5 ^= 0 "THEN"
    "BEGIN" "REAL" AUX; AUX:= P2/E5;
        TAU2:= TAU2 + AUX*E4; B2:= B2 + AUX*E6
    "END" B.C.1;
 
    "PROCEDURE" FORWARD BABUSHKA;
    "IF" L=1 "THEN"
    "BEGIN" CHI[0]:= CH:= TL:= TAU1; T[0]:= TL;
         GI[0]:=  G:= YL:=   B1; Y[0]:= YL;
        SUB[0]:= A12; PP:= A12/(CH - A12);
        CH:= TAU2 - CH*PP; G:= B2 - G*PP; TL:= TAU2; YL:= B2
    "END" "ELSE"
    "BEGIN" CHI[L1]:= CH:= CH + TAU1;
         GI[L1]:=  G:=  G + B1;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 16
                                                                 ;
 
 
        SUB[L1]:= A12; PP:= A12/(CH - A12);
        CH:= TAU2 - CH*PP; G:= B2 - G*PP;
        T[L1]:= TL + TAU1; TL:= TAU2;
        Y[L1]:= YL +   B1; YL:=   B2
    "END" FORWARD BABUSHKA 1;
 
    "PROCEDURE" BACKWARD BABUSHKA;
    "BEGIN" PP:= YL; Y[N]:= G/CH;
        G:= PP; CH:= TL; L:= N;
        "FOR" L:= L - 1 "WHILE" L >= 0 "DO"
        "BEGIN" PP:= SUB[L]; PP:= PP/(CH - PP);
            TL:= T[L]; CH:= TL - CH*PP;
            YL:= Y[L]; G:= YL - G*PP;
            Y[L]:=(GI[L] + G - YL)/(CHI[L] + CH - TL)
        "END"
    "END" BACKWARD BABUSHKA;
 
    L:= 0; XL:= X[0];
    E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
    "FOR" L:= L + 1 "WHILE" L <= N "DO"
    "BEGIN" L1:= L - 1; XL1:= XL; XL:= X[L]; H:= XL - XL1;
        "IF" ORDER = 2 "THEN" ELEMENT MAT VEC EVALUATION 1 "ELSE"
        "IF" ORDER = 4 "THEN" ELEMENT MAT VEC EVALUATION 2 "ELSE"
                              ELEMENT MAT VEC EVALUATION 3;
        "IF" L=1 "OR" L=N "THEN" BOUNDARY CONDITIONS;
        FORWARD BABUSHKA
    "END";
    BACKWARD BABUSHKA;
  "END" FEM LAG SYM;
         "EOP"
 
0"CODE" 33301;
  "PROCEDURE" FEM LAG(X, Y, N, R, F, ORDER, E);
  "VALUE" N, ORDER; "INTEGER" N, ORDER;
  "REAL" "PROCEDURE" R, F;
  "ARRAY" X, Y, E;
  "BEGIN" "INTEGER" L, L1;
    "REAL" XL1, XL, H, A12, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
      E1, E2, E3, E4, E5, E6;
    "ARRAY" T, SUB, CHI, GI[0: N-1];
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 1;
    "BEGIN" "OWN" "REAL" F2, R2; "REAL" R1, F1, H2;
        "IF" L=1 "THEN"
        "BEGIN" F2:= F(XL1); R2:= R(XL1) "END";
        A12:= - 1/H; H2:= H/2;
        R1:= R2; R2:= R(XL); F1:= F2; F2:= F(XL);
        B1:= H2*F1; B2:= H2*F2; TAU1:= H2*R1; TAU2:= H2*R2
    "END" ELEMENT MAT VEC EVALUATION 1;
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 17
                                                                 ;
 
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 2;
     "BEGIN" "OWN" "REAL" R3, F3;
        "REAL" R1, R2, F1, F2, X2, H6, H15,
        B3, TAU3, C12, A13, A22, A23;
        "IF" L=1 "THEN"
        "BEGIN" R3:= R(XL1); F3:= F(XL1) "END";
        X2:= (XL1 + XL)/2; H6:= H/6; H15:= H/1.5;
        R1:= R3; R2:= R(X2); R3:= R(XL);
        F1:= F3; F2:= F(X2); F3:= F(XL);
        B1:= H6*F1; B2:= H15*F2; B3:= H6*F3;
        TAU1:= H6*R1; TAU2:= H15*R2; TAU3:= R3*H6;
        A12:= A23:= -8/H/3; A13:= - A12/8; A22:= -2*A12 + TAU2;
        "COMMENT" STATIC CONDENSATION;
        C12:= - A12/A22; A12:= A13 + C12*A12;
        B2:= C12*B2; B1:= B1 + B2; B2:= B3 + B2;
        TAU2:= C12*TAU2; TAU1:= TAU1 + TAU2; TAU2:= TAU3 + TAU2
    "END" ELEMENT MAT VEC EVALUATION2;
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 3;
    "BEGIN" "OWN" "REAL" R4, F4;
        "REAL" R1, R2, R3, F1, F2, F3, X2, X3, H12, H24,
        DET, C12, C13, C42, C43, A13, A14, A22, A23, A24,
        A33, A34, B3, B4, TAU3, TAU4;
        "IF" L=1 "THEN"
        "BEGIN" R4:= R(XL1); F4:= F(XL1) "END";
        X2:= XL1 + 0.27639320225*H; X3:= XL - X2 + XL1;
        R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
        F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
        H12:= H/12; H24:= H/2.4;
        B1:= F1*H12; B2:= F2*H24; B3:= F3*H24; B4:= F4*H12;
        TAU1:= R1*H12; TAU2:= R2*H24; TAU3:= R3*H24; TAU4:= R4*H12;
        A12:= A34:= -4.8784183052078/H; A13:= A24:= 0.7117516385412/H;
        A14:= -0.16666666666667/H; A23:= 25*A14;
        A22:= -2*A23 + TAU2; A33:= -2*A23 + TAU3;
        "COMMENT" STATIC CONDENSATION;
        DET:= A22*A33 - A23*A23;
        C12:= (A13*A23 - A12*A33)/DET;
        C13:= (A12*A23 - A13*A22)/DET;
        C42:= (A23*A34 - A24*A33)/DET;
        C43:= (A24*A23 - A34*A22)/DET;
        TAU1:= TAU1 + C12*TAU2 + C13*TAU3;
        TAU2:= TAU4 + C42*TAU2 + C43*TAU3;
        A12:= A14 + C42*A12 + C43*A13;
        B1:= B1 + C12*B2 + C13*B3;
        B2:= B4 + C42*B2 + C43*B3
    "END" ELEMENT MAT VEC EVALUATION3;
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 18
                                                                 ;
 
 
    "PROCEDURE" BOUNDARY CONDITIONS;
    "IF" L=1 "AND" E2 = 0 "THEN"
    "BEGIN" TAU1:= 1; B1:= E3/E1; B2:= B2 - A12*B1;
        TAU2:= TAU2 - A12; A12:= 0 "END"
    "ELSE" "IF" L=1 "AND" E2 ^= 0 "THEN"
    "BEGIN" TAU1:= TAU1 - E1/E2;
        B1:= B1 - E3/E2
    "END" "ELSE" "IF" L=N "AND" E5 = 0 "THEN"
    "BEGIN" TAU2:= 1; B2:= E6/E4; B1:= B1 - A12*B2;
        TAU1:= TAU1 - A12; A12:= 0
    "END" "ELSE" "IF" L=N "AND" E5 ^= 0 "THEN"
    "BEGIN" TAU2:= TAU2 + E4/E5;
        B2:= B2 + E6/E5
    "END" BOUNDARY CONDITIONS;
 
    "PROCEDURE" FORWARD BABUSHKA;
    "IF" L=1 "THEN"
    "BEGIN" CHI[0]:= CH:= TL:= TAU1; T[0]:= TL;
        GI[0]:= G:= YL:= B1; Y[0]:= YL;
        SUB[0]:= A12; PP:= A12/(CH - A12); CH:= TAU2 - CH*PP;
        G:= B2 - G*PP; TL:= TAU2; YL:= B2
    "END" "ELSE"
    "BEGIN" CHI[L1]:= CH:= CH + TAU1;
        GI[L1]:= G:= G + B1; SUB[L1]:= A12; PP:= A12/(CH - A12);
        CH:= TAU2 - CH*PP; G:= B2 - G*PP;
        T[L1]:= TL + TAU1; TL:= TAU2;
        Y[L1]:= YL +   B1; YL:=   B2
    "END" FORWARD BABUSHKA 1;
 
    "PROCEDURE" BACKWARD BABUSHKA;
    "BEGIN" PP:= YL; Y[N]:= G/CH;
        G:= PP; CH:= TL; L:= N;
        "FOR" L:= L - 1 "WHILE" L >= 0 "DO"
        "BEGIN" PP:= SUB[L]; PP:= PP/(CH - PP);
            TL:= T[L]; CH:= TL - CH*PP;
            YL:= Y[L]; G:= YL - G*PP;
            Y[L]:=((GI[L] + G) - YL)/((CHI[L] + CH) - TL)
        "END"
    "END" BACKWARD BABUSHKA;
 
    L:= 0; XL:= X[0];
    E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
    "FOR" L:= L + 1 "WHILE" L <= N "DO"
    "BEGIN" L1:= L - 1; XL1:= XL; XL:= X[L]; H:= XL - XL1;
        "IF" ORDER = 2 "THEN" ELEMENT MAT VEC EVALUATION 1 "ELSE"
        "IF" ORDER = 4 "THEN" ELEMENT MAT VEC EVALUATION 2 "ELSE"
                              ELEMENT MAT VEC EVALUATION 3;
        "IF" L=1 "OR" L=N "THEN" BOUNDARY CONDITIONS;
        FORWARD BABUSHKA
    "END";
    BACKWARD BABUSHKA;
  "END" FEM LAGR;
         "EOP"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 19
 
 
 
 "CODE" 33308;
  "PROCEDURE" FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
  "VALUE" N, NC, ORDER; "INTEGER" N, NC, ORDER;
  "REAL" "PROCEDURE" R, F;
  "ARRAY" X, Y, E;
  "BEGIN" "INTEGER" L, L1;
    "REAL" XL1, XL, H, A12, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
      TAU3, B3, A13, A22, A23, C32, C12,
        E1, E2, E3, E4, E5, E6;
    "ARRAY" T, SUB, CHI, GI[0:N-1];
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 1;
    "BEGIN" "REAL"  XM, VL, VR,WL, WR, PR, RM, FM, XL2, XLXR, XR2;
      "IF" NC = 0 "THEN" VL:= VR:= 0.5 "ELSE" "IF" NC = 1 "THEN"
      "BEGIN" VL:= (XL1*2 + XL)/6; VR:= (XL1 + XL*2)/6 "END" "ELSE"
      "BEGIN" XL2:= XL1*XL1/12; XLXR:=XL1*XL/6; XR2:=XL*XL/12;
        VL:= 3*XL2 + XLXR + XR2;
        VR:= 3*XR2 + XLXR + XL2
      "END";
 
      WL:= H*VL; WR:=H*VR; PR:= VR/(VL +VR);
      XM:= XL1 + H*PR; FM:= F(XM); RM:=R(XM);
      TAU1:= WL*RM; TAU2:=WR*RM;
      B1:= WL*FM; B2:= WR*FM; A12:= - (VL + VR)/H + H*(1 - PR)*PR*RM
    "END" ELEM. M.V. EV.;
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 2;
    "BEGIN" "REAL" XLM, XRM, VLM, VRM, WLM, WRM, FLM, FRM,
      RLM, RRM, PL1, PL2, PL3, PR1, PR2, PR3, QL1, QL2, QL3,
      RLMPL1, RLMPL2, RLMPL3, RRMPR1, RRMPR2, RRMPR3,
      VLMQL1, VLMQL2, VLMQL3, VRMQR1, VRMQR2, VRMQR3,
      QR1, QR2,QR3;
 
      "IF" NC = 0 "THEN"
      "BEGIN" XLM:=XL1 + H*0.2113248654052; XRM:= XL1 + XL - XLM;
        VLM:= VRM:= 0.5;
        PL1:= PR3:= 0.45534180126148; PL3:= PR1:= -0.12200846792815;
        PL2:= PR2:= 1 - PL1 - PL3;
        QL1:= - 2.15470053837925; QL3:= -0.15470053837925;
        QL2:= - QL1 - QL3; QR1:= - QL3; QR3:= - QL1; QR2:= - QL2;
      "END" "ELSE" "IF" NC = 1 "THEN"
      "BEGIN" "REAL" A, A2, A3, A4, B, B2, B3, B4, P4H,
        P2, P3, P4, AUX1, AUX2;
        A:= XL1; A2:= A*A; A3:= A*A2; A4:= A*A3;
        B:= XL; B2:= B*B; B3:= B*B2; B4:= B*B3;
        P2:= 10*(A2 + 4*A*B + B2); P3:= 6*(A3 + 4*(A2*B + A*B2) + B3);
        P4:= SQRT(6*(A4 + 10*(A*B3 + A3*B) + 28*A2*B2 + B4));
        P4H:= P4*H; XLM:= (P3 - P4H)/P2; XRM:= (P3 + P4H)/P2;
        AUX1:= (A + B)/4; AUX2:= H*(A2 + 7*A*B + B2)/6/P4;
        VLM:= AUX1 - AUX2; VRM:= AUX1 + AUX2;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 20
                                                                       ;
 
 
      "END" "ELSE"
      "BEGIN" "REAL" A, A2, A3, A4, A5, A6, A7, A8,
        B, B2, B3, B4, B5, B6, B7, B8, AB4, A2B3, A3B2, A4B,
        P4, P5, P8, P8H, AUX1, AUX2;
        A:= XL1; A2:= A*A; A3:= A*A2; A4:= A*A3; A5:= A*A4; A6:= A*A5;
          A7:= A*A6; A8:= A*A7;
        B:= XL; B2:= B*B; B3:= B*B2; B4:= B*B3; B5:= B*B4; B6:= B*B5;
          B7:= B*B6; B8:= B*B7;
        AB4:= A*B4; A2B3:= A2*B3; A3B2:= A3*B2; A4B:=A4*B;
        P4:= 15*(A4 + 4*(A3*B + A*B3) + 10*A2*B2 + B4);
        P5:= 10*(A5 + 4*(A4B + AB4) + 10*(A3B2 + A2B3) + B5);
        P8:= SQRT(10*(A8 + 10*(A7*B + A*B7) + 55*(A2*B6 + A6*B2)
            + 164*(A5*B3 +A3*B5) + 290*A4*B4 + B8));
        AUX1:= (A2 +A*B + B2)/6; P8H:= P8*H;
        AUX2:= (H*(A5 + 7*(A4B + AB4) + 28*(A3B2 + A2B3) + B5))/4.8/P8;
        XLM:= (P5 - P8H)/P4; XRM:= (P5 + P8H)/P4;
        VLM:= AUX1 - AUX2; VRM:= AUX1 + AUX2
      "END";
 
      "IF" NC > 0 "THEN"
      "BEGIN" "REAL" AUX, PLM, PRM;
        PLM:= (XLM - XL1)/H; PRM:= (XRM - XL1)/H;
        AUX:= 2*PLM - 1; PL1:= AUX*(PLM - 1); PL3:= AUX*PLM;
          PL2:= 1 - PL1 - PL3;
        AUX:= 2*PRM - 1; PR1:= AUX*(PRM - 1); PR3:= AUX*PRM;
          PR2:= 1 - PR1 - PR3;
        AUX:= 4*PLM; QL1:= AUX - 3; QL3:= AUX - 1; QL2:= - QL1 - QL3;
        AUX:= 4*PRM; QR1:= AUX - 3; QR3:= AUX - 1; QR2:= - QR1 - QR3;
      "END";
 
      WLM:= H*VLM; WRM:= H*VRM; VLM:= VLM/H; VRM:= VRM/H;
      FLM:= F(XLM)*WLM; FRM:= WRM*F(XRM);
      RLM:= R(XLM)*WLM; RRM:= WRM*R(XRM);
      TAU1:= PL1*RLM + PR1*RRM;
      TAU2:= PL2*RLM + PR2*RRM;
      TAU3:= PL3*RLM + PR3*RRM;
      B1:= PL1*FLM + PR1*FRM;
      B2:= PL2*FLM + PR2*FRM;
      B3:= PL3*FLM + PR3*FRM;
      VLMQL1:= QL1*VLM; VRMQR1:= QR1*VRM;
      VLMQL2:= QL2*VLM; VRMQR2:= QR2*VRM;
      VLMQL3:= QL3*VLM; VRMQR3:= QR3*VRM;
      RLMPL1:= RLM*PL1; RRMPR1:= RRM*PR1;
      RLMPL2:= RLM*PL2; RRMPR2:= RRM*PR2;
      RLMPL3:= RLM*PL3; RRMPR3:= RRM*PR3;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 21
                                                                       ;
 
 
      A12:= VLMQL1*QL2 + VRMQR1*QR2 + RLMPL1*PL2 + RRMPR1*PR2;
      A13:= VLMQL1*QL3 + VRMQR1*QR3 + RLMPL1*PL3 + RRMPR1*PR3;
      A22:= VLMQL2*QL2 + VRMQR2*QR2 + RLMPL2*PL2 + RRMPR2*PR2;
      A23:= VLMQL2*QL3 + VRMQR2*QR3 + RLMPL2*PL3 + RRMPR2*PR3;
      "COMMENT" STATIC CONDENSATION;
      C12:= - A12/A22; C32:= - A23/A22; A12:= A13 + C32*A12;
      B1:= B1 + C12*B2; B2:= B3 + C32*B2;
      TAU1:= TAU1 + C12*TAU2; TAU2:= TAU3 + C32*TAU2
    "END" ELEMENT MAT VEC EVALUATION 2;
 
    "PROCEDURE" BOUNDARY CONDITIONS;
    "IF" L=1 "AND" E2 = 0 "THEN"
    "BEGIN"  TAU1:= 1; B1:= E3/E1;B2:= B2 - A12*B1;
        TAU2:= TAU2 - A12; A12:= 0 "END"
    "ELSE" "IF" L=1 "AND" E2 ^= 0 "THEN"
    "BEGIN" "REAL" AUX;
        AUX:= ("IF" NC = 0 "THEN" 1 "ELSE" X[0]**NC)/E2;
        B1:= B1 - E3*AUX; TAU1:= TAU1 - E1*AUX
    "END" "ELSE" "IF" L=N "AND" E5 = 0 "THEN"
    "BEGIN" TAU2:= 1; B2:= E6/E4;
        B1:= B1 - A12*B2; TAU1:= TAU1 - A12; A12:= 0
    "END" "ELSE" "IF" L=N "AND" E5 ^= 0 "THEN"
    "BEGIN" "REAL" AUX;
        AUX:= ("IF" NC = 0 "THEN" 1 "ELSE" X[N]**NC)/E5;
        TAU2:= TAU2 + AUX*E4; B2:= B2 + AUX*E6
    "END" B.C.1;
 
    "PROCEDURE" FORWARD BABUSHKA;
    "IF" L=1 "THEN"
    "BEGIN" CHI[0]:= CH:= TL:= TAU1; T[0]:= TL;
         GI[0]:=  G:= YL:=   B1; Y[0]:= YL;
        SUB[0]:= A12; PP:= A12/(CH - A12);
        CH:= TAU2 - CH*PP; G:= B2 - G*PP; TL:= TAU2; YL:= B2
    "END" "ELSE"
    "BEGIN" CHI[L1]:= CH:= CH + TAU1;
         GI[L1]:=  G:=  G + B1;
        SUB[L1]:= A12; PP:= A12/(CH - A12);
        CH:= TAU2 - CH*PP; G:= B2 - G*PP;
        T[L1]:= TL + TAU1; TL:= TAU2;
        Y[L1]:= YL +   B1; YL:=   B2
    "END" FORWARD BABUSHKA
1SECTION : 5.2.1.2.1.2.1.1    (DECEMBER 1979)                    PAGE 22
                                                                       ;
 
 
    "PROCEDURE" BACKWARD BABUSHKA;
    "BEGIN" PP:= YL; Y[N]:= G/CH;
        G:= PP; CH:= TL; L:= N;
        "FOR" L:= L - 1 "WHILE" L >= 0 "DO"
        "BEGIN" PP:= SUB[L]; PP:= PP/(CH - PP);
            TL:= T[L]; CH:= TL - CH*PP;
            YL:= Y[L]; G:= YL - G*PP;
            Y[L]:=(GI[L] + G - YL)/(CHI[L] + CH - TL)
        "END"
    "END" BACKWARD BABUSHKA;
 
    L:= 0; XL:= X[0];
    E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
    "FOR" L:= L + 1 "WHILE" L <= N "DO"
    "BEGIN" L1:= L - 1; XL1:= XL; XL:= X[L]; H:= XL - XL1;
        "IF" ORDER = 2 "THEN" ELEMENT MAT VEC EVALUATION 1 "ELSE"
                              ELEMENT MAT VEC EVALUATION 2;
        "IF" L=1 "OR" L=N "THEN" BOUNDARY CONDITIONS;
        FORWARD BABUSHKA
    "END";
    BACKWARD BABUSHKA;
  "END" FEM LAG SPHER;
         "EOP"
 
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 1
  
  
  
 AUTHOR: M. BAKKER.
  
  
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
  
  
 RECEIVED: 751231.
  
  
 BRIEF DESCRIPTION:
  
      THIS SECTION CONTAINS A PROCEDURE FOR THE SOLUTION
      OF SECOND ORDER SKEW-ADJOINT LINEAR TWO POINT
      BOUNDARY VALUE PROBLEMS;
  
      FEM LAG SKEW;
  
      THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
  
          - Y'' + Q(X)*Y' + R(X)*Y = F(X), A < X < B,
  
        WITH BOUNDARY CONDITIONS
  
            E[1]*Y(A) + E[2]*Y'(A) = E[3],
  
            E[4]*Y(B) + E[5]*Y'(B) = E[6].
  
  
 KEY WORDS AND FRASES:
  
      SECOND ORDER DIFFERENTIAL EQUATIONS,
      TWO POINT BOUNDARY VALUE PROBLEMS,
      SKEW-ADJOINT BOUNDARY VALUE PROBLEMS,
      GALERKIN'S METHOD,
      GLOBAL METHODS.
  
  
 LANGUAGE: ALGOL 60.
  
  
 REFERENCES:
  
      [1]  STRANG, G. AND G.J. FIX,
           AN ANALYSIS OF THE FINITE ELEMENT METHOD,
           PRENTICE-HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1973.
      [2]  BAKKER, M., EDITOR,
           COLLOQUIUM ON DISCRETIZATION METHODS, CHAPTER 3 (DUTCH),
           MATHEMATISCH CENTRUM, MC-SYLLABUS, TO APPEAR.
      [3]  HEMKER, P.W.,
           GALERKIN'S METHOD AND LOBATTO POINTS,
           MATHEMATISCH CENTRUM, REPORT 24/75 (1975).
      [4]  BABUSKA, I.,
           NUMERICAL STABILITY IN PROBLEMS OF LINEAR ALGEBRA,
           S.I.A.M.  J. NUM. ANAL., VOL.9, P. 53-77 (1972).
  
  
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 2
  
  
  
 SUBSECTION: FEM LAG SKEW.
  
  
 CALLING SEQUENCE:
  
      THE HEADING OF THE PROCEDURE READS:
  
      "PROCEDURE" FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
      "VALUE" N, ORDER; "INTEGER" N, ORDER;
      "ARRAY" X, Y, E;
      "REAL" "PROCEDURE" Q, R, F;
      "CODE" 33302;
  
      THE MEANING OF THE FORMAL PARAMETERS IS:
  
      N:  <ARITHMETIC EXPRESSION>;
          THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
  
      X:  <ARRAY IDENTIFIER>;
          "ARRAY" X[0:N];
          ENTRY: A = X[0] < X[1] < ... < X[N] = B IS  A
          PARTITION OF THE INTERVAL [A,B];
  
      Y:  <ARRAY IDENTIFIER>;
          "ARRAY" Y[0:N];
          EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
          SOLUTION AT X[I] TO THE DIFFERENTIAL EQUATION
  
          (1) - Y'' + Q(X)*Y' + R(X)*Y = F(X), A < X < B,
  
          WITH BOUNDARY CONDITIONS
  
                 E[1]*Y(A) + E[2]*Y'(A) = E[3],
          (2)
                 E[4]*Y(B) + E[5]*Y'(B) = E[6];
  
  
      Q:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF Q READS:
          "REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
          Q(X) IS THE COEFFICIENT OF Y' IN (1);
  
      R:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF R READS:
          "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
          R(X) IS THE COEFFICIENT OF Y IN (1);
  
      F:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF F READS:
          "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
          F(X) IS THE RIGHT HAND SIDE OF (1);
  
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 3
  
  
  
      ORDER: <ARITHMETIC EXPRESSION>;
          ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
          APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
          THEN  ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
          ORDER CAN CAN BE CHOSEN EQUAL TO 2, 4 OR 6 ONLY;
  
      E:  <ARRAY IDENTIFIER>;
          "ARRAY" E[1:6];
          E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (2);
          E[1] AND E[4] ARE NOT ALLOWWED TO VANISH BOTH.
  
  
 PROCEDURES USED: NONE.
  
  
 REQUIRED CENTRAL MEMORY:
  
      FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
  
  
 RUNNING TIME:
  
      LET K = ORDER/2; THEN
  
      (A)   K*N + 1 EVALUATIONS OF Q(X), R(X) AND F(X) ARE NEEDED;
  
      (B)   ABOUT 17*2**(K-1)*N MULTIPLICATIONS/DIVISIONS ARE NEEDED.
  
  
 DATA AND RESULTS:
      THE PROCEDURE FEM LAG SKEW HAS SOME RESTRICTIONS IN ITS USE:
      (I)  Q(X) IS NOT ALLOWED TO HAVE VERY LARGE VALUES IN SOME SENSE:
        THE PRODUCT Q(X)*(X[J] - X[J-1]) SHOULD NOT BE TOO LARGE
        ON THE CLOSED INTERVAL <X[J-1],X[J]>, OTHERWISE
        THE BOUNDARY VALUE PROBLEM MAY DEGENERATE TO A SINGULAR
        PERTURBATION OR BOUNDARY LAYER PROBLEM, FOR WHICH EITHER
        SPECIAL METHODS OR A SUITABLY CHOSEN GRID ARE NEEDED;
      (II)  Q(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY
        DIFFERENTIABLE ON THE DOMAIN OF THE BOUNDARY VALUE PROBLEM;
        THEY ARE, HOWEVER, THE DERIVATIVES ARE ALLOWED TO HAVE
        DISCONTINUITIES AT THE GRID POINTS, IN WHICH CASE THE ORDER OF
        ACCURACY (2, 4 OR 6) IS PRESERVED;
      (III)  IF Q(X) AND R(X) SATISFY THE INEQUALITY R(X) >= Q'(X)/2,
        THE EXISTENCE OF A UNIQUE SOLUTION IS GUARANTEED, OTHERWISE
        THIS REMAINS AN OPEN QUESTION;
      (IV) THE USER SHOULD NOT EXPECT GREATER ACCURACY THAN 12 DECIMALS
        DUE TO THE LOSS OF DIGITS DURING THE EVALUATION OF THE MATRIX
        AND THE VECTOR OF THE LINEAR SYSTEM TO BE SOLVED AND DURING ITS
        REDUCTION TO A TRIDIAGONAL SYSTEM; WHEN THE SOLUTION OF THE
        PROBLEM IS NOT TOO WILD, THIS 12-DIGITS ACCURACY CAN BE OBTAINED
        WITH A MODERATE MESH SIZE (E.G. < 0.1) ALREADY, PROVIDED
        A SIXTH ORDER METHOD IS USED.
  
  
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 4
  
  
  
 METHOD AND PERFORMANCE:
  
      PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
      CONTINUOUS PIECEWISE POLYNOMIAL FUNCTIONS (SEE [1], [2]);
      THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS
      ON THE INTERVAL <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS THAN
      OR EQUAL TO K (K = ORDER//2) ON EACH SEGMENT <X[J-1],X[J]>
      (J = 1, ..., N); THIS
      PIECEWISE POLYNOMIAL IS ENTIRELY DETERMINED BY THE VALUES IT HAS
      AT THE KNOTS X[J] AND ON (K-1) INTERIOR KNOTS ON EACH SEGMENT
      <X[J-1],X[J]>; THESE VALUES ARE OBTAINED BY THE SOLUTION OF AN
      (ORDER + 1)-DIAGONAL LINEAR SYSTEM WITH A SPECIALLY STRUCTURED
      MATRIX (SEE [2]); THE ENTRIES OF THE MATRIX AND THE VECTOR ARE
      INNER PRODUCTS WHICH ARE APPROXIMATED BY PICEWISE (K+1)-POINT
      LOBATTO QUADRATURE (SEE [3]); THE EVALUATION OF THE MATRIX AND
      THE VECTOR IS DONE SEGMENT BY SEGMENT: ON EACH SEGMENT
      THE CONTRIBUTIONS TO THE ENTRIES OF THE MATRIX AND THE
      VECTOR ARE COMPUTED AND EMBEDDED IN THE GLOBAL MATRIX AND
      VECTOR; SINCE THE FUNCTION VALUES ON THE INTERIOR
      POINTS OF EACH SEGMENT ARE NOT COUPLED WITH THE FUNCTION
      VALUES OUTSIDE THAT SEGMENT, THE RESULTING LINEAR SYSTEM
      CAN BE REDUCED TO A TRIDIAGONAL SYSTEM BY MEANS OF STATIC
      CONDENSATION (SEE [2]); SINCE THE FINAL TRIDIAGONAL SYSTEM
      IS OF FINITE DIFFERENCE TYPE, IT IS SOLVED BY MEANS
      OF BABUSKA'S METHOD (SEE [4]).
  
  
 EXAMPLE OF USE:
  
      WE SOLVE THE BOUNDARY VALUE PROBLEM
  
  
      - Y'' + Y'*COS(X) + Y*EXP(X) = SIN(X)*(1 + EXP(X)) + COS(X)**2,
          0 < X < PI = 3.14159265358979, Y(0) = Y(PI) = 0;
  
      FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
  
          E[1] = E[4] = 1; E[2] = E[3] = E[5] = E[6] = 0;
  
      THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
      THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
      I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4,6
      THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
  
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 5
  
  
  
      "BEGIN" "INTEGER" N; "FOR" N:= 10, 20 "DO"
      "BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X, Y[0:N], E[1:6];
  
          "REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
          Q:= COS(X);
  
          "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
          R:= EXP(X);
  
          "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
          F:= SIN(X)*(1 + EXP(X)) + COS(X)**2;
  
          "PROCEDURE" FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
          "CODE" 33302;
          E[1]:= E[4]:= 1; E[2]:= E[3]:= E[5]:= E[6]:= 0;
          PI:= 3.14159265358979;
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
          OUTPUT(61,"("//,6B"("N=")"D")",N);
          "FOR" ORDER:= 2, 4, 6 "DO"
          "BEGIN" "REAL" RHO, D;
              FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
              RHO:= 0;
              "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
              "BEGIN" D:= ABS(Y[I] - SIN(X[I]));
                  "IF" RHO < D "THEN" RHO:= D
              "END";
              OUTPUT(61,"("/,16B"("ORDER=")"DD,4B"("MAX.ERROR= ")",
              D.DD"+ZD")",ORDER,RHO)
          "END"
      "END"
      "END"
  
      RESULTS:
  
      N=10
                ORDER=2  MAX. ERROR= 2.95" -3
                ORDER=4  MAX. ERROR= 2.56" -5
                ORDER=6  MAX. ERROR= 4.26" -8
  
      N=20
                ORDER=2  MAX. ERROR= 7.55" -4
                ORDER=4  MAX. ERROR= 1.68" -6
                ORDER=6  MAX. ERROR= 6.76"-10
  
      NOTICE THAT THE MAXIMUM ERROR DECREASES BY ABOUT
      2**(-ORDER)  WHEN THE MESH SIZE IS HALVED.
  
  
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 6
  
  
  
 SOURCE TEXT(S):
  
0"CODE" 33302;
  "PROCEDURE" FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
  "INTEGER" N, ORDER;
  "REAL" "PROCEDURE" Q, R, F;
  "ARRAY" X, Y, E;
  "BEGIN" "INTEGER" L, L1;
    "REAL" XL1, XL, H, A12, A21, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
      E1, E2, E3, E4, E5, E6;
    "ARRAY" T, SUPER, SUB, CHI, GI[0:N-1];
  
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 1;
    "BEGIN" "OWN" "REAL" Q2, R2, F2;
       "REAL" Q1, R1, F1, H2, S12;
       "IF" L=1 "THEN"
       "BEGIN" Q2:= Q(XL1); R2:= R(XL1); F2:= F(XL1) "END";
       H2:= H/2; S12:= - 1/H;
       Q1:= Q2; Q2:= Q(XL);
       R1:= R2; R2:= R(XL);
       F1:= F2; F2:= F(XL);
       B1:= H2*F1; B2:= H2*F2;
       TAU1:= H2*R1; TAU2:= H2*R2;
       A12:= S12 + Q1/2; A21:= S12 - Q2/2
    "END" ELEMENT MAT VEC EV.;
  
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 2;
    "BEGIN" "OWN" "REAL" Q3, R3, F3;
       "REAL" Q1, Q2, R1, R2, F1, F2, S12, S13, S22, X2, H6, H15,
         C12, C32, A13, A31, A22, A23, A32, B3, TAU3;
       "IF" L=1 "THEN"
       "BEGIN" Q3:= Q(XL1); R3:= R(XL1); F3:= F(XL1) "END";
  
       X2:= (XL1 + XL)/2; H6:= H/6; H15:= H/1.5;
       Q1:= Q3; Q2:= Q(X2); Q3:= Q(XL);
       R1:= R3; R2:= R(X2); R3:= R(XL);
       F1:= F3; F2:= F(X2); F3:= F(XL);
       B1:= H6*F1; B2:= H15*F2; B3:= H6*F3;
       TAU1:= H6*R1; TAU2:= H15*R2; TAU3:= H6*R3;
       S12:= - 1/H/0.375; S13:= - S12/8; S22:= - 2*S12;
       A12:= S12 + Q1/1.5; A13:= S13 - Q1/6;
       A21:= S12 - Q2/1.5; A23:= S12 + Q2/1.5; A22:= S22 + TAU2;
       A31:= S13 + Q3/6; A32:= S12 - Q3/1.5;
       "COMMENT" STATIC CONDENSATION;
       C12:= - A12/A22; C32:= - A32/A22;
       A12:= A13 + C12*A23; A21:= A31 + C32*A21;
       B1:= B1 + C12*B2; B2:= B3 + C32*B2;
       TAU1:= TAU1 + C12*TAU2; TAU2:= TAU3 + C32*TAU2
    "END" ELEMENT MAT VEC EVALUATION 2;
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 7
                                                                  ;
  
  
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 3;
    "BEGIN" "OWN" "REAL" Q4, R4, F4;
       "REAL" Q1, Q2, Q3, R1, R2, R3, F1, F2, F3,
        S12, S13, S14, S22, S23, X2, X3, H12, H24,
        DET, C12, C13, C42, C43, A13, A14, A22, A23,
        A24, A31, A32, A33, A34, A41, A42, A43,
        B3, B4, TAU3, TAU4;
  
       "IF" L=1 "THEN"
       "BEGIN" Q4:= Q(XL1); R4:= R(XL1); F4:= F(XL1) "END";
       X2:= XL1 + 0.27639320225*H; X3:= XL - X2 + XL1;
       H12:= H/12; H24:= H/2.4;
       Q1:= Q4; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(XL);
       R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
       F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
       S12:= -4.8784183052080/H; S13:= 0.7117516385414/H;
       S14:= -.16666666666667/H; S23:= 25*S14; S22:= -2*S23;
       B1:= H12*F1; B2:= H24*F2; B3:= H24*F3; B4:= H12*F4;
       TAU1:= H12*R1; TAU2:= H24*R2; TAU3:= H24*R3; TAU4:= H12*R4;
       A12:= S12 + 0.67418082864578*Q1;
       A13:= S13 - 0.25751416197912*Q1;
       A14:= S14 + Q1/12;
       A21:= S12 - 0.67418082864578*Q2;
       A22:= S22 + TAU2;
       A23:= S23 + 0.93169499062490*Q2;
       A24:= S13 - 0.25751416197912*Q2;
       A31:= S13 + 0.25751416197912*Q3;
       A32:= S23 - 0.93169499062490*Q3;
       A33:= S22 + TAU3;
       A34:= S12 + 0.67418082864578*Q3;
       A41:= S14 - Q4/12;
       A42:= S13 + 0.25751416197912*Q4;
       A43:= S12 - 0.67418082864578*Q4;
       "COMMENT" STATIC CONDENSATION;
       DET:= A22*A33 - A23*A32;
       C12:= (A13*A32 - A12*A33)/DET;
       C13:= (A12*A23 - A13*A22)/DET;
       C42:= (A32*A43 - A42*A33)/DET;
       C43:= (A42*A23 - A43*A22)/DET;
       TAU1:= TAU1 + C12*TAU2 + C13*TAU3 ;
       TAU2:= TAU4 + C42*TAU2 + C43*TAU3;
       A12:= A14 + C12*A24 + C13*A34;
       A21:= A41 + C42*A21 + C43*A31;
       B1:= B1 + C12*B2 + C13*B3;
       B2:= B4 + C42*B2 + C43*B3
    "END" ELEMENT MAT VEC EVALUATION 3;
1SECTION : 5.2.1.2.1.2.1.2    (JANUARY 1976)                      PAGE 8
                                                                  ;
  
  
    "PROCEDURE" BOUNDARY CONDITIONS;
    "IF" L=1 "AND" E2 = 0 "THEN"
    "BEGIN" TAU1:= 1; B1:= E3/E1; A12:= 0 "END"
    "ELSE" "IF" L=1 "AND" E2 ^= 0 "THEN"
    "BEGIN" TAU1:= TAU1 - E1/E2; B1:= B1 - E3/E2
    "END" "ELSE" "IF" L=N "AND" E5 = 0 "THEN"
    "BEGIN" TAU2:= 1; A21:= 0; B2:= E6/E4;
    "END" "ELSE" "IF" L=N "AND" E5 ^= 0 "THEN"
    "BEGIN" TAU2:= TAU2 + E4/E5; B2:= B2 + E6/E5
    "END" B.C.1;
  
    "PROCEDURE" FORWARD BABUSKA;
    "IF" L=1 "THEN"
    "BEGIN" CHI[0]:= CH:= TL:= TAU1; T[0]:= TL;
         GI[0]:= G:= YL:= B1; Y[0]:= YL;
         SUB[0]:= A21; SUPER[0]:= A12;
         PP:= A21/(CH - A12); CH:= TAU2 - CH*PP;
         G:= B2 - G*PP; TL:= TAU2; YL:= B2
    "END" "ELSE"
    "BEGIN" CHI[L1]:= CH:= CH + TAU1;
         GI[L1]:= G:= G + B1;
         SUB[L1]:= A21; SUPER[L1]:= A12;
         PP:= A21/(CH - A12); CH:= TAU2 - CH*PP;
         G:= B2 - G*PP; T[L1]:= TL + TAU1; TL:= TAU2;
         Y[L1]:= YL +  B1; YL:= B2
    "END" FORWARD BABUSKA;
  
    "PROCEDURE" BACKWARD BABUSKA;
    "BEGIN"PP:= YL; Y[N]:= G/CH;
         G:= PP; CH:= TL; L:= N;
         "FOR" L:= L - 1 "WHILE" L >= 0 "DO"
         "BEGIN" PP:= SUPER[L]/(CH - SUB[L]);
           TL:= T[L]; CH:= TL - CH*PP;
           YL:= Y[L]; G:= YL - G*PP;
           Y[L]:=(GI[L] + G - YL)/(CHI[L] + CH - TL)  ;
         "END"
    "END" BACKWARD BABUSKA;
  
    L:= 0; XL:= X[0];
    E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
    "COMMENT" ELEMENTWISE ASSEMBLAGE OF MATRIX AND VECTOR
      COMBINED WITH FORWARD BABUSKA SUBSTITUTION;
    "FOR" L:= L + 1 "WHILE" L <= N "DO"
    "BEGIN" XL1:= XL; L1:= L - 1; XL:= X[L]; H:= XL - XL1;
       "IF" ORDER = 2 "THEN" ELEMENT MAT VEC EVALUATION 1 "ELSE"
       "IF" ORDER = 4 "THEN" ELEMENT MAT VEC EVALUATION 2 "ELSE"
       ELEMENT MAT VEC EVALUATION 3;
       "IF" L=1 "OR" L=N "THEN" BOUNDARY CONDITIONS;
       FORWARD BABUSKA
    "END";
    BACKWARD BABUSKA;
  "END" FEM LAGR;
        "EOP"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 1
  
  
  
 AUTHOR: M. BAKKER.
  
  
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
  
  
 RECEIVED: 751231.
  
  
 BRIEF DESCRIPTION:
  
      THIS SECTION CONTAINS A PROCEDURE FOR THE SOLUTION
      OF FOURTH ORDER SELF-ADJOINT LINEAR TWO POINT
      BOUNDARY VALUE PROBLEMS;
  
      FEM HERM SYM;
  
      THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
  
         (P(X)*Y'')'' - (Q(X)*Y')' + R(X)*Y = F(X),  A < X < B,
  
        WITH BOUNDARY CONDITIONS
  
            Y(A) = E[1], Y'(A) = E[2],
  
            Y(B) = E[3], Y'(B) = E[4].
  
  
 KEY WORDS AND FRASES:
  
      FOURTH ORDER DIFFERENTIAL EQUATIONS,
      TWO POINT BOUNDARY VALUE PROBLEMS,
      SELF-ADJOINT BOUNDARY VALUE PROBLEMS,
      GALERKIN'S METHOD,
      DIRICHLET BOUNDARY CONDITIONS,
      GLOBAL METHODS.
  
  
 LANGUAGE: ALGOL 60.
  
  
 REFERENCES:
  
      [1]  STRANG, G. AND G.J. FIX,
           AN ANALYSIS OF THE FINITE ELEMENT METHOD,
           PRENTICE-HALL, ENGLE WOOD CLIFFS, NEW JERSEY, 1973.
  
      [2]  BAKKER, M., EDITOR,
           COLLOQUIUM ON DISCRETIZATION METHODS, CHAPTER 3 (DUTCH),
           MATHEMATISCH CENTRUM, MC-SYLLABUS, TO APPEAR.
  
      [3]  HEMKER, P.W.,
           GALERKIN'S METHOD AND LOBATTO POINTS,
           MATHEMATISCH CENTRUM, REPORT 24/75 (1975).
  
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 2
  
  
  
 CALLING SEQUENCE:
  
      THE HEADING OF THE PROCEDURE READS:
  
      "PROCEDURE" FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
      "VALUE" N, ORDER; "INTEGER" N, ORDER;
      "ARRAY" X, Y, E;
      "REAL" "PROCEDURE" P, Q, R, F;
      "CODE" 33303;
  
      THE MEANING OF THE FORMAL PARAMETERS IS:
  
      N:  <ARITHMETIC EXPRESSION>;
          THE UPPER BOUND OF THE ARRAY X; N > 1;
  
      X:  <ARRAY IDENTIFIER>;
          "ARRAY" X[0:N];
          ENTRY: A = X[0] < X[1] < ... < X[N] = B IS  A
          PARTITION OF THE INTERVAL [A,B];
  
      Y:  <ARRAY IDENTIFIER>;
          "ARRAY" Y[1:2*N-2];
          EXIT: Y[2*I-1] IS AN APPROXIMATION TO Y(X[I]),
          Y[2*I] IS AN APPROXIMATION TO Y'(X[I]),
          WHERE Y(X) IS THE SOLUTION OF THE DIFFERENTIAL EQUATION
  
          (1) (P(X)*Y'')'' - (Q(X)*Y')' + R(X)*Y = F(X)  , A< X < B,
  
          WITH BOUNDARY CONDITIONS
  
                 Y(A) = E[1], Y'(A) = E[2],
          (2)
                 Y(B) = E[3], Y'(B) = E[4];
  
  
      P:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF P READS:
          "REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X;
          P(X) IS THE COEFFICIENT OF Y'' IN (1);
          P(X) SHOULD BE STRICTLY POSITIVE;
  
      Q:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF Q READS:
          "REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
          Q(X) IS THE COEFFICIENT OF Y' IN (1);
          Q(X) SHOULD BE NONNEGATIVE;
  
      R:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF R READS:
          "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
          R(X) IS THE COEFFICIENT OF Y IN (1);
          R(X) SHOULD BE NONNEGATIVE;
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 3
  
  
  
      F:  <PROCEDURE IDENTIFIER>;
          THE HEADING OF F READS:
          "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
          F(X) IS THE RIGHT HAND SIDE OF (1);
  
      ORDER: <ARITHMETIC EXPRESSION>;
          ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
          APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
          THEN
              ABS(Y[2*I-1]-Y(X[I])) <= C1 * H**ORDER,
              ABS(Y[2*I]-Y'(X[I])   <= C2 * H**ORDER, I = 1,...,N-1;
          ORDER CAN ONLY BE CHOSEN EQUAL TO 4, 6, 8;
  
      E:  <ARRAY IDENTIFIER>;
          "ARRAY" E[1:4];
          E[1], ... , E[4] DESCRIBE THE BOUNDARY CONDITIONS (2).
  
  
 PROCEDURES USED: CHLDECSOLBND = CP 34333
  
  
 REQUIRED CENTRAL MEMORY:
  
      ONE AUXILIARY ARRAY OF 8*(N-1) REALS IS USED.
  
  
 RUNNING TIME:
  
      LET K = ORDER/2; THEN
  
      (A)   K*N + 1 EVALUATIONS OF P(X), Q(X), R(X) AND F(X) ARE NEEDED;
  
      (B)   ABOUT (ORDER-3)*50*N MULTIPLICATIONS/DIVISIONS ARE NEEDED;
  
      (C)   ONE CALL OF CHLDECSOLBND IS DONE.
  
  
 DATA AND RESULTS:
  
      THE PROCEDURE FEM HERM SYM HAS SOME RESTRICTIONS:
      (I)  P(X) SHOULD BE POSITIVE ON THE CLOSED
        INTERVAL <X[0],X[N]> AND Q(X) AND R(X) SHOULD BE
        NONNEGATIVE THERE;
      (II)  P(X), Q(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY
        SMOOTH ON THE  INTERVAL <X[0],X[N]> EXCEPT AT THE KNOTS,
        WHERE DISCONTINUITIES OF THE DERIVATIVES ARE ALLOWED;
        IN THAT CASE THE ORDER OF ACCURACY IS PRESERVED;
      (III)  THE USER SHOULD NOT EXPECT HIGHER ACCURACY THAN 12
        DECIMALS DUE TO THE LOSS OF DIGITS DURING THE EVALUATION OF THE
        MATRIX AND VECTOR AND DURING THE REDUCTION TO A PENTADIAGONAL
        SYSTEM; THIS ACCURACY CAN BE REACHED VERY EASILY WHEN AN EIGTH
        ORDER METHOD IS USED
  
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 4
  
  
  
 METHOD AND PERFORMANCE:
  
      PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
      CONTINUOUSLY DIFFERENTIABLE PIECEWISE POLYNOMIAL FUNCTIONS
      (SEE [1], [2]) : THE SOLUTION IS APPROXIMATED BY A FUNCTION
      WHICH IS CONTINUOUSLY DIFFERENTIABLE ON THE CLOSED INTERVAL
      <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS THAN OR EQUAL TO
      K (K = 1 + ORDER//2) ON EACH CLOSED SEGMENT <X[J-1],X[J]>
      (J = 1, ..., N);
      THIS FUNCTION IS ENTIRELY DETERMINED BY THE VALUES OF
      THE ZEROETH AND FIRST DERIVATIVE AT THE KNOTS X[J] AND BY
      THE VALUES IT HAS AT (K-3) INTERIOR KNOTS ON EACH CLOSED
      SEGMENT <X[J-1],X[J]>; THE VALUES OF THE FUNCTION AND ITS
      DERIVATIVE AT THE KNOTS ARE OBTAINED BY THE SOLUTION OF AN
      (ORDER + 1)-DIAGONAL LINEAR SYSTEM OF (K-1)*N - 2 UNKNOWNS;
      THE ENTRIES OF THE MATRIX AND THE VECTOR ARE INNER PRODUCTS
      WHICH ARE APPROXIMATED BY PIECEWISE K-POINT LOBATTO
      QUADRATURE (SEE [3]); THE EVALUATION OF THE MATRIX AND
      VECTOR IS PERFORMED SEGMENT BY SEGMENT;
      IF K > 3 THE RESULTING LINEAR SYSTEM CAN BE REDUCED
      TO A PENTADIAGONAL SYSTEM BY MEANS OF STATIC
      CONDENSATION; THIS IS POSSIBLE BECAUSE THE FUNCTION
      VALUES AT THE INTERIOR KNOTS ON EACH SEGMENT <X[J-1],X[J]>
      DO NOT DEPEND ON FUNCTION VALUES OUTSIDE THAT SEGMENT;
      THE FINAL PENTADIAGONAL SYSTEM, SINCE THE MATRIX IS POSITIVE
      DEFINITE AND SYMMETRIC, IS SOLVED BY MEANS OF CHOLESKY'S
      DECOMPOSITION METHOD (SEE SECTION  3.1.2.1.1.2.1.3).
  
  
 EXAMPLE OF USE:
  
      WE SOLVE THE BOUNDARY VALUE PROBLEM
  
      Y'''' - (Y'*COS(X))' + Y*EXP(X) = SIN(X)*(1 + EXP(X) + COS(X)*2),
          0 < X < PI;
          Y(0) = Y(PI) = 0; Y'(0) = 1; Y'(PI) = -1;
                            PI = 3.14159265358979;
  
      THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
      THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
      I = 0, ..., N; WE CHOOSE N = 5, 10 AND WE COMPUTE
      THE MAXIMUM DEVIATIONS FROM Y(X[I]) AND Y'(X[I])
      FOR ORDER = 4, 6, 8;
      THE PROGRAM READS AS FOLLOWS:
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 5
  
  
  
      "BEGIN" "INTEGER" N; "FOR" N:= 5, 10 "DO"
      "BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X[0:N],
                                           Y[1:2*N-2], E[1:4];
  
          "REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X; P:= 1;
  
          "REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
          Q:= COS(X);
  
          "REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
          R:= EXP(X);
  
          "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
          F:= SIN(X)*(1 + EXP(X)+ 2*COS(X));
  
          "PROCEDURE" FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
          "CODE" 33303;
          E[1]:= E[3]:= 0; E[2]:= 1; E[4]:= - 1;
          PI:= 3.14159265358979;
          "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
          OUTPUT(61,"("//,6B"("N=")"ZD")",N);
          "FOR" ORDER:= 4, 6, 8 "DO"
          "BEGIN" "REAL" RHO1, RHO2, D1, D2;
              FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
              RHO1:= RHO2:= 0;
              "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
              "BEGIN" D1:= ABS(Y[2*I-1] - SIN(X[I]));
                  "IF" RHO1 < D1 "THEN" RHO1:= D1;
                  D2:= ABS(Y[2*I] - COS(X[I]));
                  "IF" RHO2 < D2 "THEN" RHO2:= D2
              "END";
              OUTPUT(61,"("/,16B"("ORDER=")"D,/,
              24B"("MAX ABS(Y[2*I-1]-Y(X[I]))= ")",D.3D"+ZD,
              /,24B"("MAX ABS(Y[2*I]-Y'(X[I]))= ")",D.3D"+ZD")",
                   ORDER,RHO1,RHO2)
          "END"
      "END"
      "END"
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 6
  
  
  
      RESULTS:
  
      N= 5
           ORDER=4
                MAX ABS(Y[2*I-1]-Y(X[I])) = 4.822" -4
                MAX ABS(Y[2*I]-Y'(X[I]))  = 4.548" -4
           ORDER=6
                MAX ABS(Y[2*I-1]-Y(X[I])) = 5.651" -6
                MAX ABS(Y[2*I]-Y'(X[I]))  = 2.035" -6
           ORDER=8
                MAX ABS(Y[2*I-1]-Y(X[I])) = 2.264" -8
                MAX ABS(Y[2*I]-Y'(X[I]))  = 1.600" -8
  
      N=10
           ORDER=4
                MAX ABS(Y[2*I-1]-Y(X[I])) = 2.657" -5
                MAX ABS(Y[2*I]-Y'(X[I]))  = 2.870" -5
           ORDER=6
                MAX ABS(Y[2*I-1]-Y(X[I])) = 8.398" -8
                MAX ABS(Y[2*I]-Y'(X[I]))  = 3.572" -8
           ORDER=8
                MAX ABS(Y[2*I-1]-Y(X[I])) = 7.981"-11
                MAX ABS(Y[2*I]-Y'(X[I]))  = 6.796"-11
  
      NOTICE THAT THE MAXIMUM ERROR IS DIVIDED BY
      2**ORDER, WHEN THE MESH SIZE IS HALVED.
  
  
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 7
  
  
  
 SOURCE TEXT(S):
  
0"CODE" 33303;
  "PROCEDURE"  FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
  "VALUE" N, ORDER; "INTEGER" N, ORDER;
  "ARRAY" X, Y, E;
  "REAL" "PROCEDURE"  P, Q, R, F;
  "BEGIN" "INTEGER" L, N2, V, W;
    "ARRAY" A[1:8*(N - 1)], EM[2:3];
    "REAL" A11, A12, A13, A14, A22, A23, A24, A33, A34, A44,
      YA, YB, ZA, ZB,
      B1, B2, B3, B4, D1, D2, E1, R1, R2, XL1, XL;
  
    "PROCEDURE" CHLDECSOLBND(A, N, W, AUX, B); "CODE"34333;
  
    "PROCEDURE" ELEMENTMATVECEVALUATION;
    "IF"ORDER=4"THEN"
    "BEGIN" "REAL" X2, H, H2, H3, P1, P2,
        Q1, Q2, R1, R2, F1, F2,
        B11, B12, B13, B14, B22, B23, B24, B33, B34, B44,
        S11, S12, S13, S14, S22, S23, S24, S33, S34, S44,
        M11, M12, M13, M14, M22, M23, M24, M33, M34, M44;
        "OWN" "REAL"P3, Q3, R3, F3;
  
        H:= XL - XL1; H2:= H*H; H3:= H*H2;
        X2:= (XL1 + XL)/2;
        "IF"L=1"THEN"
        "BEGIN"P3:= P(XL1); Q3:= Q(XL1); R3:= R(XL1); F3:= F(XL1)
        "END";
  
        "COMMENT" ELEMENT BENDING MATRIX;
        P1:= P3; P2:= P(X2); P3:= P(XL);
        B11:= 6*(P1 + P3); B12:= 4*P1 + 2*P3;
        B13:= - B11; B14:= B11 - B12;
        B22:= (4*P1 + P2 + P3)/1.5; B23:= - B12; B24:= B12 - B22;
        B33:= B11; B34:= - B14; B44:= B14 - B24;
  
        "COMMENT" ELEMENT STIFFNESS MATRIX;
        Q1:= Q3; Q2:= Q(X2); Q3:= Q(XL);
        S11:= 1.5*Q2; S12:= Q2/4; S13:= - S11; S14:= S12;
        S24:= Q2/24; S22:= Q1/6 + S24; S23:= - S12;
        S33:= S11; S34:= - S12; S44:= S24 + Q3/6;
  
        "COMMENT" ELEMENT MASS MATRIX;
        R1:= R3; R2:= R(X2); R3:= R(XL);
        M11:= (R1 + R2)/6; M12:= R2/24; M13:= R2/6; M14:= - M12;
        M22:= R2/96; M23:= - M14; M24:= - M22;
        M33:= (R2 + R3)/6; M34:= M14; M44:= M22;
  
        "COMMENT" ELEMENT LOAD VECTOR;
        F1:= F3; F2:= F(X2); F3:= F(XL);
        B1:= H*(F1 + 2*F2)/6; B3:= H*(F3 + 2*F2)/6;
        B2:= H2*F2/12; B4:= - B2;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 8
                                                                  ;
  
  
        A11:= B11/H3 + S11/H + M11*H; A12:= B12/H2 + S12 + M12*H2;
        A13:= B13/H3 + S13/H + M13*H; A14:= B14/H2 + S14 + M14*H2;
        A22:= B22/H + S22*H + M22*H3; A23:= B23/H2 + S23 + M23*H2;
        A24:= B24/H + S24*H + M24*H3; A34:= B34/H2 + S34 + M34*H2;
        A33:= B33/H3 + S33/H + M33*H; A44:= B44/H + S44*H + M44*H3
    "END" "ELSE" "IF"ORDER=6"THEN"
    "BEGIN" "OWN" "REAL"P4, Q4, R4, F4;
        "REAL"H, H2, H3, X2, X3,
        P1, P2, P3, Q1, Q2, Q3,
        R1, R2, R3, F1, F2, F3,
        B11, B12, B13, B14, B15, B22, B23, B24, B25,
        B33, B34, B35, B44, B45, B55,
        S11, S12, S13, S14, S15, S22, S23, S24, S25,
        S33, S34, S35, S44, S45, S55,
        M11, M12, M13, M14, M15, M22, M23, M24, M25,
        M33, M34, M35, M44, M45, M55,
        A15, A25, A35, A45, A55, C1, C2, C3, C4, B5;
        "IF"L=1"THEN"
        "BEGIN"P4:= P(XL1); Q4:= Q(XL1); R4:= R(XL1); F4:= F(XL1)
        "END";
  
        H:= XL - XL1; H2:= H*H; H3:= H*H2;
        X2:= 0.27639320225*H + XL1; X3:= XL1 + XL - X2;
  
        "COMMENT" ELEMENT BENDING MATRIX;
        P1:= P4; P2:= P(X2); P3:= P(X3); P4:= P(XL);
        B11:= + 4.0333333333333"+1*P1 + 1.1124913866738"-1*P2
              + 1.4422084194664"+1*P3 + 8.3333333333333"+0*P4;
        B12:= + 1.4666666666667"+1*P1 - 3.3191425091659"-1*P2
              + 2.7985809175818"+0*P3 + 1.6666666666667"+0*P4;
        B13:= + 1.8333333333333"+1*(P1+P4)
              + 1.2666666666667"+0*(P2+P3);
        B15:= - (B11 + B13); B14:= - (B12 + B13 + B15/2);
        B22:= + 5.3333333333333"+0*P1 + 9.9027346441674"-1*P2
              + 5.4305986891624"-1*P3 + 3.3333333333333"-1*P4;
        B23:= + 6.6666666666667"+0*P1 - 3.7791278464167"+0*P2
              + 2.4579451308295"-1*P3 + 3.6666666666667"+0*P4;
        B25:= - (B12 + B23); B24:= - (B22 + B23 + B25/2);
        B33:= + 8.3333333333333"+0*P1 + 1.4422084194666"+1*P2
              + 1.1124913866726"-1*P3 + 4.0333333333333"+1*P4;
        B35:= - (B13 + B33); B34:= - (B23 + B33 + B35/2);
        B45:= - (B14 + B34); B44:= - (B24 + B34 + B45/2);
        B55:= - (B15 + B35);
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                      PAGE 9
                                                                  ;
  
  
        "COMMENT" ELEMENT STIFFNESS MATRIX;
        Q1:= Q4; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(XL);
        S11:= + 2.8844168389330"+0*Q2 + 2.2249827733448"-2*Q3;
        S12:= + 2.5671051872498"-1*Q2 + 3.2894812749994"-3*Q3;
        S13:= + 2.5333333333333"-1*(Q2+Q3);
        S14:= - 3.7453559925005"-2*Q2 - 2.2546440074988"-2*Q3;
        S15:= - (S13 + S11);
        S22:= + 8.3333333333333"-2*Q1 + 2.2847006554164"-2*Q2
              + 4.8632677916445"-4*Q3;
        S23:= + 2.2546440075002"-2*Q2 + 3.7453559924873"-2*Q3;
        S24:= - 3.3333333333333"-3*(Q2+Q3);
        S25:= - (S12 + S23);
        S33:= + 2.2249827733471"-2*Q2 + 2.8844168389330"+0*Q3;
        S34:= - 3.2894812750127"-3*Q2 - 2.5671051872496"-1*Q3;
        S35:= - (S13 + S33);
        S44:= + 4.8632677916788"-4*Q2
              + 2.2847006554161"-2*Q3 + 8.3333333333338"-2*Q4;
        S45:= - (S14 + S34);
        S55:= - (S15 + S35);
  
        "COMMENT" ELEMENT MASS MATRIX;
        R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
        M11:= + 8.3333333333333"-2*R1 + 1.0129076086083"-1*R2
              + 7.3759058058380"-3*R3;
        M12:= + 1.3296181273333"-2*R2 + 1.3704853933353"-3*R3;
        M13:= - 2.7333333333333"-2*(R2+R3);
        M14:= + 5.0786893258335"-3*R2 + 3.5879773408333"-3*R3;
        M15:= + 1.3147987115999"-1*R2 - 3.5479871159991"-2*R3;
        M22:= + 1.7453559925000"-3*R2 + 2.5464400750059"-4*R3;
        M23:= - 3.5879773408336"-3*R2 - 5.0786893258385"-3*R3;
        M24:= + 6.6666666666667"-4*(R2+R3);
        M25:= + 1.7259029213333"-2*R2 - 6.5923625466719"-3*R3;
        M33:= + 7.3759058058380"-3*R2
              + 1.0129076086083"-1*R3 + 8.3333333333333"-2*R4;
        M34:= - 1.3704853933333"-3*R2 - 1.3296181273333"-2*R3;
        M35:= - 3.5479871159992"-2*R2 + 1.3147987115999"-1*R3;
        M44:= + 2.5464400750008"-4*R2 + 1.7453559924997"-3*R3;
        M45:= + 6.5923625466656"-3*R2 - 1.7259029213330"-2*R3;
        M55:= + .17066666666667"+0*(R2+R3);
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                     PAGE 10
                                                                 ;
  
  
        "COMMENT" ELEMENT LOAD VECTOR;
        F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
        B1:= + 8.3333333333333"-2*F1 + 2.0543729868749"-1*F2
             - 5.5437298687489"-2*F3;
        B2:= + 2.6967233145832"-2*F2 - 1.0300566479175"-2*F3;
        B3:= - 5.5437298687489"-2*F2
             + 2.0543729868749"-1*F3 + 8.3333333333333"-2*F4;
        B4:= + 1.0300566479165"-2*F2 - 2.6967233145830"-2*F3;
        B5:= + 2.6666666666667"-1*(F2+F3);
  
        A11:= H2*(H2*M11 + S11) + B11; A12:= H2*(H2*M12 + S12) + B12;
        A13:= H2*(H2*M13 + S13) + B13; A14:= H2*(H2*M14 + S14) + B14;
        A15:= H2*(H2*M15 + S15) + B15; A22:= H2*(H2*M22 + S22) + B22;
        A23:= H2*(H2*M23 + S23) + B23; A24:= H2*(H2*M24 + S24) + B24;
        A25:= H2*(H2*M25 + S25) + B25; A33:= H2*(H2*M33 + S33) + B33;
        A34:= H2*(H2*M34 + S34) + B34; A35:= H2*(H2*M35 + S35) + B35;
        A44:= H2*(H2*M44 + S44) + B44; A45:= H2*(H2*M45 + S45) + B45;
        A55:= H2*(H2*M55 + S55) + B55;
  
        "COMMENT" STATIC CONDENSATION;
        C1:= A15/A55; C2:= A25/A55; C3:= A35/A55; C4:= A45/A55;
        B1:= (B1 - C1*B5)*H; B2:= (B2 - C2*B5)*H2;
        B3:= (B3 - C3*B5)*H; B4:= (B4 - C4*B5)*H2;
        A11:= (A11 - C1*A15)/H3; A12:= (A12 - C1*A25)/H2;
        A13:= (A13 - C1*A35)/H3; A14:= (A14 - C1*A45)/H2;
        A22:= (A22 - C2*A25)/H; A23:= (A23 - C2*A35)/H2;
        A24:= (A24 - C2*A45)/H; A33:= (A33 - C3*A35)/H3;
        A34:= (A34 - C3*A45)/H2; A44:= (A44 - C4*A45)/H;
    "END" "ELSE"
    "BEGIN" "OWN" "REAL"P5, Q5, R5, F5;
        "REAL" X2, X3, X4, H, H2, H3,
        P1, P2, P3, P4, Q1, Q2, Q3, Q4,
        R1, R2, R3, R4, F1, F2, F3, F4,
        B11, B12, B13, B14, B15, B16, B22, B23, B24, B25, B26,
        B33, B34, B35, B36, B44, B45, B46, B55, B56, B66,
        S11, S12, S13, S14, S15, S16, S22, S23, S24, S25, S26,
        S33, S34, S35, S36, S44, S45, S46, S55, S56, S66,
        M11, M12, M13, M14, M15, M16, M22, M23, M24, M25, M26,
        M33, M34, M35, M36, M44, M45, M46, M55, M56, M66,
        C15, C16, C25, C26, C35, C36, C45, C46, B5, B6,
        A15, A16, A25, A26, A35, A36, A45, A46, A55, A56, A66, DET;
        "IF"L=1"THEN"
        "BEGIN"P5:= P(XL1); Q5:= Q(XL1); R5:= R(XL1); F5:= F(XL1)
        "END";
        H:= XL - XL1; H2:= H*H; H3:= H*H2;
        X2:= XL1 + H*.172673164646; X3:= XL1 + H/2; X4:= XL1 + XL - X2;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                     PAGE 11
                                                                 ;
  
  
        "COMMENT" ELEMENT BENDING MATRIX;
        P1:= P5; P2:= P(X2); P3:= P(X3); P4:= P(X4); P5:= P(XL);
        B11:= + 105.8*P1 + 9.8*P5     + 7.3593121303513"-2*P2
              + 2.2755555555556"+1*P3 + 7.0565656088553"+0*P4;
        B12:= + 27.6*P1  + 1.4*P5     - 3.41554824811"-1*P2
              + 2.8444444444444"+0*P3 + 1.0113960946522"+0*P4;
        B13:= - 32.2*(P1 + P5)        - 7.2063492063505"-1*(P2 + P4)
              + 2.2755555555556"+1*P3;
        B14:= + 4.6*P1 + 8.4*P5       + 1.0328641222944"-1*P2
              - 2.8444444444444"+0*P3 - 3.3445562534992"+0*P4;
        B15:= - (B11 + B13); B16:= - (B12 + B13 + B14 + B15/2);
        B22:= + 7.2*P1 + 0.2*P5       + 1.5851984028581"+0*P2
              + 3.5555555555556"-1*P3 + 1.4496032730059"-1*P4;
        B23:= - 8.4*P1 - 4.6*P5       + 3.3445562534992"+0*P2
              + 2.8444444444444"+0*P3 - 1.0328641222944"-1*P4;
        B24:= + 1.2*(P1 + P5)         - 4.7936507936508"-1*(P2 + P4)
              - 3.5555555555556"-1*P3;
        B25:= - (B12 + B23); B26:= - (B22 + B23 + B24 + B25/2);
        B33:= + 7.0565656088553"+0*P2 + 2.2755555555556"+1*P3
              + 7.3593121303513"-2*P4 + 105.8*P5 + 9.8*P1;
        B34:= - 1.4*P1 - 27.6*P5      - 1.0113960946522"+0*P2
              - 2.8444444444444"+0*P3 + 3.4155482481100"-1*P4;
        B35:= - (B13 + B33); B36:= - (B23 + B33 + B34 + B35/2);
        B44:= +7.2*P5 + P1/5          + 1.4496032730059"-1*P2
              + 3.5555555555556"-1*P3 + 1.5851984028581"+0*P4;
        B45:= - (B14 + B34); B46:= - (B24 + B34 + B44 + B45/2);
        B55:= - (B15 + B35); B56:= - (B16 + B36);
        B66:= - (B26 + B36 + B46 + B56/2);
  
        "COMMENT" ELEMENT STIFFNESS MATRIX;
        Q1:= Q5; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(X4); Q5:= Q(XL);
        S11:= + 3.0242424037951"+0*Q2 + 3.1539909130065"-2*Q4;
        S12:= + 1.2575525581744"-1*Q2 + 4.1767169716742"-3*Q4;
        S13:= - 3.0884353741496"-1*(Q2+Q4);
        S14:= + 4.0899041243062"-2*Q2 + 1.2842455355577"-2*Q4;
        S15:= - (S13 + S11);
        S16:= + 5.9254861177068"-1*Q2 + 6.0512612719116"-2*Q4;
        S22:= + 5.2292052865422"-3*Q2 + 5.5310763862796"-4*Q4 + Q1/20;
        S23:= - 1.2842455355577"-2*Q2 - 4.0899041243062"-2*Q4;
        S24:= + 1.7006802721088"-3*(Q2+Q4);
        S25:= - (S12 + S23);
        S26:= + 2.4639593097426"-2*Q2 + 8.0134681270641"-3*Q4;
        S33:= + 3.1539909130065"-2*Q2 + 3.0242424037951"+0*Q4;
        S34:= - 4.1767169716742"-3*Q2 - 1.2575525581744"-1*Q4;
        S35:= - (S13 + S33);
        S36:= - 6.0512612719116"-2*Q2 - 5.9254861177068"-1*Q4;
        S44:= + 5.5310763862796"-4*Q2 + 5.2292052865422"-3*Q4 + Q5/20;
        S45:= - (S14 + S34);
        S46:= + 8.0134681270641"-3*Q2 + 2.4639593097426"-2*Q4;
        S55:= - (S15 + S35); S56:= -(S16 + S36);
        S66:= + 1.1609977324263"-1*(Q2+Q4) + 3.5555555555556"-1*Q3;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                     PAGE 12
                                                                 ;
  
  
        "COMMENT" ELEMENT MASS MATRIX;
        R1:= R5; R2:= R(X2); R3:= R(X3); R4:= R(X4); R5:= R(XL);
        M11:= + 9.7107020727310"-2*R2 + 1.5810259199180"-3*R4 + R1/20;
        M12:= + 8.2354889460254"-3*R2 + 2.1932154960071"-4*R4;
        M13:= + 1.2390670553936"-2*(R2+R4);
        M14:= - 1.7188466249968"-3*R2 - 1.0508326752939"-3*R4;
        M15:= + 5.3089789712119"-2*R2 + 6.7741558661060"-3*R4;
        M16:= - 1.7377712856076"-2*R2 + 2.2173630018466"-3*R4;
        M22:= + 6.9843846173145"-4*R2 + 3.0424512029349"-5*R4;
        M23:= + 1.0508326752947"-3*R2 + 1.7188466249936"-3*R4;
        M24:= - 1.4577259475206"-4*(R2+R4);
        M25:= + 4.5024589679127"-3*R2 + 9.3971790283374"-4*R4;
        M26:= - 1.4737756452780"-3*R2 + 3.0759488725998"-4*R4;
        M33:= + 1.5810259199209"-3*R2 + 9.7107020727290"-2*R4 + R5/20;
        M34:= - 2.1932154960131"-4*R2 - 8.2354889460254"-3*R4;
        M35:= + 6.7741558661123"-3*R2 + 5.3089789712112"-2*R4;
        M36:= - 2.2173630018492"-3*R2 + 1.7377712856071"-2*R4;
        M44:= + 3.0424512029457"-5*R2 + 6.9843846173158"-4*R4;
        M45:= - 9.3971790283542"-4*R2 - 4.5024589679131"-3*R4;
        M46:= + 3.0759488726060"-4*R2 - 1.4737756452778"-3*R4;
        M55:= + 2.9024943310657"-2*(R2+R4) + 3.5555555555556"-1*R3;
        M56:= + 9.5006428402050"-3*(R4-R2);
        M66:= + 3.1098153547125"-3*(R2+R4);
  
        "COMMENT" ELEMENT LOAD VECTOR;
        F1:= F5; F2:= F(X2); F3:= F(X3); F4:= F(X4); F5:= F(XL);
        B1:= + 1.6258748099336"-1*F2 + 2.0745852339969"-2*F4 + F1/20;
        B2:= + 1.3788780589233"-2*F2 + 2.8778860774335"-3*F4;
        B3:= + 2.0745852339969"-2*F2 + 1.6258748099336"-1*F4 + F5/20;
        B4:= - 2.8778860774335"-3*F2 - 1.3788780589233"-2*F4;
        B5:= + (F2 + F4)/11.25       + 3.5555555555556"-1*F3;
        B6:= + 2.9095718698132"-2*(F4-F2);
  
        A11:= H2*(H2*M11 + S11) + B11; A12:= H2*(H2*M12 + S12) + B12;
        A13:= H2*(H2*M13 + S13) + B13; A14:= H2*(H2*M14 + S14) + B14;
        A15:= H2*(H2*M15 + S15) + B15; A16:= H2*(H2*M16 + S16) + B16;
        A22:= H2*(H2*M22 + S22) + B22; A23:= H2*(H2*M23 + S23) + B23;
        A24:= H2*(H2*M24 + S24) + B24; A25:= H2*(H2*M25 + S25) + B25;
        A26:= H2*(H2*M26 + S26) + B26; A33:= H2*(H2*M33 + S33) + B33;
        A34:= H2*(H2*M34 + S34) + B34; A35:= H2*(H2*M35 + S35) + B35;
        A36:= H2*(H2*M36 + S36) + B36; A44:= H2*(H2*M44 + S44) + B44;
        A45:= H2*(H2*M45 + S45) + B45; A46:= H2*(H2*M46 + S46) + B46;
        A55:= H2*(H2*M55 + S55) + B55; A56:= H2*(H2*M56 + S56) + B56;
        A66:= H2*(H2*M66 + S66) + B66;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.2.2.1   (JANUARY 1976)                     PAGE 13
                                                                 ;
  
  
        "COMMENT" STATIC CONDENSATION;
        DET:= - A55*A66 + A56*A56;
        C15:= (A15*A66 - A16*A56)/DET; C16:= (A16*A55 - A15*A56)/DET;
        C25:= (A25*A66 - A26*A56)/DET; C26:= (A26*A55 - A25*A56)/DET;
        C35:= (A35*A66 - A36*A56)/DET; C36:= (A36*A55 - A35*A56)/DET;
        C45:= (A45*A66 - A46*A56)/DET; C46:= (A46*A55 - A45*A56)/DET;
        A11:= (A11 + C15*A15 + C16*A16)/H3;
        A12:= (A12 + C15*A25 + C16*A26)/H2;
        A13:= (A13 + C15*A35 + C16*A36)/H3;
        A14:= (A14 + C15*A45 + C16*A46)/H2;
        A22:= (A22 + C25*A25 + C26*A26)/H;
        A23:= (A23 + C25*A35 + C26*A36)/H2;
        A24:= (A24 + C25*A45 + C26*A46)/H;
        A33:= (A33 + C35*A35 + C36*A36)/H3;
        A34:= (A34 + C35*A45 + C36*A46)/H2;
        A44:= (A44 + C45*A45 + C46*A46)/H;
        B1:= (B1 + C15*B5 + C16*B6)*H; B2:= (B2 + C25*B5 + C26*B6)*H2;
        B3:= (B3 + C35*B5 + C36*B6)*H; B4:= (B4 + C45*B5 + C46*B6)*H2;
    "END"EL.MATVECEVAL.;
  
  
    L:= 1; W:= V:= 0; N2:= N + N - 2; XL1:= X[0]; XL:= X[1];
    YA:= E[1]; ZA:= E[2]; YB:= E[3]; ZB:= E[4];
    ELEMENTMATVECEVALUATION; EM[2]:= "-12;
    R1:= B3 - A13*YA - A23*ZA; D1:= A33; D2:= A44;
    R2:= B4 - A14*YA - A24*ZA; E1:= A34;
    "FOR"L:= L + 1"WHILE"L<N"DO"
    "BEGIN" XL1:= XL; XL:= X[L];
      ELEMENTMATVECEVALUATION;
      A[W + 1]:= D1 + A11; A[W + 4]:= E1 + A12;
      A[W + 7]:= A13; A[W + 10]:= A14;
      A[W + 5]:= D2 + A22; A[W + 8]:= A23;
      A[W + 11]:= A24; A[W + 14]:= 0;
      Y[V + 1]:= R1 + B1; Y[V + 2]:= R2 + B2;
      R1:= B3; R2:= B4; V:= V + 2; W:= W + 8;
      D1:= A33; D2:= A44; E1:= A34
    "END";
    L:= N; XL1:= XL; XL:= X[L]; ELEMENTMATVECEVALUATION;
    Y[N2      - 1]:= R1 + B1 - A13*YB - A14*ZB;
    Y[N2]:= R2 + B2 - A23*YB - A24*ZB;
      A[W + 1]:= D1 + A11; A[W + 4]:= E1 + A12; A[W + 5]:= D2 + A22;
    CHLDECSOLBND(A, N2, 3, EM, Y)
  "END" FEMHERM;
        "EOP"
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 1      
        
        
        
 AUTHOR: J.J.G. ADMIRAAL.   
        
        
 INSTITUTE: UNIVERSITY OF AMSTERDAM.  
        
        
 RECEIVED: 751101.
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS TWO MAIN AND NINE AUXILIARY PROCEDURES:  
        
     THE TWO MAIN PROCEDURES ARE:     
     A. QZIVAL FINDS N PAIRS OF SCALARS (ALFA[M],BETA[M]),WHERE BETA[M]       
        IS REAL,SUCH THAT THE MATRIX BETA[M] * A - ALFA[M] * B IS   
        SINGULAR. 
     B. QZI FINDS N PAIRS OF SCALARS (ALFA[M], BETA[M]),WHERE BETA[M] IS      
        REAL,SUCH THAT THE MATRIX BETA[M] * A - ALFA[M] * B IS SINGULAR;      
        MOREOVER THE GENERALIZED EIGENVECTORS (THE HOMOGENOUS SOLUTION
        OF ( BETA[M] * A - ALFA[M] * B ) * X = 0 )  ARE CALCULATED. 
        
     THE AUXILIARY PROCEDURES ARE:    
     A. HSHDECMUL:
        THIS PROCEDURE CALCULATES REAL MATRICES Q AND R SUCH
        THAT Q.A=R WHERE A IS A GIVEN REAL SQUARE MATRIX, Q IS A    
        PRODUCT OF HOUSEHOLDER MATRICES AND R AN UPPERTRIANGULAR MATRIX.      
        MOREOVER Q.B IS FORMED WITH B, A GIVEN MATRIX.    
     B. HESTGL3:  
        GIVEN THE REAL SQUARE MATRICES A,B AND X, WITH B AN 
        UPPER TRIANGULAR MATRIX,HESTGL3 CALCULATES THE MATRICES Q,Z,H,R,      
        WHERE Q,Z ARE ORTHOGONAL, H UPPER HESSENBERG AND R AN UPPER 
        TRIANGULAR MATRIX SUCH THAT Q.A.Z = H AND Q.B.Z = R.
        FURTHER: A:= Q.A.Z ; B:= Q.B.Z AND X:= Q.X.Z.     
     C. HESTGL2:  
        SEE HESTGL3,BUT HERE THE MATRIX X HAS BEEN LEFT OUT.
     D. HSH2COL:  
        THIS PROCEDURE CALCULATES A HOUSEHOLDER MATRIX Q  
        SUCH THAT BY PREMULTIPLYING A GIVEN COLUMN VECTOR V BY Q    
        A ZERO ELEMENT IS FORMED IN V.
        HERE THE VECTOR V IS A COLUMN OF A MATRIX.
        FURTHER: A:= Q.A AND B:= Q.B  
        WHERE A,B ARE TWO GIVEN REAL MATRICES.  
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 2      
        
        
        
     E. HSH3COL:  
        THIS PROCEDURE CALCULATES A HOUSEHOLDER MATRIX    
        Q SUCH THAT BY PREMULTIPLYING A GIVEN COLUMN VECTOR V BY Q TWO
        SUCCESSIVE ZERO ELEMENTS ARE FORMED IN V. 
        HERE THE VECTOR V IS A COLUMN OF A MATRIX.
        FURTHER: A:= Q.A AND B:= Q.B  
        WHERE A AND B ARE TWO GIVEN REAL MATRICES.
     F. HSH2ROW3: 
        THIS PROCEDURE CALCULATES A HOUSEHOLDER MATRIX    
        Z SUCH THAT BY POSTMULTIPLYING A GIVEN ROWVECTOR  
        V BY Z A ZERO ELEMENT IS FORMED IN V.   
        HERE THE VECTOR V IS A ROW OF A MATRIX. 
        FURTHER: A:= A.Z; B:= B.Z AND X:= X.Z   
        WHERE A,B,X ARE THREE GIVEN REAL MATRICES.
     G. HSH2ROW2: 
        SEE HSH2ROW3, BUT HERE THE MATRIX X HAS BEEN LEFT OUT.      
     H. HSH3ROW3: 
        THIS PROCEDURE CALCULATES A HOUSEHOLDER MATRIX    
        Z SUCH THAT BY POSTMULTIPLYING A GIVEN ROWVECTOR  
        V BY Z, TWO SUCCESSIVE ZERO ELEMENTS ARE FORMED   
        IN V. HERE THE VECTOR V IS A ROW OF A MATRIX.     
        FURTHER: A:= A.Z; B:= B.Z AND X:= X.Z   
        WHERE A,B AND X ARE THREE GIVEN REAL MATRICES.    
     I: HSH3ROW2: 
        SEE HSH3ROW3, BUT HERE THE MATRIX X HAS BEEN LEFT OUT.      
        
 KEYWORDS:
           HOUSEHOLDER'S TRANSFORMATION,
           GENERALIZED EIGENVALUES,   
           GENERALIZED EIGENVECTORS,  
           UPPER HESSENBERG MATRIX,   
           UPPER TRIANGULAR MATRIX.   
        
 REFERENCES:      
        
     [1]. C.B. MOLER AND G.W. STEWART.
          AN ALGORITHM FOR THE GENERALIZED MATRIX EIGENVALUE
          PROBLEM  A * X = LAMBDA * B * X.      
          REPORT STANFORD UNIVERSITY  
          STAN-CS-232-71;   
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 3      
        
        
        
 SUBSECTION: QZIVAL 
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" QZIVAL(N,A,B,ALFR,ALFI,BETA,ITER,EM);    
     "VALUE" N; "INTEGER" N; "ARRAY" A,B,ALFR,ALFI,BETA,EM; 
     "INTEGER" "ARRAY" ITER;
     "CODE" 34600;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:    <ARITHMETIC EXPRESSION>;   
           THE NUMBER OF ROWS AND COLUMNS OF THE MATRICES A,B       
     A:    <ARRAY IDENTIFIER>;
           "ARRAY" A[1:N,1:N];
           ENTRY: THE GIVEN MATRIX;   
           EXIT: A QUASI UPPER-TRIANGULAR MATRIX
                 (SEE METHOD AND PERFORMANCE);  
     B:    <ARRAY IDENTIFIER>;
           "ARRAY" B[1:N,1:N];
           ENTRY: THE GIVEN MATRIX;   
           EXIT: AN UPPER-TRIANGULAR MATRIX;    
     ALFR: <ARRAY IDENTIFIER>;
           "ARRAY" ALFR[1:N]; 
           EXIT : THE REAL PARTS OF ALFA[1:N]   
           (SEE METHOD AND PERFORMANCE);
     ALFI: <ARRAY IDENTIFIER>;
           "ARRAY" ALFI[1:N]; 
           EXIT : THE IMAGINARY PARTS OF ALFA[1:N];       
     BETA: <ARRAY IDENTIFIER>;
           "ARRAY" BETA[1:N]; 
           EXIT : THE REAL SCALARS BETA[N]      
     ITER: <ARRAY IDENTIFIER>;
           "INTEGER" "ARRAY" ITER[1:N]; 
           TROUBLE INDICATOR AND ITERATION COUNTER;       
           IF ITER[1]=0 THEN NO TROUBLE IS SIGNALIZED,    
           FURTHER SEE METHOD AND PERFORMANCE;  
     EM:   <ARRAY IDENTIFIER>;
           "ARRAY" EM[0:1]; 
            ENTRY: EM[0]: THE SMALLEST POSITIVE MACHINE NUMBER;     
                   EM[1]: THE RELATIVE PRECISION OF ELEMENTS OF A AND B;      
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 4      
        
        
        
 PROCEDURES USED: 
     TAMMAT    = CP 34014   
     ELMCOL    = CP 34023   
     HSHDECMUL = CP 34602   
     HESTGL2   = CP 34604   
     HSH2COL   = CP 34605   
     HSH3COL   = CP 34606   
     HSH2ROW2  = CP 34608   
     HSH3ROW2  = CP 34610   
     CHSH2     = CP 34611   
     HSHVECMAT = CP 31070   
     HSHVECTAM = CP 31073   
        
 REQUIRED CENTRAL MEMORY : IN HSHDECMUL, AN ARRAY OF N REALS IS DECLARED.     
        
 RUNNING TIME: PROPORTIONAL TO N ** 3 
        
 METHOD AND PERFORMANCE;    
     THE PROCEDURE QZIVAL SOLVES THE GENERALIZED MATRIX   
     EIGENVALUE PROBLEM A * X = LAMBDA * B * X BY MEANS   
     OF QZ ITERATION (SEE REF[1]);    
     QZIVAL FINDS N PAIRS OF SCALARS (ALFA[M],BETA[M])    
     SUCH THAT BETA[M] * A - ALFA[M] * B IS SINGULAR.     
     THE EIGENVALUES OF A * X - LAMBDA * B * X CAN BE OBTAINED      
     BY DIVIDING ALFA[M] BY BETA[M],EXCEPT BETA[M] MIGHT BE ZERO.   
     IN THIS ALGORITHM ONLY UNITARY TRANSFORMATIONS ARE   
     APPLIED; A FORTIORI NO INVERSES ARE CALCULATED, SO   
     EITHER A OR B (OR BOTH) MAY BE SINGULAR.   
     BETA[M] IS REAL, ALFA[M] IS COMPLEX.       
     REAL AND IMAGINARY PARTS ARE GIVEN IN ALFR[M] AND ALFI[M].     
     THE OCCURRENCE OF COMPLEX PAIRS IS ALWAYS IN 
     SUCCESSIVE ELEMENTS, SUCH THAT ALFA[M]/BETA[M] AND   
     ALFA[M+1]/BETA[M+1] ARE COMPLEX CONJUGATE, BUT ALFA[M] 
     AND ALFA[M+1] ARE NOT NECESSARILY CONJUGATE. 
     ONLY REAL ARITHMETIC IS USED IN THE PROCEDURE.       
     IF A AND B WERE REDUCED TO TRIANGULAR FORM BY UNITARY
     TRANSFORMATIONS,ALFA AND BETA WOULD BE THE DIAGONALS.
     A AND B ARE ACTUALLY REDUCED TO QUASI-TRIANGULAR FORM HAVING ONLY
     1-BY-1 AND 2-BY-2 BLOCKS ON THE DIAGONAL OF A.       
     IF ALFA[M] IS NOT REAL, THEN BETA[M] IS NOT ZERO.    
     ITER IS THE TROUBLE INDICATOR AND ITERATION COUNTER. 
     IF ITER[1]=0 THEN EVERYTHING IS O.K.       
     ITER[M] IS THE NUMBER OF ITERATIONS NEEDED FOR THE M-TH EIGENVALUE.      
     IF ITER[1] THROUGH ITER[M]= -1 THEN THE ITERATION FOR THE M-TH 
     EIGENVALUE DID NOT CONVERGE AND ALFA[1] THROUGH ALFA[M] AND BETA[1]      
     THROUGH BETA[M] ARE PROBABLY INACCURATE.   
        
 EXAMPLE OF USE:  
        
  "BEGIN" "ARRAY" A,B[1:4,1:4],ALFR,ALFI,BETA[1:4],EM[0:1]; 
      "INTEGER" "ARRAY" ITER[1:4];"INTEGER" K,L;
      "PROCEDURE" QZIVAL(N,A,B,ALFR,ALFI,BETA,ITER,EM);"CODE" 34600;
      A[1,1]:=2; A[1,2]:=3; A[1,3]:=-3; A[1,4]:=4;
1SECTION : 3.4.1.2           (JANUARY 1976)                      PAGE 5      
        
        
        
      A[2,1]:=1; A[2,2]:=-1; A[2,3]:=5; A[2,4]:=1;
      A[3,1]:=0; A[3,2]:=2; A[3,3]:=6; A[3,4]:=8; 
      A[4,1]:=1; A[4,2]:=1; A[4,3]:=0; A[4,4]:=4; 
      B[1,1]:=1; B[1,2]:=5; B[1,3]:=9; B[1,4]:=0; 
      B[2,1]:=2; B[2,2]:=6; B[2,3]:=10; B[2,4]:=2;
      B[3,1]:=3; B[3,2]:=7; B[3,3]:=11; B[3,4]:=-1;       
      B[4,1]:=4; B[4,2]:=8; B[4,3]:=12; B[4,4]:=3;
      OUTPUT(61,"(""("A")",/,4(4(+ZDBB),/),/")",A);       
      OUTPUT(61,"(""("B")",/,4(4(+ZDBB),/),/")",B);       
      EM[0]:="-294;EM[1]:="-15;       
      QZIVAL(4,A,B,ALFR,ALFI,BETA,ITER,EM);     
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      OUTPUT(61,"(""("ITER[")",D,"("]=")",ZD,/")",K,ITER[K]);       
      OUTPUT(61,"(""("ALFA(REAL PART)")"8B,"("ALFA(IMAGINARY PART)")" 
      3B,"("BETA")",/")");  
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      OUTPUT(61,"("3(N),/")",ALFR[K],ALFI[K],BETA[K]);    
      OUTPUT(61,"("/"("LAMBDA(REAL PART)")"6B,  
      "("LAMBDA(IMAGINARY PART)")"/")");
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      "BEGIN" "IF" BETA[K]=0 "THEN"   
          OUTPUT(61,"(""("INFINITE")"15B,"("INDEFINITE")"/,")")     
          "ELSE" OUTPUT(61,"("2(N),/")",ALFR[K]/BETA[K],ALFI[K]/BETA[K])      
      "END"       
  "END" 
        
  A     
   +2   +3   -3   +4
   +1   -1   +5   +1
   +0   +2   +6   +8
   +1   +1   +0   +4
        
  B     
   +1   +5   +9   +0
   +2   +6  +10   +2
   +3   +7  +11   -1
   +4   +8  +12   +3
        
  ITER[1]= 0      
  ITER[2]= 0      
  ITER[3]= 0      
  ITER[4]= 5      
        
  ALFA(REAL PART)        ALFA(IMAGINARY PART)   BETA      
  -4.4347115652167"+000  +0.0000000000000"+000  +0.0000000000000"+000 
  -5.7288406521003"+000  +0.0000000000000"+000  +2.8441121744896"+000 
  -8.6671777386054"-001  +2.7607904944916"+000  +8.7617886336960"+000 
  -4.7262205157527"-001  -1.5054617625576"+000  +4.7778119295757"+000 
        
  LAMBDA(REAL PART)      LAMBDA(IMAGINARY PART) 
  INFINITE               INDEFINITE   
  -2.0142808372628"+000  +0.0000000000000"+000  
  -9.8920187429234"-002  +3.1509439566644"-001  
  -9.8920187429236"-002  -3.1509439566645"-001  
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 6      
        
        
        
 SUBSECTION: QZI  
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" QZI(N,A,B,X,ALFR,ALFI,BETA,ITER,EM);     
     "VALUE" N; "INTEGER" N; "ARRAY" A,B,X,ALFR,ALFI,BETA,EM;       
     "INTEGER" "ARRAY" ITER;
     "CODE" 34601;
        
     THE MEANING OF THE FORMAL PARAMETERS IS;   
     N:    <ARITHMETIC EXPRESSION>;   
           THE NUMBER OF ROWS AND COLUMNS OF THE MATRICES A,B AND X;
     A:    <ARRAY IDENTIFIER> 
           "ARRAY" A[1:N,1:N];
           ENTRY: THE GIVEN MATRIX A; 
           EXIT: A QUASI UPPER TRIANGULAR MATRIX; 
                 (SEE METHOD AND PERFORMANCE);  
     B:    <ARRAY IDENTIFIER>;
           "ARRAY" B[1:N,1:N];
           ENTRY: THE GIVEN MATRIX B; 
           EXIT: AN UPPER-TRIANGULAR MATRIX;    
     X:    <ARRAY IDENTIFIER>;
           "ARRAY" X[1:N,1:N];
           ENTR: THE N*N UNIT MATRIX; 
           EXIT: THE MATRIX OF EIGENVECTORS,    
           THE EIGENVECTORS ARE STORED IN THE ARRAY X AS FOLLOWS:   
           IF ALFI[M]=0 THEN X[.,M] IS THE M-TH REAL EIGENVECTOR;   
           OTHERWISE, FOR EACH PAIR OF CONSECUTIVE COLUMNS
           X[.,M] AND X[.,M+1] ARE THE REAL     
           AND IMAGINARY PARTS OF THE M-TH COMPLEX EIGENVECTOR.     
           X[.,M] AND -X[.,M+1] ARE THE REAL AND IMAGINARY PARTS    
           OF THE M+1 -ST COMPLEX EIGENVECTOR.  
           THE EIGENVECTORS ARE NORMALIZED SUCH THAT THE LARGEST    
           COMPONENT IS 1 OR 1 + 0 * I. 
     ALFR: <ARRAY IDENTIFIER>;
           "ARRAY" ALFR[1:N]; 
           EXIT: THE REAL PARTS OF ALFA[1:N];   
     ALFI: <ARRAY IDENTIFIER>;
           "ARRAY" ALFI[1:N]; 
           EXIT: THE IMAGINARY PARTS OF ALFA[1:N];
     BETA: <ARRAY IDENTIFIER>;
           "ARRAY" BETA[1:N]; 
     ITER: <ARRAY IDENTIFIER>;
           "INTEGER" "ARRAY" ITER[1:N]; 
           TROUBLE INDICATOR AND ITERATION COUNTER;       
           IF ITER[1]=0 THEN NO TROUBLE IS SIGNALIZED,    
           FOR FURTHER INFORMATION SEE
           METHOD AND PEFORMANCE OF PROCEDURE QZIVAL (THIS SECTION).
     EM:   <ARRAY IDENTIFIER>;
           "ARRAY" EM[0:1]; 
           ENTRY: EM[0]: THE SMALLEST POSITIVE MACHINE NUMBER;      
                  EM[1]: THE RELATIVE PRECISION OF ELEMENTS OF A AND B;       
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 7      
        
        
        
 PROCEDURES USED: 
        
     MATMAT    = CP 34013   
     TAMMAT    = CP 34014   
     ELMCOL    = CP 34023   
     HSHDECMUL = CP 34602   
     HESTGL3   = CP 34603   
     HSH2COL   = CP 34605   
     HSH2ROW3  = CP 34607   
     HSH3ROW3  = CP 34609   
     HSH3COL   = CP 34606   
     CHSH2     = CP 34611   
     COMDIV    = CP 34342   
     HSHVECMAT = CP 31070   
     HSHVECTAM = CP 31073   
        
 RUNNING TIME: PROPORTIONAL TO N ** 3;
        
        
 REQUIRED CENTRAL MEMORY : IN HSHDECMUL, AN ARRAY OF N REALS IS DECLARED.     
        
        
 METHOD AND PERFORMANCE;    
        
     THE PROCEDURE QZI APPLIES THE SAME METHOD AS QZIVAL. 
        
        
 EXAMPLE OF USE:  
        
  "BEGIN" "ARRAY" A,B,X[1:4,1:4],ALFR,ALFI,BETA[1:4],EM[0:1];       
      "INTEGER" "ARRAY" ITER[1:4];"INTEGER" K,L;
      "PROCEDURE" QZI(N,A,B,X,ALFR,ALFI,BETA,ITER,EM);"CODE" 34601; 
      A[1,1]:=2; A[1,2]:=3; A[1,3]:=-3; A[1,4]:=4;
      A[2,1]:=1; A[2,2]:=-1; A[2,3]:=5; A[2,4]:=1;
      A[3,1]:=0; A[3,2]:=2; A[3,3]:=6; A[3,4]:=8; 
      A[4,1]:=1; A[4,2]:=1; A[4,3]:=0; A[4,4]:=4; 
      B[1,1]:=1; B[1,2]:=5; B[1,3]:=9; B[1,4]:=0; 
      B[2,1]:=2; B[2,2]:=6; B[2,3]:=10; B[2,4]:=2;
      B[3,1]:=3; B[3,2]:=7; B[3,3]:=11; B[3,4]:=-1;       
      B[4,1]:=4; B[4,2]:=8; B[4,3]:=12; B[4,4]:=3;
      "FOR" K:=1,2,3,4 "DO" "FOR" L:=1,2,3,4 "DO" 
      X[K,L]:="IF" K=L "THEN" 1 "ELSE" 0;       
      OUTPUT(61,"(""("A")",/,4(4(+ZDBB),/),/")",A);       
      OUTPUT(61,"(""("B")",/,4(4(+ZDBB),/),/")",B);       
      EM[0]:="-294;EM[1]:="-15;       
      QZI(4,A,B,X,ALFR,ALFI,BETA,ITER,EM);      
1SECTION : 3.4.1.2           (JANUARY 1976)                      PAGE 8      
        
        
        
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      OUTPUT(61,"(""("ITER[")",D,"("]=")",ZD,/")",K,ITER[K]);       
      OUTPUT(61,"("/"("EIGENVECTORS")",/,4(4(+D.8D"+2D2B),/),/")",X); 
      OUTPUT(61,"(""("ALFA(REAL PART)")"8B,"("ALFA(IMAGINARY PART)")" 
      9B,"("BETA")",/")");  
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      OUTPUT(61,"("3(N),/")",ALFR[K],ALFI[K],BETA[K]);    
      OUTPUT(61,"("/"("LAMBDA(REAL PART)")"6B,  
      "("LAMBDA(IMAGINARY PART)")"/")");
      "FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
      "BEGIN" "IF" BETA[K]=0 "THEN"   
          OUTPUT(61,"(""("INFINITE")"15B,"("INDEFINITE")"/")")      
          "ELSE"  
          OUTPUT(61,"("2(N),/")",ALFR[K]/BETA[K],ALFI[K]/BETA[K])   
      "END"       
  "END" 
        
  A     
   +2   +3   -3   +4
   +1   -1   +5   +1
   +0   +2   +6   +8
   +1   +1   +0   +4
        
  B     
   +1   +5   +9   +0
   +2   +6  +10   +2
   +3   +7  +11   -1
   +4   +8  +12   +3
        
  ITER[1]= 0      
  ITER[2]= 0      
  ITER[3]= 0      
  ITER[4]= 5      
        
  EIGENVECTORS    
  -5.00000000"-01  +1.00000000"+00  -6.29204867"-01  +6.52026261"-01
  +1.00000000"+00  -3.82541766"-02  +1.00000000"+00  +0.00000000"+00
  -5.00000000"-01  -3.04677732"-02  +1.65896051"-01  +1.09306265"-01
  -4.35116786"-15  -7.63328122"-01  -5.84845537"-01  +1.77430910"-01
        
  ALFA(REAL PART)        ALFA(IMAGINARY PART)   BETA      
  -4.4347115652167"+000  +0.0000000000000"+000  +0.0000000000000"+000 
  -5.7288406521003"+000  +0.0000000000000"+000  +2.8441121744896"+000 
  -8.6671777386054"-001  +2.7607904944916"+000  +8.7617886336960"+000 
  -4.7262205157527"-001  -1.5054617625576"+000  +4.7778119295757"+000 
        
  LAMBDA(REAL PART)      LAMBDA(IMAGINARY PART) 
  INFINITE               INDEFINITE   
  -2.0142808372628"+000  +0.0000000000000"+000  
  -9.8920187429234"-002  +3.1509439566644"-001  
  -9.8920187429236"-002  -3.1509439566645"-001  
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                     PAGE 9      
        
        
        
 SUBSECTION: HSHDECMUL.     
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" HSHDECMUL(N,A,B,DWARF);
     "VALUE" N,DWARF; "INTEGER"N; "REAL"DWARF; "ARRAY" A,B; 
     "CODE" 34602;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N: <ARITHMETIC EXPRESSION>;      
        THE ORDER OF THE GIVEN MATRICES;
     A: <ARRAY IDENTIFIER>; 
        "REAL" "ARRAY" A[1:N,1:N];    
        ENTRY: THE GIVEN MATRIX A;    
        EXIT: THE TRANSFORMED MATRIX Q.A (SEE BRIEF DESCRIPTION);   
     B: <ARRAY IDENTIFIER>; 
        "REAL" "ARRAY" B[1:N,1:N];    
        ENTRY: THE GIVEN MATRIX B;    
        EXIT: THE UPPER TRIANGULAR MATRIX Q.B (SEE BRIEF DESCRIPTION);
     DWARF: < ARITHMETIC EXPRESSION>; 
        THE SMALLEST POSITIVE MACHINE NUMBER.   
        
        
 PROCEDURES USED: 
     TAMMAT    = CP 34014;  
     HSHVECMAT = CP 31070.  
        
        
 REQUIRED CENTRAL MEMORY : AN ARRAY OF N REALS IS DECLARED. 
        
        
 METHOD AND PERFORMANCE:    
     FOR EACH MATRIX COLUMN A[I] A HOUSEHOLDERMATRIX Q IS FORMED,   
     SUCH THAT Q.A[I] HAS ONLY ZERO ELEMENTS BELOW THE DIAGONAL ELEMENT.      
     WHEN SUCCESSIVELY FOR I = 1,2,....,N-1 THESE TRANSFORMATIONS HAVE
     BEEN PERFORMED,THE MATRIX A HAS BEEN CHANGED INTO AN UPPER     
     TRIANGULAR MATRIX.     
     THE SAME TRANSFORMATIONS ARE PERFORMED ON THE MATRIX B 
        
        
 EXAMPLE OF USE:  
        
     THE PROCEDURE HSHDECMUL IS USED IN QZI AND QZIVAL (THIS SECTION).
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 10      
        
        
        
 SUBSECTION: HESTGL3:       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HESTGL3(N,A,B,X);    
     "VALUE" N; "INTEGER" N; "ARRAY" A,B,X;     
     "CODE" 34603;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N: <ARITHMETIC EXPRESSION>;      
        THE ORDER OF THE GIVEN MATRICES;
     A: <ARRAY IDENTIFIER>; 
        "ARRAY" A[1:N,1:N]; 
        ENTRY: THE GIVEN MATRIX A;    
        EXIT: THE UPPER HESSENBERG MATRIX Q.A.Z (SEE BRIEF DESCRIPTION);      
     B: <ARRAY IDENTIFIER>; 
        "ARRAY" B[1:N,1:N]; 
        ENTRY: THE GIVEN UPPER TRIANGULAR MATRIX B;       
        EXIT: THE UPPER TRIANGULAR MATRIX Q.B.Z (SEE BRIEF DESCRIPTION);      
     X: <ARRAY IDENTIFIER>; 
        "ARRAY" X[1:N,1:N]; 
        ENTRY: THE GIVEN MATRIX X;    
        EXIT: THE TRANSFORMED MATRIX Q.X.Z (SEE BRIEF DESCRIPTION); 
        
        
 PROCEDURES USED: 
     HSH2COL  = CP 34605    
     HSH2ROW3 = CP 34607    
        
        
 METHODE AND PERFORMANCE:   
     THE REDUCTION OF THE MATRIX A TO UPPER HESSENBERG FORM 
     WHILE PRESERVING THE TRIANGULARITY OF THE MATRIX B   
     IS THE RESULT OF A NUMBER OF STEPS, WHICH DO THE FOLLOWING     
     ACTIONS: INTRODUCING A ZERO ELEMENT IN A AND RESTORING 
     THE DISTURBED ZERO IN B, WITHOUT DISTURBING THE ZERO 
     INTRODUCED IN A. THIS IS DONE BY PRE-AND POSTMULTIPLICATIONS OF
     HOUSEHOLDER MATRICES.  
     THE MATRIX X SHARES THE TRANSFORMATION.    
     FOR FURTHER DETAILS SEE [1]      
        
        
 EXAMPLE OF USE:  
        
     THE PROCEDURE HESTGL3 IS USED IN QZI (THIS SECTION). 
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 11      
        
        
        
 SUBSECTION: HESTGL2:       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HESTGL2(N,A,B); "VALUE" N; "INTEGER" N; "ARRAY" A,B; 
     "CODE" 34604;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N: <ARITHMETIC EXPRESSION>;      
        THE ORDER OF THE GIVEN MATRICES;
     A: <ARRAY IDENTIFIER>; 
        "ARRAY" A[1:N,1:N]; 
        ENTRY: THE GIVEN MATRIX A;    
        EXIT: THE UPPER HESSENBERG MATRIX Q.A.Z (SEE BRIEF DESCRIPTION);      
     B: <ARRAY IDENTIFIER>; 
        "ARRAY" B[1:N,1:N]; 
        ENTRY: THE GIVEN UPPER TRIANGULAR MATRIX B
        EXIT: THE UPPER TRIANGULAR MATRIX Q.B.Z (SEE BRIEF DESCRIPTION);      
        
        
 PROCEDURES USED: 
     HSH2COL  = CP 34605    
     HSH2ROW2 = CP 34608    
        
        
 METHODE AND PERFORMANCE:   
     SEE HESTGL3, BUT HERE THE MATRIX X HAS BEEN LEFT OUT.
        
        
 EXAMPLE OF USE:  
     THE PROCEDURE HESTGL2 IS USED IN QZIVAL (THIS SECTION).
        
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 12      
        
        
        
 SUBSECTION HSH2COL:
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B); "VALUE" LA,LB,U,I,A1,A2;       
     "INTEGER" LA,LB,U,I; "REAL" A1,A2; "ARRAY" A,B;      
     "CODE" 34605;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     LA: <ARITHMETIC EXPRESSION>;     
         THE LOWER BOUND OF THE RUNNING COLUMN SUBSCRIPT OF A;      
     LB: <ARITHMETIC EXPRESSION>;     
         THE LOWER BOUND OF THE RUNNING COLUMN SUBSCRIPT OF B;      
     U:  <ARITHMETIC EXPRESSION>;     
         THE UPPER BOUND OF THE RUNNING COLUMN SUBSCRIPTS 
         OF A AND B.
     I:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPTS OF A AND B.   
         I+1 IS THE UPPERBOUND.       
     A1: <ARITHMETIC EXPRESSION>;     
         THE I-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED.
     A2: <ARITHMETIC EXPRESSION>;     
         THE (I+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[I:I+1,LA:U];       
         ENTRY THE GIVEN MATRIX A;    
         EXIT: THE TRANSFORMED MATRIX Q.A (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "ARRAY" B[I:I+1,LB:U];       
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX Q.B (SEE BRIEF DESCRIPTION);  
        
        
 PROCEDURES USED: 
     HSHVECMAT = CP 31070   
        
        
 METHOD AND PERFORMANCE:    
        WHEN THE CALCULATED HOUSEHOLDER MATRIX Q PREMULTIPLIES      
        A MATRIX M, ONLY ROWS I AND I+1 ARE CHANGED.      
        IF THE ELEMENTS I AND I+1 IN A COLUMN OF M ARE ZERO, THEY   
        REMAIN ZERO IN Q.M. 
        THEREFORE ONLY THE SUBMATRICES A[I:I+1,LA:U] AND  
        B[I:I+1,LB:U] ARE CHANGED, WHERE Q.A AND Q.B ARE  
        OVERWRITTEN IN A RESP B.      
        
        
 EXAMPLE OF USE: THE PROCEDURE HSH2COL IS USED IN QZI AND QZIVAL,   
                 (THIS SECTION).      
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 13      
        
        
        
 SUBSECTION HSH3COL:
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HSH3COL(LA,LB,U,I,A1,A2,A3,A,B); 
     "VALUE" LA,LB,U,I,A1,A2,A3;"INTEGER" LA,LB,I,U;"REAL" A1,A2,A3;
     "ARRAY" A,B; 
     "CODE" 34606;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     LA: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A;       
     LB: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF B;       
     U:  <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A AND B; 
     I:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF A AND B,    
         I+2 IS THE UPPERBOUND;       
     A1: <ARITHMETIC EXPRESSION>;     
         THE I-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;
     A2: <ARITHMETIC EXPRESSION>;     
         THE (I+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A3: <ARITHMETIC EXPRESSION>;     
         THE (I+2)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED.    
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[I:I+2,LA:U];       
         ENTRY: THE GIVEN MATRIX A;   
         EXIT: THE TRANSFORMED MATRIX Q.A (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "ARRAY" B[I:I+2,LB:U];       
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX Q.B (SEE BRIEF DESCRIPTION);  
        
        
 PROCEDURES USED: HSHVECMAT = CP 31070; 
        
        
 METHOD AND PERFORMANCE:    
        WHEN THE CALCULATED HOUSEHOLDER MATRIX Q PREMULTIPLIES A MATRIX       
        M, ONLY ROWS I, (I+1) AND (I+2) ARE CHANGED.      
        IF THE ELEMENTS I, I+1 AND I+2 ARE ZERO, THEN THEY REMAIN ZERO
        IN Q.M.   
        THEREFORE ONLY THE SUBMATRICES A[I:I+2,LA:U] AND B[I:I+2,LB:U]
        ARE CHANGED, WHERE Q.A AND Q.B ARE OVERWRITTEN IN A RESP B. 
        
        
 EXAMPLE OF USE: THE PROCEDURE HSH3COL IS USED IN QZI AND QZIVAL    
                 (THIS SECTION).      
        
1SECTION : 3.4.1.2            (FEBRUARY 1979)                    PAGE 14      
        
        
        
 SUBSECTION HSH2ROW3:       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HSH2ROW3(L,UA,UB,UX,J,A1,A2,A,B,X);      
     "VALUE" L,UA,UB,UX,J,A1,A2; "INTEGER" L,UA,UB,UX,J;  
     "REAL" A1,A2;"ARRAY" A,B,X;      
     "CODE" 34607;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     L:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF A,B AND X;  
     UA: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF A;
     UB: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF B;
     UX: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF X;
     J:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A,B AND X; 
         J+1 IS THE UPPERBOUND;       
     A1: <ARITHMETIC EXPRESSION>;     
         THE J-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;
     A2: <ARITHMETIC EXPRESSION>;     
         THE (J+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[L:UA,J:J+1];       
         ENTRY: THE GIVEN MATRIX A;   
         EXIT: THE TRANSFORMED MATRIX A.Z (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "ARRAY" B[L:UB,J:J+1];       
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX B.Z (SEE BRIEF DESCRIPTION);  
     X:  <ARRAY IDENTIFIER>;
         "ARRAY" X[L:UX,J:J+1];       
         ENTRY: THE GIVEN MATRIX X;   
         EXIT: THE TRANSFORMED MATRIX X.Z (SEE BRIEF DESCRIPTION);  
        
 PROCEDURES USED: HSHVECTAM = CP 31073; 
        
 METHOD AND PERFORMANCE:    
     WHEN THE CALCULATED HOUSEHOLDER MATRIX Z POSTMULTIPLIES
     A MATRIX M, ONLY COLUMNS J AND J+1 ARE CHANGED.      
     IF THE ELEMENTS J AND J+1 IN A ROW OF M ARE ZERO, THEN 
     THEY REMAIN ZERO IN M.Z
     THEREFORE ONLY THE SUBMATRICES A[L:UA,J:J+1],B[L:UB,J:J+1]     
     AND X[L:UX,J:J+1] ARE CHANGED, WHERE A.Z, B.Z AND X.Z ARE      
     OVERWRITTEN IN RESP. A,B AND X.  
        
 EXAMPLE OF USE: THE PROCEDURE HSH2ROW3 IS USED IN QZI (THIS SECTION).
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 15      
        
        
        
 SUBSECTION HSH2ROW2:       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE HSH2ROW2(LA,LB,UA,UB,J,A1,A2,A,B); "VALUE" LA,LB,UA,
     UB,J,A1,A2; "INTEGER" LA,LB,UA,UB,J; "REAL" A1,A2;"ARRAY" A,B; 
     "CODE" 34608;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     LA: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF A;
     LB: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF B;
     UA: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF A;
     UB: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF B;
     J:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A AND B; 
         J+1 IS THE UPPERBOUND;       
     A1: < ARITHMETIC EXPRESSION>;    
         THE J-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;
     A2: <ARITHMETIC EXPRESSION>;     
         THE (J+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" A[LA:UA,J:J+1]; 
         ENTRY: THE GIVEN MATRIX A;   
         EXIT: THE TRANSFORMED MATRIX A.Z (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" B[LB:UB,J:J+1]; 
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX B.Z (SEE BRIEF DESCRIPTION);  
        
        
 PROCEDURES USED: HSHVECTAM = CP 31073; 
        
        
 METHOD AND PERFORMANCE:    
     SEE HSH2ROW3, BUT HERE THE MATRIX X HAS BEEN LEFT OUT. 
        
        
 EXAMPLE OF USE:  
         THE PROCEDURE HSH2ROW2 IS USED IN QZIVAL,
         (THIS SECTION).    
        
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 16      
        
        
        
 SUBSECTION: HSH3ROW3.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HSH3ROW3(L,U,UX,J,A1,A2,A3,A,B,X);       
     "VALUE" L,U,UX,J,A1,A2,A3; "INTEGER" L,U,UX,J; "REAL" A1,A2,A3;
     "ARRAY" A,B,X; 
     "CODE" 34609;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     L:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF A AND B AND X;
     U:  <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF A AND B;    
     UX: <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF X;
     J:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A,B AND X; 
         J+2 IS THE UPPERBOUND;       
     A1: <ARITHMETIC EXPRESSION>;     
         THE J-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;
     A2: <ARITHMETIC EXPRESSION>;     
         THE (J+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A3: <ARITHMETIC EXPRESSION>;     
         THE (J+2)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" A[L:U,J:J+2]; 
         ENTRY: THE GIVEN MATRIX A;   
         EXIT: THE TRANSFORMED MATRIX A.Z (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" B[L:U,J:J+2]; 
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX B.Z (SEE BRIEF DESCRIPTION);  
     X:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" X[L:UX,J:J+2];
         ENTRY: THE GIVEN MATRIX X;   
         EXIT: THE TRANSFORMED MATRIX X.Z (SEE BRIEF DESCRIPTION);  
        
 PROCEDURES USED: HSHVECTAM = CP 31073; 
        
 METHOD AND PERFORMANCE:    
     WHEN THE CALCULATED HOUSEHOLDER MATRIX Z POSTMULTIPLIES A MATRIX M,      
     ONLY COLUMNS J,J+1 AND J+2 ARE CHANGED.    
     IF THE ELEMENTS J, J+1 AND J+2 IN A ROW OF M ARE ZERO, THEN THEY 
     REMAIN ZERO IN M.Z.    
     THEREFORE ONLY THE SUBMATRICES A[L:U,J:J+2], B[L:U,J:J+2] AND  
     X[L:UX,J:J+2] ARE CHANGED, WHERE A.Z , B.Z AND X.Z ARE OVERWRITTEN       
     ON RESP. A,B AND X;    
        
 EXAMPLE OF USE: THE PROCEDURE HSH3ROW3 IS USED IN QZI (THIS SECTION).
        
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 17      
        
        
        
 SUBSECTION: HSH3ROW2.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" HSH3ROW2(LA,LB,U,J,A1,A2,A3,A,B);
     "VALUE" LA,LB,U,J,A1,A2,A3; "INTEGER" LA,LB,U,J;     
     "REAL" A1,A2,A3; "ARRAY" A,B;    
     "CODE" 34610;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     LA: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF A;
     LB: <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING ROW SUBSCRIPT OF B;
     U:  <ARITHMETIC EXPRESSION>;     
         THE UPPERBOUND OF THE RUNNING ROW SUBSCRIPT OF A AND B;    
     J:  <ARITHMETIC EXPRESSION>;     
         THE LOWERBOUND OF THE RUNNING COLUMN SUBSCRIPT OF A AND B, 
         J+2 IS THE UPPERBOUND;       
     A1: <ARITHMETIC EXPRESSION>;     
         THE J-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;
     A2: <ARITHMETIC EXPRESSION>;     
         THE (J+1)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A3: <ARITHMETIC EXPRESSION>;     
         THE (J+2)-TH COMPONENT OF THE VECTOR TO BE TRANSFORMED;    
     A:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" A[LA:U,J:J+2];
         ENTRY: THE GIVEN MATRIX A;   
         EXIT: THE TRANSFORMED MATRIX A.Z (SEE BRIEF DESCRIPTION);  
     B:  <ARRAY IDENTIFIER>;
         "REAL" "ARRAY" B[LB:U,J:J+2];
         ENTRY: THE GIVEN MATRIX B;   
         EXIT: THE TRANSFORMED MATRIX B.Z (SEE BRIEF DESCRIPTION);  
        
        
 PROCEDURES USED: HSHVECTAM = CP 31073; 
        
        
 METHOD AND PERFORMANCE:    
        SEE HSH3ROW3, BUT HERE THE MATRIX X HAS BEEN LEFT OUT.      
        
        
 EXAMPLE OF USE: HSH3ROW2 IS USED IN QZIVAL (THIS SECTION). 
        
1SECTION : 3.4.1.2            (JANUARY 1979)                     PAGE 18      
        
        
        
 SOURCE TEXTS:    
0"CODE" 34600;    
 "PROCEDURE" QZIVAL(N,A,B,ALFR,ALFI,BETA,ITER,EM);
 "VALUE" N;"INTEGER" N;"ARRAY" A,B,ALFR,ALFI,BETA,EM;     
 "INTEGER" "ARRAY" ITER;    
 "BEGIN" "REAL" DWARF,EPS,EPSA,EPSB;  
 "PROCEDURE" ELMCOL(L,U,I,J,A,B,X);"CODE" 34023;
 "PROCEDURE" HSHDECMUL(N,A,B,DWARF);"CODE" 34602; 
 "PROCEDURE" HESTGL2(N,A,B);"CODE" 34604;       
 "PROCEDURE" HSH2ROW2(LA,LB,UA,UB,J,A1,A2,A,B);"CODE" 34608;
 "PROCEDURE" HSH3ROW2(LA,LB,U,J,A1,A2,A3,A,B);"CODE" 34610; 
 "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B);"CODE" 34605;   
 "PROCEDURE" HSH3COL(LA,LB,U,I,A1,A2,A3,A,B);"CODE" 34606;
 "PROCEDURE" CHSH2(A1R,A1I,A2R,A2I,C,SR,SI);"CODE" 34611; 
 "PROCEDURE" HSHVECMAT(LR,UR,LC,UC,X,U,A);"CODE" 31070;   
 "PROCEDURE" HSHVECTAM(LR,UR,LC,UC,X,U,A);"CODE" 31073;   
 "PROCEDURE" QZIT(N,A,B,EPS,EPSA,EPSB,ITER);"VALUE"N,EPS; 
 "REAL" EPS,EPSA,EPSB;"INTEGER" N;"INTEGER" "ARRAY" ITER;"ARRAY" A,B; 
 "BEGIN" "REAL" ANORM,BNORM,ANI,BNI,CONST,A10,A20,A30,B11,
     B22,B33,B44,A11,A12,A21,A22,A33,A34,A43,A44,B12,B34,OLD1,OLD2; 
     "INTEGER" I,Q,M,M1,Q1,J,K,K1,K2,K3,KM1;"BOOLEAN" STATIONARY;   
     ANORM:=BNORM:=0;"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"   
     "BEGIN" BNI:=0;ITER[I]:=0;ANI:="IF" I>1"THEN"ABS(A[I,I-1])"ELSE" 0;      
         "FOR" J:=I "STEP" 1 "UNTIL" N "DO"     
         "BEGIN" ANI:=ANI+ABS(A[I,J]);BNI:=BNI+ABS(B[I,J])
         "END";"IF" ANI>ANORM "THEN" ANORM:=ANI;"IF" BNI>BNORM"THEN"
         BNORM:=BNI 
     "END";"IF" ANORM=0 "THEN" ANORM:=EPS;"IF"BNORM=0 "THEN" BNORM:=EPS;      
     EPSA:=EPS*ANORM;EPSB:=EPS*BNORM; 
     "FOR" M:=N,M "WHILE" M>=3 "DO"   
   "BEGIN"
     "FOR" I:=M+1,I-1 "WHILE"("IF" I>1 "THEN"ABS(A[I,I-1])>EPSA "ELSE"
     "FALSE") "DO" Q:=I-1;  
     "IF" Q>1 "THEN" A[Q,Q-1]:=0;     
  L: "IF" Q>=M-1 "THEN" M:=Q-1 "ELSE" 
     "BEGIN"      
         "IF" ABS(B[Q,Q])<=EPSB "THEN"
         "BEGIN" B[Q,Q]:=0;Q1:=Q+1;   
             HSH2COL(Q,Q,M,Q,A[Q,Q],A[Q1,Q],A,B);A[Q1,Q]:=0;
             Q:=Q1;"GOTO" L 
         "END" "ELSE" M1:=M-1;Q1:=Q+1;CONST:=0.75;ITER[M]:=ITER[M]+1; 
         STATIONARY:="IF" ITER[M]=1 "THEN" "TRUE" "ELSE"  
         ABS(A[M,M-1])>=CONST*OLD1"AND"ABS(A[M-1,M-2])>=CONST*OLD2; 
         "IF" ITER[M]>30"AND"STATIONARY "THEN"  
         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" ITER[I]:=-1;    
             "GOTO" OUT     
         "END";   
         "IF" ITER[M]=10"AND"STATIONARY "THEN"  
         "BEGIN" A10:=0;A20:=1;A30:=1.1605      
         "END" "ELSE"       
         "BEGIN" B11:=B[Q,Q];B22:="IF" ABS(B[Q1,Q1])<EPSB "THEN"EPSB
             "ELSE" B[Q1,Q1]; 
             B33:="IF" ABS(B[M1,M1])<EPSB "THEN" EPSB "ELSE"B[M1,M1]; 
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 19      
                                                                 ;  
        
        
             B44:="IF" ABS(B[M,M])<EPSB "THEN" EPSB "ELSE" B[M,M] ; 
             A11:=A[Q,Q]/B11;A12:=A[Q,Q1]/B22;A21:=A[Q1,Q]/B11;     
             A22:=A[Q1,Q1]/B22;A33:=A[M1,M1]/B33;A34:=A[M1,M]/B44;  
             A43:=A[M,M1]/B33;A44:=A[M,M]/B44;B12:=B[Q,Q1]/B22;     
             B34:=B[M1,M]/B44;
             A10:=((A33-A11)*(A44-A11)-A34*A43+A43*B34*A11)/A21     
                  +A12-A11*B12;       
             A20:=(A22-A11-A21*B12)-(A33-A11)-(A44-A11)+A43*B34;    
             A30:=A[Q+2,Q1]/B22       
         "END";OLD1:=ABS(A[M,M-1]);OLD2:=ABS(A[M-1,M-2]); 
         "FOR" K:=Q "STEP" 1 "UNTIL" M1 "DO"    
         "BEGIN" K1:=K+1;K2:=K+2;K3:="IF" K+3>M "THEN" M "ELSE" K+3;
             KM1:="IF" K-1<Q "THEN" Q "ELSE" K-1; 
             "IF" K^=M1 "THEN"
             "BEGIN" "IF" K=Q "THEN"  
                 "BEGIN"    
                  HSH3COL(KM1,KM1,M,K,A[K,KM1],A[K1,KM1],A[K2,KM1],A,B);      
                   A[K1,KM1]:=A[K2,KM1]:=0      
                 "END";     
                 HSH3ROW2(Q,Q,K3,K,B[K2,K2],B[K2,K1],B[K2,K],A,B);  
                 B[K2,K]:=B[K2,K1]:=0 ; 
             "END" "ELSE"   
             "BEGIN" HSH2COL(KM1,KM1,M,K,A[K,KM1],A[K1,KM1],A,B);   
                 A[K1,KM1]:=0 
             "END"; 
             HSH2ROW2(Q,Q,K3,K3,K,B[K1,K1],B[K1,K],A,B);B[K1,K]:=0  
         "END"    
     "END";       
     OUT: 
   "END"
 "END" QZIT;      
         "EOP"    
 "PROCEDURE" QZVAL(N,A,B,EPSA,EPSB,ALFR,ALFI,BETA);"VALUE" N;       
 "REAL" EPSA,EPSB;"INTEGER" N;"ARRAY" ALFR,ALFI,BETA,A,B; 
 "BEGIN" "INTEGER" M,L,J;"REAL" AN,BN,A11,A12,A21,A22,B11,B12,B22,E,C,D,      
     ER,EI,A11R,A11I,A12R,A12I,A21R,A21I,A22R,A22I,CZ,SZR,SZI,      
     CQ,SQR,SQI,SSR,SSI,TR,TI,BDR,BDI,R;
     "FOR" M:=N,M "WHILE" M>0 "DO"    
     "IF"("IF" M>1 "THEN" A[M,M-1]=0 "ELSE" "TRUE") "THEN"
     "BEGIN" ALFR[M]:=A[M,M];BETA[M]:=B[M,M];ALFI[M]:=0;M:=M-1      
     "END" "ELSE" 
     "BEGIN" L:=M-1;"IF" ABS(B[L,L])<=EPSB "THEN" 
         "BEGIN" B[L,L]:=0;HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);      
             A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];     
             BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;   
         "END" "ELSE" "IF" ABS(B[M,M])<=EPSB "THEN"       
         "BEGIN"B[M,M]:=0;HSH2ROW2(1,1,M,M,L,A[M,M],A[M,L],A,B);    
             A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];     
             BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;   
         "END""ELSE"
       "BEGIN"    
         AN:=ABS(A[L,L])+ABS(A[L,M])+ABS(A[M,L])+ABS(A[M,M]);       
         BN:=ABS(B[L,L])+ABS(B[L,M])+ABS(B[M,M]); 
         A11:=A[L,L]/AN;A12:=A[L,M]/AN;A21:=A[M,L]/AN;A22:=A[M,M]/AN; 
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 20      
                                                                 ;  
        
        
         B11:=B[L,L]/BN;B12:=B[L,M]/BN;B22:=B[M,M]/BN;    
         E:=A11/B11;
         C:=((A22-E*B22)/B22-(A21*B12)/(B11*B22))/2;      
         D:=C*C+(A21*(A12-E*B12))/(B11*B22);    
         "IF" D>=0 "THEN"   
         "BEGIN"E:=E+("IF"C<0"THEN"C-SQRT(D)"ELSE"C+SQRT(D));       
             A11:=A11-E*B11;A12:=A12-E*B12;A22:=A22-E*B22;
             "IF" ABS(A11)+ABS(A12)>=ABS(A21)+ABS(A22) "THEN"       
             HSH2ROW2(1,1,M,M,L,A12,A11,A,B)"ELSE"
             HSH2ROW2(1,1,M,M,L,A22,A21,A,B);   
             "IF"AN>=ABS(E)*BN"THEN"  
             HSH2COL(L,L,N,L,B[L,L],B[M,L],A,B) "ELSE"    
             HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);
             A[M,L]:=B[M,L]:=0;       
             ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];BETA[L]:=B[L,L];       
             BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0; 
         "END""ELSE"
         "BEGIN"  
             ER:=E+C;EI:=SQRT(-D);A11R:=A11-ER*B11;A11I:=EI*B11;    
             A12R:=A12-ER*B12;A12I:=EI*B12;A21R:=A21;A21I:=0;       
             A22R:=A22-ER*B22;A22I:=EI*B22;     
             "IF"ABS(A11R)+ABS(A11I)+ABS(A12R)+ABS(A12I)>=
             ABS(A21R)+ABS(A22R)+ABS(A22I)"THEN"
             CHSH2(A12R,A12I,-A11R,-A11I,CZ,SZR,SZI)"ELSE"
             CHSH2(A22R,A22I,-A21R,-A21I,CZ,SZR,SZI);     
             "IF"AN>=(ABS(ER)+ABS(EI))*BN"THEN" 
             CHSH2(CZ*B11+SZR*B12,SZI*B12,SZR*B22,SZI*B22,CQ,SQR,SQI) 
             "ELSE"CHSH2(CZ*A11+SZR*A12,SZI*A12,CZ*A21+SZR*A22,SZI*A22,       
             CQ,SQR,SQI);SSR:=SQR*SZR+SQI*SZI;SSI:=SQR*SZI-SQI*SZR; 
             TR:=CQ*CZ*A11+CQ*SZR*A12+SQR*CZ*A21+SSR*A22; 
             TI:=CQ*SZI*A12-SQI*CZ*A21+SSI*A22; 
             BDR:=CQ*CZ*B11+CQ*SZR*B12+SSR*B22; 
             BDI:=CQ*SZI*B12+SSI*B22; 
             R:=SQRT(BDR*BDR+BDI*BDI);BETA[L]:=BN*R;      
             ALFR[L]:=AN*(TR*BDR+TI*BDI)/R;     
             ALFI[L]:=AN*(TR*BDI-TI*BDR)/R;     
             TR:=SSR*A11-SQR*CZ*A12-CQ*SZR*A21+CQ*CZ*A22; 
             TI:=-SSI*A11-SQI*CZ*A12+CQ*SZI*A21;
             BDR:=SSR*B11-SQR*CZ*B12+CQ*CZ*B22; 
             BDI:=-SSI*B11-SQI*CZ*B12;
             R:=SQRT(BDR*BDR+BDI*BDI);BETA[M]:=BN*R;      
             ALFR[M]:=AN*(TR*BDR+TI*BDI)/R;     
             ALFI[M]:=AN*(TR*BDI-TI*BDR)/R;     
         "END"    
       "END";M:=M-2 
     "END"
 "END" QZVAL;     
 DWARF:=EM[0];EPS:=EM[1];   
 HSHDECMUL(N,A,B,DWARF);    
 HESTGL2(N,A,B);  
 QZIT(N,A,B,EPS,EPSA,EPSB,ITER);      
 QZVAL(N,A,B,EPSA,EPSB,ALFR,ALFI,BETA); 
 "END" QZIVAL;    
          "EOP"   
        
        
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 21      
        
        
        
0"CODE" 34601;    
 "PROCEDURE" QZI(N,A,B,X,ALFR,ALFI,BETA,ITER,EM); 
 "VALUE" N;"INTEGER" N;"ARRAY"A,B,X,ALFR,ALFI,BETA,EM;    
 "INTEGER" "ARRAY" ITER;    
 "BEGIN" "REAL" DWARF,EPS,EPSA,EPSB;  
 "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;     
 "PROCEDURE" HSHDECMUL(N,A,B,DWARF);"CODE" 34602; 
 "PROCEDURE" HESTGL3(N,A,B,X);"CODE" 34603;     
 "PROCEDURE" HSH2ROW3(L,UA,UB,UX,J,A1,A2,A,B,X);"CODE" 34607;       
 "PROCEDURE" HSH3ROW3(L,U,UX,J,A1,A2,A3,A,B,X);"CODE" 34609;
 "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B);"CODE" 34605;   
 "PROCEDURE" HSH3COL(LA,LB,U,I,A1,A2,A3,A,B);"CODE" 34606;
 "PROCEDURE" CHSH2(A1R,A1I,A2R,A2I,C,SR,SI);"CODE" 34611; 
 "PROCEDURE" COMDIV(XR,XI,YR,YI,ZR,ZI);"CODE" 34342;      
 "PROCEDURE" QZIT(N,A,B,X,EPS,EPSA,EPSB,ITER);"VALUE"N,EPS; 
 "REAL" EPS,EPSA,EPSB;"INTEGER" N;"INTEGER" "ARRAY" ITER;"ARRAY" A,B,X;       
 "BEGIN" "REAL" ANORM,BNORM,ANI,BNI,CONST,A10,A20,A30,B11,
     B22,B33,B44,A11,A12,A21,A22,A33,A34,A43,A44,B12,B34,OLD1,OLD2; 
     "INTEGER" I,Q,M,M1,Q1,J,K,K1,K2,K3,KM1;"BOOLEAN" STATIONARY;   
     ANORM:=BNORM:=0;"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"   
     "BEGIN" BNI:=0;ITER[I]:=0;ANI:="IF" I>1"THEN"ABS(A[I,I-1])"ELSE" 0;      
         "FOR" J:=I "STEP" 1 "UNTIL" N "DO"     
         "BEGIN" ANI:=ANI+ABS(A[I,J]);BNI:=BNI+ABS(B[I,J])
         "END";"IF" ANI>ANORM "THEN" ANORM:=ANI;"IF" BNI>BNORM"THEN"
         BNORM:=BNI 
     "END";"IF" ANORM=0 "THEN" ANORM:=EPS;"IF"BNORM=0 "THEN" BNORM:=EPS;      
     EPSA:=EPS*ANORM;EPSB:=EPS*BNORM; 
     "FOR" M:=N,M "WHILE" M>=3 "DO"   
   "BEGIN"
     "FOR" I:=M+1,I-1 "WHILE"("IF" I>1 "THEN"ABS(A[I,I-1])>EPSA "ELSE"
     "FALSE") "DO" Q:=I-1;  
     "IF" Q>1 "THEN" A[Q,Q-1]:=0;     
  L: "IF" Q>=M-1 "THEN" M:=Q-1 "ELSE" 
     "BEGIN"      
         "IF" ABS(B[Q,Q])<=EPSB "THEN"
         "BEGIN" B[Q,Q]:=0;Q1:=Q+1;   
             HSH2COL(Q,Q,N,Q,A[Q,Q],A[Q1,Q],A,B);A[Q1,Q]:=0;
             Q:=Q1;"GOTO" L 
         "END" "ELSE" M1:=M-1;Q1:=Q+1;CONST:=0.75;ITER[M]:=ITER[M]+1; 
         STATIONARY:="IF" ITER[M]=1 "THEN" "TRUE" "ELSE"  
         ABS(A[M,M-1])>=CONST*OLD1"AND"ABS(A[M-1,M-2])>=CONST*OLD2; 
         "IF" ITER[M]>30"AND"STATIONARY "THEN"  
         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" ITER[I]:=-1;    
             "GOTO" OUT     
         "END";   
         "IF" ITER[M]=10"AND"STATIONARY "THEN"  
         "BEGIN" A10:=0;A20:=1;A30:=1.1605      
         "END" "ELSE"       
         "BEGIN" B11:=B[Q,Q];B22:="IF" ABS(B[Q1,Q1])<EPSB "THEN"EPSB
             "ELSE" B[Q1,Q1]; 
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 22      
                                                                 ;  
        
        
             B33:="IF" ABS(B[M1,M1])<EPSB "THEN" EPSB "ELSE"B[M1,M1]; 
             B44:="IF" ABS(B[M,M])<EPSB "THEN" EPSB "ELSE" B[M,M] ; 
             A11:=A[Q,Q]/B11;A12:=A[Q,Q1]/B22;A21:=A[Q1,Q]/B11;     
             A22:=A[Q1,Q1]/B22;A33:=A[M1,M1]/B33;A34:=A[M1,M]/B44;  
             A43:=A[M,M1]/B33;A44:=A[M,M]/B44;B12:=B[Q,Q1]/B22;     
             B34:=B[M1,M]/B44;
             A10:=((A33-A11)*(A44-A11)-A34*A43+A43*B34*A11)/A21     
                  +A12-A11*B12;       
             A20:=(A22-A11-A21*B12)-(A33-A11)-(A44-A11)+A43*B34;    
             A30:=A[Q+2,Q1]/B22       
         "END";OLD1:=ABS(A[M,M-1]);OLD2:=ABS(A[M-1,M-2]); 
         "FOR" K:=Q "STEP" 1 "UNTIL" M1 "DO"    
         "BEGIN" K1:=K+1;K2:=K+2;K3:="IF" K+3>M "THEN" M "ELSE" K+3;
             KM1:="IF" K-1<Q "THEN" Q "ELSE" K-1; 
             "IF" K^=M1 "THEN"
             "BEGIN" "IF" K=Q "THEN"  
                 HSH3COL(KM1,KM1,N,K,A10,A20,A30,A,B) "ELSE"
                 "BEGIN"    
                  HSH3COL(KM1,KM1,N,K,A[K,KM1],A[K1,KM1],A[K2,KM1],A,B);      
                   A[K1,KM1]:=A[K2,KM1]:=0      
                 "END";     
                 HSH3ROW3(1,K3,N,K,B[K2,K2],B[K2,K1],B[K2,K],A,B,X);
                 B[K2,K]:=B[K2,K1]:=0 ; 
             "END" "ELSE"   
             "BEGIN" HSH2COL(KM1,KM1,N,K,A[K,KM1],A[K1,KM1],A,B);   
                 A[K1,KM1]:=0 
             "END"; 
             HSH2ROW3(1,K3,K3,N,K,B[K1,K1],B[K1,K],A,B,X);B[K1,K]:=0
         "END"    
     "END"
   "END"; OUT:    
 "END" QZIT;      
 "PROCEDURE" QZVAL(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);"VALUE" N;     
 "REAL" EPSA,EPSB;"INTEGER" N;"ARRAY" ALFR,ALFI,BETA,A,B,X; 
 "BEGIN" "INTEGER" M,L,J;"REAL" AN,BN,A11,A12,A21,A22,B11,B12,B22,E,C,D,      
     ER,EI,A11R,A11I,A12R,A12I,A21R,A21I,A22R,A22I,CZ,SZR,SZI,      
     CQ,SQR,SQI,SSR,SSI,TR,TI,BDR,BDI,R;
     "FOR" M:=N,M "WHILE" M>0 "DO"    
     "IF"("IF" M>1 "THEN" A[M,M-1]=0 "ELSE" "TRUE")"THEN" 
     "BEGIN" ALFR[M]:=A[M,M];BETA[M]:=B[M,M];ALFI[M]:=0;M:=M-1      
     "END" "ELSE" 
     "BEGIN" L:=M-1;"IF" ABS(B[L,L])<=EPSB "THEN" 
         "BEGIN" B[L,L]:=0;HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);      
             A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];     
             BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;   
         "END" "ELSE" "IF" ABS(B[M,M])<=EPSB "THEN"       
         "BEGIN"B[M,M]:=0;HSH2ROW3(1,M,M,N,L,A[M,M],A[M,L],A,B,X);  
             A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];     
             BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;   
         "END""ELSE"
       "BEGIN"    
         AN:=ABS(A[L,L])+ABS(A[L,M])+ABS(A[M,L])+ABS(A[M,M]);       
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (FEBRUARY 1979)                    PAGE 23      
                                                                 ;  
        
        
         BN:=ABS(B[L,L])+ABS(B[L,M])+ABS(B[M,M]); 
         A11:=A[L,L]/AN;A12:=A[L,M]/AN;A21:=A[M,L]/AN;A22:=A[M,M]/AN; 
         B11:=B[L,L]/BN;B12:=B[L,M]/BN;B22:=B[M,M]/BN;    
         E:=A11/B11;
         C:=((A22-E*B22)/B22-(A21*B12)/(B11*B22))/2;      
         D:=C*C+(A21*(A12-E*B12))/(B11*B22);    
         "IF" D>=0 "THEN"   
         "BEGIN"E:=E+("IF"C<0"THEN"C-SQRT(D)"ELSE"C+SQRT(D));       
             A11:=A11-E*B11;A12:=A12-E*B12;A22:=A22-E*B22;
             "IF" ABS(A11)+ABS(A12)>=ABS(A21)+ABS(A22) "THEN"       
             HSH2ROW3(1,M,M,N,L,A12,A11,A,B,X)"ELSE"      
             HSH2ROW3(1,M,M,N,L,A22,A21,A,B,X); 
             "IF"AN>=ABS(E)*BN"THEN"  
             HSH2COL(L,L,N,L,B[L,L],B[M,L],A,B) "ELSE"    
             HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);
             A[M,L]:=B[M,L]:=0;       
             ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];BETA[L]:=B[L,L];       
             BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0; 
         "END""ELSE"
         "BEGIN"  
             ER:=E+C;EI:=SQRT(-D);A11R:=A11-ER*B11;A11I:=EI*B11;    
             A12R:=A12-ER*B12;A12I:=EI*B12;A21R:=A21;A21I:=0;       
             A22R:=A22-ER*B22;A22I:=EI*B22;     
             "IF"ABS(A11R)+ABS(A11I)+ABS(A12R)+ABS(A12I)>=
             ABS(A21R)+ABS(A22R)+ABS(A22I)"THEN"
             CHSH2(A12R,A12I,-A11R,-A11I,CZ,SZR,SZI)"ELSE"
             CHSH2(A22R,A22I,-A21R,-A21I,CZ,SZR,SZI);     
             "IF"AN>=(ABS(ER)+ABS(EI))*BN"THEN" 
             CHSH2(CZ*B11+SZR*B12,SZI*B12,SZR*B22,SZI*B22,CQ,SQR,SQI) 
             "ELSE"CHSH2(CZ*A11+SZR*A12,SZI*A12,CZ*A21+SZR*A22,SZI*A22,       
             CQ,SQR,SQI);SSR:=SQR*SZR+SQI*SZI;SSI:=SQR*SZI-SQI*SZR; 
             TR:=CQ*CZ*A11+CQ*SZR*A12+SQR*CZ*A21+SSR*A22; 
             TI:=CQ*SZI*A12-SQI*CZ*A21+SSI*A22; 
             BDR:=CQ*CZ*B11+CQ*SZR*B12+SSR*B22; 
             BDI:=CQ*SZI*B12+SSI*B22; 
             R:=SQRT(BDR*BDR+BDI*BDI);BETA[L]:=BN*R;      
             ALFR[L]:=AN*(TR*BDR+TI*BDI)/R;     
             ALFI[L]:=AN*(TR*BDI-TI*BDR)/R;     
             TR:=SSR*A11-SQR*CZ*A12-CQ*SZR*A21+CQ*CZ*A22; 
             TI:=-SSI*A11-SQI*CZ*A12+CQ*SZI*A21;
             BDR:=SSR*B11-SQR*CZ*B12+CQ*CZ*B22; 
             BDI:=-SSI*B11-SQI*CZ*B12;
             R:=SQRT(BDR*BDR+BDI*BDI);BETA[M]:=BN*R;      
             ALFR[M]:=AN*(TR*BDR+TI*BDI)/R;     
             ALFI[M]:=AN*(TR*BDI-TI*BDR)/R;     
         "END"    
       "END";M:=M-2 
     "END"
 "END" QZVAL      
         "EOP"    
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 24      
                                                                 ;  
        
        
 "PROCEDURE" QZVEC(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);"VALUE"N,EPSA,EPSB;      
 "INTEGER" N;"REAL" EPSA,EPSB;"ARRAY" A,B,ALFR,ALFI,BETA,X; 
 "BEGIN""INTEGER" M,MR,MI,L,L1,J,K;"REAL" BETM,ALFM,SL,SK,D,TKK,TKL,TLK,      
     TLL,ALMI,ALMR,TR,TI,SLR,SLI,SKR,SKI,DR,DI,TKKR,TKKI,TKLR,TKLI,TLKR,      
     TLKI,TLLR,TLLI,S,R;    
     "FOR" M:=N "STEP" -1 "UNTIL" 1 "DO"
     "IF" ALFI[M]=0 "THEN"  
     "BEGIN" "COMMENT" M-TH REAL VECTOR;
         ALFM:=ALFR[M];BETM:=BETA[M];B[M,M]:=1;L1:=M;     
         "FOR" L:=M-1 "STEP" -1 "UNTIL" 1 "DO"  
         "BEGIN" SL:=0;     
             "FOR" J:=L1 "STEP" 1 "UNTIL" M "DO"
             SL:=SL+(BETM*A[L,J]-ALFM*B[L,J])*B[J,M];     
             "IF"("IF"L^=1"THEN"BETM*A[L,L-1]=0"ELSE""TRUE")"THEN"  
             "BEGIN""COMMENT" 1-1 BLOCK;
                 D:=BETM*A[L,L]-ALFM*B[L,L];    
                 "IF" D=0 "THEN" D:=(EPSA+EPSB)/2;B[L,M]:=-SL/D     
             "END""ELSE"    
             "BEGIN" "COMMENT" 2-2 BLOCK;K:=L-1;SK:=0;    
                 "FOR"J:=L1 "STEP" 1 "UNTIL" M "DO"       
                 SK:=SK+(BETM*A[K,J]-ALFM*B[K,J])*B[J,M]; 
                 TKK:=BETM*A[K,K]-ALFM*B[K,K];  
                 TKL:=BETM*A[K,L]-ALFM*B[K,L];  
                 TLK:=BETM*A[L,K];    
                 TLL:=BETM*A[L,L]-ALFM*B[L,L];  
                 D:=TKK*TLL-TKL*TLK;"IF" D=0 "THEN" D:=(EPSA+EPSB)/2; 
                 B[L,M]:=(TLK*SK-TKK*SL)/D;     
                B[K,M]:="IF"ABS(TKK)>=ABS(TLK)"THEN"-(SK+TKL*B[L,M])/TKK      
                 "ELSE" -(SL+TLL*B[L,M])/TLK;L:=L-1       
             "END";L1:=L    
         "END"    
     "END""ELSE"  
     "BEGIN" "COMMENT" COMPLEX VECTOR;
         ALMR:=ALFR[M-1];ALMI:=ALFI[M-1];BETM:=BETA[M-1];MR:=M-1;MI:=M;       
         B[M-1,MR]:=ALMI*B[M,M]/(BETM*A[M,M-1]);
         B[M-1,MI]:=(BETM*A[M,M]-ALMR*B[M,M])/(BETM*A[M,M-1]);      
         B[M,MR]:=0;B[M,MI]:=-1;L1:=M-1;
         "FOR" L:=M-2 "STEP" -1 "UNTIL" 1 "DO"  
         "BEGIN" SLR:=SLI:=0; 
             "FOR" J:=L1 "STEP" 1 "UNTIL" M "DO"
             "BEGIN" TR:=BETM*A[L,J]-ALMR*B[L,J]; 
                 TI:=-ALMI*B[L,J];    
                 SLR:=SLR+TR*B[J,MR]-TI*B[J,MI];
                 SLI:=SLI+TR*B[J,MI]+TI*B[J,MR] 
             "END"; 
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 25      
                                                                 ;  
        
        
             "IF"("IF"L^=1"THEN"BETM*A[L,L-1]=0"ELSE""TRUE")"THEN"  
             "BEGIN"DR:=BETM*A[L,L]-ALMR*B[L,L];
                 DI:=-ALMI*B[L,L];    
                 COMDIV(-SLR,-SLI,DR,DI,B[L,MR],B[L,MI]); 
             "END""ELSE"    
             "BEGIN" K:=L-1;SKR:=SKI:=0;
                 "FOR" J:=L1 "STEP" 1 "UNTIL" M "DO"      
                 "BEGIN" TR:=BETM*A[K,J]-ALMR*B[K,J];     
                     TI:=-ALMI*B[K,J];
                     SKR:=SKR+TR*B[J,MR]-TI*B[J,MI];      
                     SKI:=SKI+TR*B[J,MI]+TI*B[J,MR]       
                 "END";     
                 TKKR:=BETM*A[K,K]-ALMR*B[K,K]; 
                 TKKI:=-ALMI*B[K,K];  
                 TKLR:=BETM*A[K,L]-ALMR*B[K,L]; 
                 TKLI:=-ALMI*B[K,L];  
                 TLKR:=BETM*A[L,K];TLKI:=0;     
                 TLLR:=BETM*A[L,L]-ALMR*B[L,L]; 
                 TLLI:=-ALMI*B[L,L];  
                 DR:=TKKR*TLLR-TKKI*TLLI-TKLR*TLKR;       
                 DI:=TKKR*TLLI+TKKI*TLLR-TKLI*TLKR;       
                 "IF"DR=0"AND"DI=0"THEN"DR:=(EPSA+EPSB)/2;
                 COMDIV(TLKR*SKR-TKKR*SLR+TKKI*SLI,TLKR*SKI-TKKR*SLI- 
                 TKKI*SLR,DR,DI,B[L,MR],B[L,MI]); 
                 "IF"ABS(TKKR)+ABS(TKKI)>=ABS(TLKR)"THEN" 
                 COMDIV(-SKR-TKLR*B[L,MR]+TKLI*B[L,MI],-SKI-TKLR*B[L,MI]      
                 -TKLI*B[L,MR],TKKR,TKKI,B[K,MR],B[K,MI])"ELSE"     
                 COMDIV(-SLR-TLLR*B[L,MR]+TLLI*B[L,MI],-SLI-TLLR*B[L,MI]      
                 -TLLI*B[L,MR],TLKR,TLKI,B[K,MR],B[K,MI]);L:=L-1    
             "END";L1:=L    
         "END";M:=M-1       
     "END";       
     "FOR" M:=N "STEP" -1 "UNTIL" 1 "DO"
     "FOR" K:=1 "STEP" 1 "UNTIL" N "DO" 
     X[K,M]:=MATMAT(1,M,K,M,X,B);     
     "FOR" M:=N "STEP" -1 "UNTIL" 1 "DO"
     "BEGIN" S:=0;"IF" ALFI[M]=0 "THEN" 
         "BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"       
             "BEGIN" R:=ABS(X[K,M]);  
                 "IF" R>=S "THEN""BEGIN" S:=R;D:=X[K,M] "END"       
             "END";"FOR" K:=1 "STEP" 1 "UNTIL" N "DO"     
             X[K,M]:=X[K,M]/D 
         "END""ELSE"
         "BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"       
             "BEGIN" R:=ABS(X[K,M-1])+ABS(X[K,M]);
                 R:=R*SQRT((X[K,M-1]/R)**2+(X[K,M]/R)**2);
                 "IF" R>=S"THEN"      
                 "BEGIN" S:=R;DR:=X[K,M-1];DI:=X[K,M] "END" 
             "END";"FOR" K:=1 "STEP" 1 "UNTIL" N "DO"     
             COMDIV(X[K,M-1],X[K,M],DR,DI,X[K,M-1],X[K,M]);M:=M-1   
         "END"    
     "END"
 "END" QZVEC;     
                                                               "COMMENT"      
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 26      
                                                                 ;  
        
        
 DWARF:=EM[0];EPS:=EM[1];   
 HSHDECMUL(N,A,B,DWARF);    
 HESTGL3(N,A,B,X);
 QZIT(N,A,B,X,EPS,EPSA,EPSB,ITER);    
 QZVAL(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);       
 QZVEC(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA)
 "END" QZI;       
          "EOP"   
        
        
        
0"CODE" 34602;    
 "PROCEDURE" HSHDECMUL(N,A,B,DWARF);"VALUE"N,DWARF;"INTEGER"N;      
 "REAL" DWARF;"ARRAY" A,B;  
 "BEGIN"  "ARRAY" V[1:N];"INTEGER" J,K,K1,N1;"REAL" R,T,C;
     "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014; 
     "PROCEDURE" HSHVECMAT(LR,UR,LC,UC,X,U,A);"CODE" 31070; 
     K:=1;N1:=N+1;
     "FOR" K1:=2 "STEP" 1 "UNTIL" N1 "DO"       
     "BEGIN" R:=TAMMAT(K1,N,K,K,B,B); 
         "IF" R>DWARF "THEN"
         "BEGIN" R:="IF" B[K,K]<0 "THEN" -SQRT(R+B[K,K]*B[K,K])     
             "ELSE" SQRT(R+B[K,K]*B[K,K]);T:=B[K,K]+R;C:=-T/R;      
             B[K,K]:=-R;V[K]:=1;      
             "FOR" J:=K1 "STEP" 1 "UNTIL" N "DO" V[J]:=B[J,K]/T;    
             HSHVECMAT(K,N,K1,N,C,V,B);HSHVECMAT(K,N,1,N,C,V,A)     
         "END";K:=K1
     "END"
 "END" HSHDECMUL; 
          "EOP"   
        
0"CODE" 34603;    
 "PROCEDURE" HESTGL3(N,A,B,X);"VALUE" N;"INTEGER" N;"ARRAY" A,B,X;  
 "BEGIN" "INTEGER" NM1,K,L,K1,L1;     
     "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B);"CODE" 34605; 
     "PROCEDURE" HSH2ROW3(L,UA,UB,UX,J,A1,A2,A,B,X);"CODE" 34607;   
     "IF" N>2 "THEN"
     "BEGIN" "FOR" K:=2 "STEP" 1 "UNTIL" N "DO" 
         "FOR" L:=1 "STEP" 1 "UNTIL" K-1 "DO" B[K,L]:=0;  
         NM1:=N-1;K:=1;     
         "FOR" K1:= 2 "STEP" 1 "UNTIL" NM1 "DO" 
         "BEGIN" L1:=N;     
             "FOR" L:=N-1 "STEP" -1 "UNTIL" K1 "DO"       
             "BEGIN"
                 HSH2COL(K,L,N,L,A[L,K],A[L1,K],A,B);A[L1,K]:=0;    
                 HSH2ROW3(1,N,L1,N,L,B[L1,L1],B[L1,L],A,B,X);       
                 B[L1,L]:=0;L1:=L     
             "END";K:=K1    
         "END"    
     "END"
 "END" HESTGL3;   
          "EOP"   
        
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 27      
        
        
        
0"CODE" 34604;    
 "PROCEDURE" HESTGL2(N,A,B);"VALUE" N;"INTEGER" N;"ARRAY" A,B;      
 "BEGIN" "INTEGER" NM1,K,L,K1,L1;     
     "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B);"CODE" 34605; 
     "PROCEDURE" HSH2ROW2(LA,LB,UA,UB,A1,A2,A,B);"CODE" 34608;      
     "IF" N>2 "THEN"
     "BEGIN" "FOR" K:=2 "STEP" 1 "UNTIL" N "DO" 
         "FOR" L:=1 "STEP" 1 "UNTIL" K-1 "DO" B[K,L]:=0;  
         NM1:=N-1;K:=1;     
         "FOR" K1:= 2 "STEP" 1 "UNTIL" NM1 "DO" 
         "BEGIN" L1:=N;     
             "FOR" L:=N-1 "STEP" -1 "UNTIL" K1 "DO"       
             "BEGIN"
                 HSH2COL(K,L,N,L,A[L,K],A[L1,K],A,B);A[L1,K]:=0;    
                 HSH2ROW2(1,1,N,L1,L,B[L1,L1],B[L1,L],A,B); 
                 B[L1,L]:=0;L1:=L     
             "END";K:=K1    
         "END"    
     "END"
 "END" HESTGL2;   
           "EOP"  
        
0"CODE" 34605;    
 "PROCEDURE" HSH2COL(LA,LB,U,I,A1,A2,A,B);"VALUE" LA,LB,U,I,A1,A2;  
 "INTEGER" LA,LB,U,I;"REAL" A1,A2;"ARRAY" A,B;  
 "IF"A2^=0 "THEN" 
 "BEGIN" "REAL" R,T,C;"ARRAY" V[I:I+1]; 
     "PROCEDURE" HSHVECMAT(LR,UR,LC,UC,X,U,A);"CODE" 31070; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2) "ELSE" SQRT(A1*A1+A2*A2); 
     T:=A1+R;C:=-T/R;V[I]:=1;V[I+1]:=A2/T;      
     HSHVECMAT(I,I+1,LA,U,C,V,A);HSHVECMAT(I,I+1,LB,U,C,V,B)
 "END" HSH2COL;   
          "EOP"   
        
0"CODE" 34606;    
 "PROCEDURE" HSH3COL(LA,LB,U,I,A1,A2,A3,A,B);   
 "VALUE"LA,LB,U,I,A1,A2,A3;"INTEGER"LA,LB,I,U;"REAL"A1,A2,A3;"ARRAY"A,B;      
 "IF" A2^=0 "OR" A3^=0 "THEN" 
 "BEGIN" "REAL" R,T,C;"ARRAY" V[I:I+2]; 
     "PROCEDURE" HSHVECMAT(LR,UR,LC,UC,X,U,A);"CODE" 31070; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2+A3*A3) 
     "ELSE" SQRT(A1*A1+A2*A2+A3*A3);  
     T:=A1+R;C:=-T/R;V[I]:=1;V[I+1]:=A2/T;V[I+2]:=A3/T;   
     HSHVECMAT(I,I+2,LA,U,C,V,A);HSHVECMAT(I,I+2,LB,U,C,V,B)
 "END" HSH3COL;   
          "EOP"   
        
1SECTION : 3.4.1.2           (JANUARY 1976)                     PAGE 28      
        
        
        
0"CODE" 34607;    
 "PROCEDURE" HSH2ROW3(L,UA,UB,UX,J,A1,A2,A,B,X);"VALUE" L,UA,UB,UX, 
 J,A1,A2;"INTEGER" L,UA,UB,UX,J;"REAL" A1,A2;"ARRAY" A,B,X; 
 "IF"A2^=0 "THEN" 
 "BEGIN""REAL" R,T,C;"INTEGER" K;"ARRAY" V[J:J+1];
     "PROCEDURE" HSHVECTAM(LR,UR,LC,UC,X,U,A);"CODE" 31073; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2) "ELSE" SQRT(A1*A1+A2*A2); 
     T:=A1+R;C:=-T/R;V[J+1]:=1;V[J]:=A2/T;      
     HSHVECTAM(L,UA,J,J+1,C,V,A);HSHVECTAM(L,UB,J,J+1,C,V,B);       
     HSHVECTAM(1,UX,J,J+1,C,V,X)      
 "END" HSH2ROW3;  
          "EOP"   
        
0"CODE" 34608;    
 "PROCEDURE" HSH2ROW2(LA,LB,UA,UB,J,A1,A2,A,B);"VALUE"LA,LB,UA,UB,  
 J,A1,A2;"INTEGER" LA,LB,UA,UB,J;"REAL" A1,A2;"ARRAY" A,B;
 "IF"A2^=0 "THEN" 
 "BEGIN" "REAL" R,T,C;"INTEGER" K;"ARRAY" V[J:J+1];       
     "PROCEDURE" HSHVECTAM(LR,UR,LC,UC,X,U,A);"CODE" 31073; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2) "ELSE" SQRT(A1*A1+A2*A2); 
     T:=A1+R;C:=-T/R;V[J+1]:=1;V[J]:=A2/T;      
     HSHVECTAM(LA,UA,J,J+1,C,V,A);HSHVECTAM(LB,UB,J,J+1,C,V,B)      
 "END" HSH2ROW2;  
          "EOP"   
        
0"CODE" 34609;    
 "PROCEDURE" HSH3ROW3(L,U,UX,J,A1,A2,A3,A,B,X); 
 "VALUE"L,U,UX,J,A1,A2,A3;"INTEGER"L,J,U,UX;"REAL"A1,A2,A3;"ARRAY"A,B,X;      
 "IF" A2^=0 "OR" A3^=0 "THEN" 
 "BEGIN" "REAL" R,T,C;"ARRAY" V[J:J+2];"INTEGER" K;       
     "PROCEDURE" HSHVECTAM(LR,UR,LC,UC,X,U,A);"CODE" 31073; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2+A3*A3) 
     "ELSE" SQRT(A1*A1+A2*A2+A3*A3);  
     T:=A1+R;C:=-T/R;V[J+2]:=1;V[J+1]:=A2/T;V[J]:=A3/T;   
     HSHVECTAM(L,U,J,J+2,C,V,A);HSHVECTAM(L,U,J,J+2,C,V,B); 
     HSHVECTAM(L,UX,J,J+2,C,V,X)      
 "END" HSH3ROW3;  
          "EOP"   
        
0"CODE" 34610;    
 "PROCEDURE" HSH3ROW2(LA,LB,U,J,A1,A2,A3,A,B);  
 "VALUE"LA,LB,U,J,A1,A2,A3;"INTEGER"LA,LB,U,J;"REAL"A1,A2,A3;"ARRAY"A,B;      
 "IF" A2^=0 "OR" A3^=0 "THEN" 
 "BEGIN" "REAL" R,T,C;"ARRAY" V[J:J+2]; 
     "PROCEDURE" HSHVECTAM(LR,UR,LC,UC,X,U,A);"CODE" 31073; 
     R:="IF" A1<0 "THEN" -SQRT(A1*A1+A2*A2+A3*A3) 
     "ELSE" SQRT(A1*A1+A2*A2+A3*A3);  
     T:=A1+R;C:=-T/R;V[J+2]:=1;V[J+1]:=A2/T;V[J]:=A3/T;   
     HSHVECTAM(LA,U,J,J+2,C,V,A);HSHVECTAM(LB,U,J,J+2,C,V,B)
 "END" HSH3ROW2;  
          "EOP"   
        
        
1SECTION : 1.1.4.3            (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR: D.T.WINTER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 751208.
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS PROCEDURES  THAT MULTIPLY A GIVEN MATRIX BY A
     (GENERALIZED) HOUSEHOLDER MATRIX, I.E. A MATRIX M = I + X * U * U',
     WHERE  I  IS THE UNIT MATRIX,   X  A REAL CONSTANT AND  U  A VECTOR
     (CALLED THE HOUSEHOLDER CONSTANT AND HOUSEHOLDER VECTOR,
     RESPECTIVELY)
 
     HSHVECMAT  PREMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX, THE
                HOUSEHOLDER VECTOR HAS BEEN STORED
                IN A ONE-DIMENSIONAL ARRAY.
     HSHCOLMAT  PREMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX,
                THE HOUSEHOLDER VECTOR HAS BEEN STORED AS A COLUMN
                IN A TWO-DIMENSIONAL ARRAY.
     HSHROWMAT  PREMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX,
                THE HOUSEHOLDER VECTOR HAS BEEN STORED AS A ROW
                IN A TWO-DIMENSIONAL ARRAY.
     HSHVECTAM  POSTMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX,
                THE HOUSEHOLDER VECTOR HAS BEEN STORED
                IN A ONE-DIMENSIONAL ARRAY.
     HSHCOLTAM  POSTMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX,
                THE HOUSEHOLDER VECTOR HAS BEEN STORED
                AS A COLUMN IN A TWO-DIMENSIONAL ARRAY.
     HSHROWTAM  POSTMULTIPLIES  A  MATRIX BY A HOUSEHOLDER MATRIX,
                THE HOUSEHOLDER VECTOR HAS BEEN STORED
                AS A ROW IN A TWO-DIMENSIONAL ARRAY.
 
 KEYWORDS:
 
     HOUSEHOLDER MATRIX
     ORTHOGONAL TRANSFORMATION
 
 
 LANGUAGE: ALGOL 60
 
1SECTION : 1.1.4.3            (JANUARY 1976)                      PAGE 2
 
 
 
 SUBSECTION: HSHVECMAT
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHVECMAT(LR, UR, LC, UC, X, U, A);
     "VALUE" LR, UR, LC, UC, X; "INTEGER" LR, UR, LC, UC;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31070;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[LR:UR];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE PREMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
 PROCEDURES USED:
     TAMVEC    = CP34012
     ELMCOLVEC = CP34022
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
 SUBSECTION: HSHCOLMAT
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHCOLMAT(LR, UR, LC, UC, I, X, U, A);
     "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31071;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     I:  <ARITHMETIC EXPRESSION>;
         THE COLUMN INDEX OF THE HOUSEHOLDER VECTOR;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[LR:UR,I:I];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE PREMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
 PROCEDURES USED:
     TAMMAT    = CP34014
     ELMCOL    = CP34023
 
1SECTION : 1.1.4.3            (JANUARY 1976)                      PAGE 3
 
 
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
 SUBSECTION: HSHROWMAT
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHROWMAT(LR, UR, LC, UC, I, X, U, A);
     "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31072;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     I:  <ARITHMETIC EXPRESSION>;
         THE ROW INDEX OF THE HOUSEHOLDER VECTOR;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[I:I,LR:UR];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE PREMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
 PROCEDURES USED:
     MATMAT    = CP34013
     ELMCOLROW = CP34027
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
 SUBSECTION: HSHVECTAM
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHVECTAM(LR, UR, LC, UC, X, U, A);
     "VALUE" LR, UR, LC, UC, X; "INTEGER" LR, UR, LC, UC;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31073;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[LC:UC];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE POSTMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
1SECTION : 1.1.4.3            (JANUARY 1976)                      PAGE 4
 
 
 
 PROCEDURES USED:
     MATVEC    = CP34011
     ELMROWVEC = CP34027
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
 SUBSECTION: HSHCOLTAM
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHCOLTAM(LR, UR, LC, UC, I, X, U, A);
     "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31074;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     I:  <ARITHMETIC EXPRESSION>;
         THE COLUMN INDEX OF THE HOUSEHOLDER VECTOR;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[LC:UC,I:I];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE POSTMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
 PROCEDURES USED:
     MATMAT    = CP34013
     ELMROWCOL = CP34028
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
1SECTION : 1.1.4.3            (JANUARY 1976)                      PAGE 5
 
 
 
 SUBSECTION: HSHROWTAM
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" HSHROWTAM(LR, UR, LC, UC, I, X, U, A);
     "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
     "REAL" X; "ARRAY" U, A;
     "CODE" 31075;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     LR,UR: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER ROW INDICES;
     LC,UC: <ARITHMETIC EXPRESSIONS>;
         THE LOWER AND UPPER COLUMN INDICES;
     I:  <ARITHMETIC EXPRESSION>;
         THE ROW INDEX OF THE HOUSEHOLDER VECTOR;
     X:  <ARITHMETIC EXPRESSION>;
         THE HOUSEHOLDER CONSTANT;
     U:  <ARRAY IDENTIFIER>; "ARRAY" U[I:I,LC:UC];
         THE HOUSEHOLDER VECTOR;
     A:  <ARRAY IDENTIFIER>; "ARRAY" A[LR:UR,LC:UC];
         THE MATRIX TO BE POSTMULTIPLIED BY THE HOUSEHOLDER MATRIX.
 
 PROCEDURES USED:
     MATTAM    = CP34015
     ELMROW    = CP34024
 
 RUNNING TIME: PROPORTIONAL TO (UC-LC+1)*(UR-LR+1)
 
 SOURCE TEXTS:
 
0"CODE" 31070;
 "PROCEDURE" HSHVECMAT(LR, UR, LC, UC, X, U, A);
 "VALUE" LR, UR, LC, UC, X; "INTEGER" LR, UR, LC, UC;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;
     "PROCEDURE" ELMCOLVEC(L, U, I, A, B, X); "CODE" 34022;
     "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO"
     ELMCOLVEC(LR, UR, LC, A, U, TAMVEC(LR, UR, LC, A, U) * X)
 "END";
          "EOP"
0"CODE" 31071;
 "PROCEDURE" HSHCOLMAT(LR, UR, LC, UC, I, X, U, A);
 "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;
     "PROCEDURE" ELMCOL(L, U, I, J, A, B, X); "CODE" 34023;
     "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO"
     ELMCOL(LR, UR, LC, I, A, U, TAMMAT(LR, UR, LC, I, A, U) * X)
 "END";
          "EOP"
1SECTION : 1.1.4.3            (JANUARY 1976)                      PAGE 6
 
 
 
0"CODE" 31072;
 "PROCEDURE" HSHROWMAT(LR, UR, LC, UC, I, X, U, A);
 "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;
     "PROCEDURE" ELMCOLROW(L, U, I, J, A, B, X); "CODE" 34029;
     "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO"
     ELMCOLROW(LR, UR, LC, I, A, U, MATMAT(LR, UR, I, LC, U, A) * X)
 "END";
          "EOP"
0"CODE" 31073;
 "PROCEDURE" HSHVECTAM(LR, UR, LC, UC, X, U, A);
 "VALUE" LR, UR, LC, UC, X; "INTEGER" LR, UR, LC, UC;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;
     "PROCEDURE" ELMROWVEC(L, U, I, A, B, X); "CODE" 34027;
     "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO"
     ELMROWVEC(LC, UC, LR, A, U, MATVEC(LC, UC, LR, A, U) * X)
 "END";
          "EOP"
0"CODE" 31074;
 "PROCEDURE" HSHCOLTAM(LR, UR, LC, UC, I, X, U, A);
 "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;
     "PROCEDURE" ELMROWCOL(L, U, I, J, A, B, X); "CODE" 34028;
     "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO"
     ELMROWCOL(LC, UC, LR, I, A, U, MATMAT(LC, UC, LR, I, A, U) * X)
 "END";
          "EOP"
0"CODE" 31075;
 "PROCEDURE" HSHROWTAM(LR, UR, LC, UC, I, X, U, A);
 "VALUE" LR, UR, LC, UC, I, X; "INTEGER" LR, UR, LC, UC, I;
 "REAL" X; "ARRAY" U, A;
 "BEGIN" "REAL" "PROCEDURE" MATTAM(L, U, I, J, A, B); "CODE" 34015;
     "PROCEDURE" ELMROW(L, U, I, J, A, B, X); "CODE" 34024;
     "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO"
     ELMROW(LC, UC, LR, I, A, U, MATTAM(LC, UC, LR, I, A, U) * X)
 "END";
          "EOP"
1SECTION : 1.5.1              (MARCH 1977)                        PAGE 1
 
 
 
 AUTHORS: D.T.WINTER(A-D,F-I), T.J.DEKKER(E,J)
 
 
 CONTRIBUTOR: J.KOOPMAN(E,J)
 
 
 INSTITUTES: MATHEMATICAL CENTRE,UNIVERSITY OF AMSTERDAM.
 
 
 RECEIVED: 770328
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS  PROCEDURES FOR THE  ELEMENTARY OPERATIONS IN
     DOUBLE PRECISION ARITHMETIC.
 
     A. DPADD  ADDS TWO SINGLE PRECISION NUMBERS TO A DOUBLE PRECISION
               SUM.
     B. DPSUB  SUBTRACTS TWO SINGLE PRECISION NUMBERS TO A DOUBLE
               PRECISION DIFFERENCE.
     C. DPMUL  MULTIPLIES   TWO  SINGLE  PRECISION  NUMBERS   TO  A
               DOUBLE PRECISION PRODUCT.
     D. DPDIV  DIVIDES  TWO SINGLE PRECISION NUMBERS  TO A DOUBLE
               PRECISION QUOTIENT.
     E. DPPOW  COMPUTES A**EXPON IN DOUBLE PRECISION,WHERE
               A IS A SINGLE PRECISION REAL NUMBER AND EXPON THE
               INTEGER EXPONENT.
     F. LNGADD ADDS TWO DOUBLE PRECISION NUMBERS.
     G. LNGSUB SUBTRACTS TWO DOUBLE PRECISION NUMBERS.
     H. LNGMUL MULTIPLIES TWO DOUBLE PRECISION NUMBERS.
     I. LNGDIV DIVIDES TWO DOUBLE PRECISION NUMBERS.
     J. LNGPOW COMPUTES (A,AA)**EXPON IN DOUBLE PRECISION,WHERE
               (A,AA) IS A DOUBLE PRECISION REAL NUMBER AND EXPON
               THE INTEGER EXPONENT.
 
 
 KEYWORDS:
 
     DOUBLE PRECISION ARITHMETIC
     EXPONENTIATION.
 
 
 LANGUAGE: COMPASS(A-D,F-I),ALGOL 60(E,J)
1SECTION : 1.5.1              (DECEMBER 1979)                          PAGE 2
 
 
 
 METHOD AND PERFORMANCE:
 
     THE PROCEDURES A-D,F-I USE THE HARDWARE FUNCTIONS FOR DOUBLE
     PRECISION THAT ARE AVAILIBLE ON THE CYBER.
     THE PROCEDURES  LNG ADD,  LNG SUB,  LNG MUL AND  LNG DIV  CHECK THE
     INPUT PARAMETERS  (A,AA)  AND  (B,BB)  FOR CORRECTNESS. A HEAD/TAIL
     PAIR IS  A CORRECT  DOUBLE  PRECISION  PARAMETER  IN THE  FOLLOWING
     CASES:
     A)  THE TAIL IS ZERO;
     B)  THE EXPONENT  IN THE BINARY  REPRESENTATION  OF THE TAIL  IS 48
         LESS THAN THE EXPONENT OF THE HEAD.
     AN OUTPUT PARAMETER  OF THESE PROCEDURES ALWAYS IS A CORRECT DOUBLE
     PRECISION NUMBER.  IF AN INPUT PARAMETER IS NOT CORRECT,  THE ERROR
     MESSAGE "DP PARAMETER TAIL ERROR" WILL BE ISSUED.
 
     BOTH PROCEDURES E AND J MAKE USE OF THE BINARY REPRESENTATION OF
     THE INTEGER EXPONENT. IF X DENOTES THE NUMBER THAT IS TO BE
     EXPONENTIATED, THE PROCEDURES E AND J RUN AS FOLLOWS:
     THE SEQUENCE X,X**2,X**4,X**8,... IS FORMED WHILE
     SIMULTANEOUSLY THE BINARY REPRESENTATION OF THE EXPONENT
     IS CHECKED; WHEN THE I-TH DIGIT EQUALS ONE,THE FACTOR
     X**(2**(I-1)) IS TAKEN INTO ACCOUNT.
 
 
 EXAMPLE OF USE:
 
     SEE THE PROCEDURE LNGREATODECI (SECTION 1.5.3).
 
1SECTION : 1.5.1              (MARCH 1977)                        PAGE 3
 
 
 
 SUBSECTION: DP ADD
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" DP ADD(A, B, C, CC);
     "VALUE" A, B; "REAL" A, B, C, CC;
     "CODE" 31101;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,B: <ARITHMETIC EXPRESSIONS>;
         THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE DOUBLE PRECISION RESULT OF A+B.
 
 
 SUBSECTION: DP SUB
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" DP SUB(A, B, C, CC);
     "VALUE" A, B; "REAL" A, B, C, CC;
     "CODE" 31102;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,B: <ARITHMETIC EXPRESSIONS>;
         THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE DOUBLE PRECISION RESULT OF A-B.
 
 
 SUBSECTION: DP MUL
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" DP MUL(A, B, C, CC);
     "VALUE" A, B; "REAL" A, B, C, CC;
     "CODE" 31103;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,B: <ARITHMETIC EXPRESSIONS>;
         THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE DOUBLE PRECISION RESULT OF A*B.
 
1SECTION : 1.5.1              (MARCH 1977)                        PAGE 4
 
 
 
 SUBSECTION: DP DIV
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" DP DIV(A, B, C, CC);
     "VALUE" A, B; "REAL" A, B, C, CC;
     "CODE" 31104;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,B: <ARITHMETIC EXPRESSIONS>;
         THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE DOUBLE PRECISION RESULT OF A/B.
 
 
 SUBSECTION: DP POW.
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE"DP POW(A ,EXPON ,C ,CC );
     "VALUE"A,EXPON;"INTEGER"EXPON;"REAL"A,C,CC;
     "CODE"31109;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A       : <ARITHMETIC EXPRESSION>;
               THE NUMBER THAT IS TO BE EXPONENTIATED;
     EXPON   : <ARITHMETIC EXPRESSION>;
               THE (INTEGER) POWER TO WHICH A WILL BE RAISED;
     C , CC  : <REAL VARIABLES>;
               EXIT: THE HEAD (C) AND TAIL (CC) OF THE DOUBLE
                     PRECISION RESULT A**EXPON.
 
 
 PROCEDURES USED:
 
     LNG POW = CP31110.
 
 
 RUNNING TIME:
 
     ROUGHLY PROPORTIONAL TO LN(EXPON).
1SECTION : 1.5.1              (MARCH 1977)                        PAGE 5
 
 
 
 SUBSECTION: LNG ADD
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" LNG ADD(A, AA, B, BB, C, CC);
     "VALUE" A, AA, B, BB; "REAL" A, AA, B, BB, C, CC;
     "CODE" 31105;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,AA,B,BB: <ARITHMETIC EXPRESSIONS>;
         THE HEADS (A AND B) AND THE TAILS (AA AND BB) OF THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE RESULT (A,AA)+(B,BB).
 
 
 SUBSECTION: LNG SUB
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" LNG SUB(A, AA, B, BB, C, CC);
     "VALUE" A, AA, B, BB; "REAL" A, AA, B, BB, C, CC;
     "CODE" 31106;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,AA,B,BB: <ARITHMETIC EXPRESSIONS>;
         THE HEADS (A AND B) AND THE TAILS (AA AND BB) OF THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE RESULT (A,AA)-(B,BB).
 
 
 SUBSECTION: LNG MUL
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" LNG MUL(A, AA, B, BB, C, CC);
     "VALUE" A, AA, B, BB; "REAL" A, AA, B, BB, C, CC;
     "CODE" 31107;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,AA,B,BB: <ARITHMETIC EXPRESSIONS>;
         THE HEADS (A AND B) AND THE TAILS (AA AND BB) OF THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE RESULT (A,AA)*(B,BB).
1SECTION : 1.5.1              (MARCH 1977)                        PAGE 6
 
 
 
 SUBSECTION: LNG DIV
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE" LNG DIV(A, AA, B, BB, C, CC);
     "VALUE" A, AA, B, BB; "REAL" A, AA, B, BB, C, CC;
     "CODE" 31108;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,AA,B,BB: <ARITHMETIC EXPRESSIONS>;
         THE HEADS (A AND B) AND THE TAILS (AA AND BB) OF THE OPERANDS;
     C,CC: <REAL VARIABLES>;
         THE HEAD AND TAIL OF THE RESULT (A,AA)/(B,BB).
 
 
 SUBSECTION: LNG POW.
 
 CALLING SEQUENCE:
 
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
     "PROCEDURE"LNG POW(A ,AA ,EXPON ,C ,CC );
     "VALUE"A,AA,EXPON;"INTEGER"EXPON;"REAL"A,AA,C,CC;
     "CODE"31110;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A,AA    : <ARITHMETIC EXPRESSIONS>;
               THE HEAD (A) AND TAIL (AA) OF THE NUMBER THAT
               IS TO BE EXPONENTIATED;
     EXPON   : <ARITHMETIC EXPRESSION>;
               THE (INTEGER) POWER TO WHICH (A,AA) WILL BE RAISED;
     C,CC    : <REAL VARIABLES>;
               EXIT: THE HEAD (C) AND TAIL (CC) OF THE DOUBLE
                     PRECISION RESULT (A,AA)**EXPON.
 
 
 PROCEDURES USED:
 
     LNG MUL = CP31107.
     LNG DIV = CP31108.
 
 
 RUNNING TIME:
 
     ROUGHLY PROPORTIONAL TO LN(EXPON).
1SECTION : 1.5.1              (DECEMBER 1979)                     PAGE 7
 
 
 
 SOURCE TEXT(S):
 
 ALL PROCEDURES, EXCEPT POW AND LNG POW, ARE WRITTEN IN COMPASS,  IT IS
 NOT POSSIBLE TO SIMULATE THESE PROCEDURES IN ALGOL 60, SO ONLY THE TEXT
 IS GIVEN FOR POW AND LNG POW.
 
 
 
 "CODE"31109;
 "PROCEDURE"DP POW(A,EXPON,C,CC);
 "VALUE"A,EXPON;"INTEGER"EXPON;"REAL"A,C,CC;
 "BEGIN""PROCEDURE"LNG POW(A,AA,EXPON,C,CC);"CODE"31110;
        LNG POW(A,0,EXPON,C,CC)
 "END" DP POW;
    "EOP"
 
 
 "CODE"31110;
 "PROCEDURE"LNG POW(A,AA,EXPON,C,CC);
 "VALUE"A,AA,EXPON;"INTEGER"EXPON;"REAL"A,AA,C,CC;
 "BEGIN""INTEGER"OLDEX,NEWEX;"REAL"D,DD;
        "PROCEDURE"LNG MUL(A,AA,B,BB,C,CC);"CODE"31107;
        "PROCEDURE"LNG DIV(A,AA,B,BB,C,CC);"CODE"31108;
        D:=A;DD:=AA;C:=1;CC:=0;NEWEX:=ABS(EXPON);
        "FOR"OLDEX:=NEWEX"WHILE"NEWEX^=0"DO"
        "BEGIN"NEWEX:=OLDEX//2;
               "IF"NEWEX+NEWEX^=OLDEX
               "THEN"LNG MUL(C,CC,D,DD,C,CC);
               "IF"NEWEX^=0
               "THEN"LNG MUL(D,DD,D,DD,D,DD)
        "END";
        "IF"EXPON<0"THEN"LNG DIV(1,0,C,CC,C,CC)
 "END" LNG POW;
    "EOP"
1SECTION : 6.1               (JANUARY 1976)                      PAGE 1
  
  
  
 AUTHOR: D.T.WINTER.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 751208.
  
  
 BRIEF DESCRIPTION:
     THIS SECTION CONTAINS TWO PROCEDURES:
     1)  PI: DELIVERS THE VALUE OF PI;
     2)  E:  DELIVERS THE VALUE OF E.
  
  
 KEYWORDS:
  
     MATHEMATICAL CONSTANTS
     PI
     E
  
  
 SUBSECTION: PI
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" PI;
     "CODE" 30006;
  
     PI:= THE CONSTANT PI IN 48 BITS PRECISION.
  
 LANGUAGE: COMPASS
  
  
 SUBSECTION: E
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" E;
     "CODE" 30007;
  
     E:= THE CONSTANT E IN 48 BITS PRECISION.
  
 LANGUAGE: COMPASS
1SECTION : 6.1                (JANUARY 1976)                      PAGE 2
  
  
  
 SOURCE TEXTS:
  
  THE SOURCE TEXTS  GIVEN HERE  ARE NOT THE ACTUAL  SOURCE TEXTS,  AS
  THESE PROCEDURES ARE WRITTEN IN COMPASS. EVEN, THE TEXTS GIVEN
  HERE DO NOT GIVE THE SAME RESULTS ON THE CDC CYBER SYSTEM, SINCE
  THE ALGOL-60 COMPILER CANNOT READ THE CONSTANTS IN 48 BIT PECISION.
  
0"CODE" 30006;
 "REAL" "PROCEDURE" PI;
 PI:= 3.14159265358979;
"EOP"  
0"CODE" 30007;
 "REAL" "PROCEDURE" E;
 E:= 2.71828182845905;
"EOP"
1SECTION : 6.2                (JANUARY 1979)                      PAGE 1
 
 
 
 AUTHOR: D.T.WINTER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE,AMSTERDAM.
 
 
 RECEIVED: 751208.
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS SEVEN PROCEDURES:
     A)  MBASE: DELIVERS THE BASE OF THE ARITHMETIC OF THE COMPUTER;
     B)  ARREB: DELIVERS THE ARITHMETIC ERROR BOUND OF THE COMPUTER;
     C)  DWARF: DELIVERS THE SMALLEST  (IN ABSOLUTE VALUE) REPRESENTABLE
             REAL NUMBER;
     D)  GIANT: DELIVERS THE LARGEST REPRESENTABLE REAL NUMBER;
     E)  INTCAP: DELIVERS THE INTEGER CAPACITY;
     F)  OVERFLOW: TESTS WHETHER A VALUE IS AN OVERFLOW VALUE;
     G)  UNDERFLOW: TESTS WHETHER A VALUE IS AN UNDERFLOW VALUE;
     FOR A DETAILED DESCRIPTION SEE METHOD AND PERFORMANCE.
 
 
 KEYWORDS:
 
     ARITHMETIC CONSTANTS
     MACHINE CONSTANTS
     OVERFLOW
     UNDERFLOW
 
 
 SUBSECTION: MBASE
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "INTEGER" "PROCEDURE" MBASE;
     "CODE" 30001;
 
     MBASE:= 2, THE BASE OF THE ARITHMETIC OF THE CYBER.
 
 LANGUAGE: COMPASS
 
 
 SUBSECTION: ARREB
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" ARREB;
     "CODE" 30002;
 
     ARREB:= 2 ** (-47), THE ARITHMETIC RELATIVE ERROR BOUND.
 
 LANGUAGE: COMPASS
 
 
1SECTION : 6.2                (DECEMBER 1979)                     PAGE 2
 
 
 
 SUBSECTION: DWARF
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" DWARF;
     "CODE" 30003;
 
     DWARF:= THE SMALLEST (IN ABSOLUTE VALUE) REPRESENTABLE REAL NUMBER.
 
 LANGUAGE: COMPASS
 
 SUBSECTION: GIANT
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" GIANT;
     "CODE" 30004;
 
     GIANT:= THE LARGEST REPRESENTABLE REAL NUMBER.
 
 LANGUAGE: COMPASS
 
 
 SUBSECTION: INTCAP
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "INTEGER" "PROCEDURE" INTCAP;
     "CODE" 30005;
 
     INTCAP:= THE INTEGER CAPACITY.
 
 LANGUAGE: COMPASS
 
 
 SUBSECTION: OVERFLOW
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "BOOLEAN" "PROCEDURE" OVERFLOW(X); "VALUE" X; "REAL" X;
     "CODE" 30008;
 
     THE MEANING OF THE FORMAL PARAMETER IS:
     X:  <REAL VARIABLE>;
         CONTAINS THE VALUE TO BE TESTED.
 
     OVERFLOW  DELIVERS  "TRUE"  IF  X  CONTAINS AN OVERFLOW VALUE,  AND
     "FALSE" OTHERWISE.
 
 LANGUAGE: COMPASS
1SECTION : 6.2                (DECEMBER 1979)                     PAGE 3
 
 
 
 SUBSECTION: UNDERFLOW
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS:
     "BOOLEAN" "PROCEDURE" UNDERFLOW(X); "VALUE" X; "REAL" X;
     "CODE" 30009;
 
     THE MEANING OF THE FORMAL PARAMETER IS:
     X:  <REAL VARIABLE>;
         CONTAINS THE VALUE TO BE TESTED.
 
     UNDERFLOW DELIVERS  "TRUE"  IF  X CONTAINS AN UNDERFLOW VALUE,  AND
     "FALSE" OTHERWISE.
 
 LANGUAGE: COMPASS
 
 METHOD AND PERFORMANCE:
 
     THE PROCEDURES DELIVER  THE FOLLOWING VALUES,  THAT ARE ESSENTIALLY
     MACHINE DEPENDENT:
     1)  MBASE: 2;
     2)  ARREB: 2**(-47);
     3)  DWARF: 2**48*2**(-1022);
     4)  GIANT: (2**48-1)*2**1022;
     5)  INTCAP: 2**48-2.
     FOR MBASE,  DWARF AND  GIANT THE VALUES  ARE CLEAR,  WE EXPLAIN THE
     OTHERS HERE:
     ARREB: THIS IS THE SMALLEST POSITIVE NUMBER SO THAT 1+ARREB^=1;
     INTCAP: THIS IS  THE LARGEST POSITIVE NUMBER  SO THAT THE FOLLOWING
         BOOLEAN EXPRESSION DELIVERS "TRUE" FOR EVERY INTEGER I:
             "IF" I<0 "OR" I>INTCAP "THEN" "TRUE" "ELSE" I-1^=I;
         THE CORRECT VALUE IS NOT  2**48-1,  AS IN THE  CYBER ARITHMETIC
         I=J IF I=2**48 AND J=2**48-1.
     WARNING: DWARF IS NOT VERY USEFUL WHEN TRAPPING UNDERFLOW VALUES:
         ABS(X) >= DWARF  NEARLY ALWAYS  DELIVERS TRUE EVEN IF ABS(X) IS
         SMALLER  THEN  DWARF  DUE  TO THE ARITHMETIC.  ONE  SHOULD USE:
         ABS(X) > DWARF  (AND  ONE  TRAPS  NON-UNDERFLOW VALUES TOO)  OR
         THE PROCEDURE UNDERFLOW.
     NOTE: AS THE ALGOL 60 ERRORMESSAGE "ARITHMETIC OVERFLOW"
         IS NOT ISSUED AT THE MOMENT THE OVERFLOW VALUE
         IS CREATED BUT WHEN SUCH A VALUE IS USED, THE
         PROCEDURE OVERFLOW IS WELL-DEFINED.
 
 
1SECTION : 6.2                (DECEMBER 1979)                     PAGE 4
 
 
 
 EXAMPLE OF USE:
 
 
     HERE WE GIVE  AN EXAMPLE  OF USE  OF THE  PROCEDURES  OVERFLOW  AND
     UNDERFLOW:
 
     "BEGIN"
         "BOOLEAN" "PROCEDURE" OVERFLOW(X); "CODE" 30009;
         "BOOLEAN" "PROCEDURE" UNDERFLOW(X); "CODE" 30008;
         "REAL" "PROCEDURE" DWARF; "CODE" 30003;
         "REAL" X, Y;
         Y:= 0; X:= 1 / Y;
         "IF" OVERFLOW(X) "THEN" OUTPUT(61, "(""("OVERFLOW")", /")");
         X:= DWARF;
         "IF" "NOT" UNDERFLOW(X) "THEN"
         OUTPUT(61, "(""("NO UNDERFLOW WITH DWARF")", /")");
         X:= X / 2;
         "IF" X ^= 0 "THEN"
         OUTPUT(61, "(""("DWARF / 2 ^= 0")", /")");
         "IF" UNDERFLOW(X) "THEN"
         OUTPUT(61, "(""("DWARF / 2 IS UNDERFLOW")", /")");
         "IF" X * 2 = 0 "THEN"
         OUTPUT(61, "(""("BECAUSE (DWARF / 2) * 2 = 0")", /")")
     "END"
 
     RESULTS:
     OVERFLOW
     NO UNDERFLOW WITH DWARF
     DWARF / 2 ^= 0
     DWARF / 2 IS UNDERFLOW
     BECAUSE (DWARF / 2) * 2 = 0
 
1SECTION : 6.2                (JANUARY 1976)                      PAGE 5
 
 
 
 SOURCE TEXTS:
 
 THESE  ARE NOT  THE ACTUAL  SOURCE TEXTS,  AS THESE  PROCEDURES ARE
 WRITTEN IN COMPASS,  MOREOVER,  THE RESULTS NEED NOT BE THAT OF THE
 ACTUAL PROCEDURES.
 
 "CODE" 30001;
 "INTEGER" "PROCEDURE" MBASE;
 MBASE:= 2;
 
 "CODE" 30002;
 "REAL" "PROCEDURE" ARREB;
 ARREB:= 2**(-47);
 
 "CODE" 30003;
 "REAL" "PROCEDURE" DWARF;
 DWARF:= 2**48*2**(-1022);
 
 "CODE" 30004;
 "REAL" "PROCEDURE" GIANT;
 GIANT:= (2**48-1)*2**1022;
 
 "CODE" 30005;
 "INTEGER" "PROCEDURE" INTCAP;
 INTCAP:= 2**48-2;
 
     THE SOURCE TEXT OF  OVERFLOW AND UNDERFLOW  ARE NOT GIVEN HERE,  AS
     THESE EVEN CANNOT BE SIMULATED IN ALGOL-60.
1SECTION : 0.0               (JULY 1979)                         PAGE 0       
1SECTION : 0.0               (JULY 1979)                         PAGE 0       
1SECTION : 0.0               (JULY 1979)                         PAGE 0       
1SECTION : 0.0               (JULY 1979)                         PAGE 0