1SECTION : 3.3.1.1.3.3        (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR/CONTRIBUTOR: J.J.G. ADMIRAAL.
 
 INSTITUTE: UNIVERSITY OF AMSTERDAM.
 
 RECEIVED: 761101.
 
 BRIEF DESCRIPTION:
     THE PROCEDURE SYMEIGIMP IMPROVES A GIVEN APPROXIMATION OF
     A REAL SYMMETRIC EIGENSYSTEM AND CALCULATES ERROR BOUNDS
     FOR THE EIGENVALUES.
 
 KEYWORDS:
     EIGENVALUES.
     EIGENVECTORS.
     SYMMETRIC MATRIX.
     RAYLEIGH QUOTIENTS.
     ERROR BOUNDS.
     IMPROVED EIGENSYSTEM.
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE IS :
 
     "PROCEDURE" SYMEIGIMP(N,A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX);
     "VALUE" N;"INTEGER" N;
     "ARRAY" A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX;
     "CODE" 36401;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF MATRIX A;
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N] CONTAINS A REAL SYMMETRIC MATRIX
             WHOSE EIGENSYSTEM HAS TO BE IMPROVED;
     VEC:    <ARRAY IDENTIFIER>;
             "ARRAY" VEC[1:N,1:N] CONTAINS A MATRIX WHOSE COLUMNS ARE
             A SYSTEM OF APPROXIMATE EIGENVECTORS OF MATRIX A;
             ENTRY: INITIAL APPROXIMATIONS;
             EXIT: IMPROVED APPROXIMATIONS;
     VAL1:   <ARRAY IDENTIFIER>;
             "ARRAY" VAL1[1:N];
             ENTRY: INITIAL APPROXIMATIONS OF THE EIGENVALUES OF A;
             EXIT: THE HEAD PARTS OF THE DOUBLE PRECISION IMPROVED
             APPROXIMATIONS OF THE EIGENVALUES OF A;
     VAL2:   <ARRAY IDENTIFIER>;
             "ARRAY" VAL2[1:N];
             EXIT: THE TAIL PARTS OF THE DOUBLE PRECISION
                   IMPROVED EIGENVALUES OF A;
1SECTION : 3.3.1.1.3.3        (DECEMBER 1979)                     PAGE 2
 
 
 
     LBOUND,
     UBOUND: <ARRAY IDENTIFIER>;
             EXIT: "ARRAY" LBOUND, UBOUND [1:N] CONTAIN THE LOWER
                   AND UPPER ERRORBOUNDS RESPECTIVELY FOR THE EIGENVALUE
                   APPROXIMATIONS IN VAL1,VAL2[1:N] SUCH THAT THE
                   I-TH EXACT EIGENVALUE LIES BETWEEN VAL1[I]+VAL2[I]
                   -LBOUND[I] AND VAL1[I]+VAL2[I]+UBOUND[I];
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[0:5];
             ENTRY: AUX[0]= THE RELATIVE PRECISION OF THE ELEMENTS OF A;
                    AUX[2]= THE RELATIVE TOLERANCE FOR THE RESIDUAL MATRIX;
                            THE ITERATION ENDS WHEN THE MAXIMUM ABSOLUTE
                            VALUE OF THE RESIDU ELEMENTS IS SMALLER THAN
                            AUX[2]*AUX[1].
                    AUX[4]= THE MAXIMUM NUMBER OF ITERATIONS ALLOWED;
             EXIT:  AUX[1]= INFINITY NORM OF THE MATRIX A;
                    AUX[3]= MAXIMUM ABSOLUTE ELEMENT OF THE RESIDUAL MATRIX;
                    AUX[5]= NUMBER OF ITERATIONS;
 
 PROCEDURES USED:
     LNGMATVEC = CP34411,
     LNGMATMAT = CP34413,
     LNGTAMMAT = CP34414,
     VECVEC    = CP34010,
     MATMAT    = CP34013,
     TAMMAT    = CP34014.
     MERGESORT = CP36405,
     VECPERM   = CP36404,
     ROWPERM   = CP36403,
     ORTHOG    = CP36402,
     QRISYM    = CP34163,
     INFNRMMAT = CP31064.
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.
 
 REQUIRED CENTRAL MEMORY:
     AUXILIARY ARRAYS ARE DECLARED TO A TOTAL OF 3*N*N + 6*N REALS
     AND N INTEGERS; MOREOVER, N INTEGERS OR N BOOLEANS ARE USED
     BY MERGESORT, VECPERM AND ROWPERM.
 
 METHOD AND PERFORMANCE: SEE[1].
1SECTION : 3.3.1.1.3.3        (DECEMBER 1979)                     PAGE 3
 
 
 
 REFERENCES:
     [1]. J.J.G. ADMIRAAL.
          ITERATIEF VERBETEREN VAN REEEL SYMMETRISCH EIGENSYSTEEM
          EN BEREKENEN VAN FOUTGRENZEN VOOR DE VERKREGEN EIGENWAARDEN.
          DOCTORAL SCRIPTION,MARCH 1976,
          UNIVERSITEIT VAN AMSTERDAM.
     [2]. R.T. GREGORY AND D.L. KARNEY.
          A COLLECTION OF MATRICES FOR TESTING COMPUTATIONAL
          ALGORITHMS,
          WILEY-INTERSCIENCE, 1969.
 
 EXAMPLE OF USE.
 
 "BEGIN" "INTEGER" I,J;"REAL" S;
    "ARRAY" A,X[1:4,1:4],VAL1,VAL2,LBOUND,UBOUND[1:4],EM,AUX[0:5];
    "PROCEDURE" SYMEIGIMP(N,A,X,VAL1,VAL2,LBOUND,UBOUND,AUX);
    "CODE" 36401;
    "INTEGER" "PROCEDURE" QRISYM(A,N,VAL,EM);"CODE" 34163;
    A[1,1]:=A[2,2]:=A[3,3]:=A[4,4]:=6;
    A[1,2]:=A[2,1]:=A[3,1]:=A[1,3]:=4;
    A[4,2]:=A[2,4]:=A[3,4]:=A[4,3]:=4;
    A[1,4]:=A[4,1]:=A[3,2]:=A[2,3]:=1;
    "FOR" I:=1 "STEP" 1 "UNTIL" 4 "DO"
    "FOR" J:=I "STEP" 1 "UNTIL" 4 "DO" X[I,J]:=X[J,I]:=A[I,J];
    OUTPUT(61,"(""("A")",/,4(4(+DB),/)")",A);
    EM[0]:="-14;EM[4]:=100;EM[2]:="-5;
    QRISYM(X,4,VAL1,EM);
    AUX[0]:=0;AUX[4]:=10;AUX[2]:="-14;
    SYMEIGIMP(4,A,X,VAL1,VAL2,LBOUND,UBOUND,AUX);
    OUTPUT(61,"("/,"("THE EXACT EIGENVALUES ARE: -1 , +5 , +5 , +15")",
    //,
    "("THE DIFFERENCES BETWEEN THE CALCULATED AND THE EXACT EIGENVALUES"
    )",//,4(N,/)")",(VAL1[1]+1)+VAL2[1],(VAL1[2]-5)+VAL2[2],(VAL1[3]-
    5)+VAL2[3],(VAL1[4]-15)+VAL2[4]);
    OUTPUT(61,"("/,"("LOWERBOUNDS  UPPERBOUNDS")",//")");
    "FOR" I:=1 "STEP" 1 "UNTIL" 4 "DO"
    OUTPUT(61,"("2(+D.D"+DD5B),/")",LBOUND[I],UBOUND[I]);
    OUTPUT(61,"("/,"("NUMBER OF ITERATIONS = ")",ZD//,
    "("INFINITY NORM OF A = ")",ZD//,
    "("MAXIMUM ABSOLUTE ELEMENT OF RESIDU = ")",D.D"+DD")",
    AUX[5],AUX[1],AUX[3])
 "END"  EXAMPLE OF USE
1SECTION : 3.3.1.1.3.3        (NOVEMBER 1976)                     PAGE 4
 
 
 
 DELIVERS:
 
 A
 +6 +4 +4 +1
 +4 +6 +1 +4
 +4 +1 +6 +4
 +1 +4 +4 +6
 
 THE EXACT EIGENVALUES ARE: -1 , +5 , +5 , +15
 
 THE DIFFERENCES BETWEEN THE CALCULATED AND THE EXACT EIGENVALUES
 
 -6.3423147029256"-022
 +5.5934784498910"-018
 +4.0389678347316"-028
 -5.5947317864427"-018
 
 LOWERBOUNDS  UPPERBOUNDS
 
 +1.2"-23     +1.2"-23
 +7.5"-09     +7.5"-09
 +1.0"-13     +1.0"-13
 +5.6"-18     +5.6"-18
 
 NUMBER OF ITERATIONS =  2
 
 INFINITY NORM OF A = 15
 
 MAXIMUM ABSOLUTE ELEMENT OF RESIDU = 2.8"-14
1SECTION : 3.3.1.1.3.3        (NOVEMBER 1976)                     PAGE 5
 
 
 
 SOURCETEXT:
 
 "CODE" 36401;
 "PROCEDURE" SYMEIGIMP(N,A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX);
 "VALUE" N;"INTEGER" N;"ARRAY" A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX;
 "BEGIN"
     "PROCEDURE" ORTHOG(N,LC,UC,X);"CODE" 36402;
     "PROCEDURE" MERGESORT(VEC1,VEC2,LOW,UPP);"CODE" 36405;
     "PROCEDURE" VECPERM(PERM,LOW,UPP,VECTOR);"CODE" 36404;
     "PROCEDURE" ROWPERM(PERM,LOW,UPP,MATRIX);"CODE" 36403;
     "INTEGER" K,I,J,I0,I1,ITER,MAXITP1;"REAL" S,HEAD,TAIL,MAX,TOL,
     MATEPS,RELERRA,RELTOLR,NORMA;"INTEGER" "ARRAY" PERM[1:N];
     "ARRAY" R,P,Y[1:N,1:N],RQ,RQT,EPS,Z,VAL3,ETA[1:N];
     "PROCEDURE" BOUNDS(I0,I1,N,LBOUND,UBOUND);"VALUE" I0,I1,N;
     "INTEGER" I0,I1,N;"ARRAY" LBOUND,UBOUND;
     "BEGIN" "INTEGER" K,I,J,I01;"REAL" EPS2,DL,DR;
         "FOR" I:=I0,I01 "WHILE" I<=I1 "DO"
         "BEGIN" J:=I01:=I;
             "FOR" J:=J+1 "WHILE" "IF" J>I1 "THEN" "FALSE" "ELSE"
             RQ[J]-RQ[J-1]<=EPS[J]+EPS[J-1] "DO" I01:=J;
             "IF" I = I01 "THEN"
             "BEGIN"
                 "IF" I<N "THEN"
                 "BEGIN"
                     "IF" I=1 "THEN" DL:=DR:=RQ[I+1]-RQ[I]-EPS[I+1]
                     "ELSE" "BEGIN" DL:=RQ[I]-RQ[I-1]-EPS[I-1];
                                    DR:=RQ[I+1]-RQ[I]-EPS[I+1]
                            "END"
                 "END" "ELSE" DL:=DR:=RQ[I]-RQ[I-1]-EPS[I-1];
                 EPS2:=EPS[I]*EPS[I];LBOUND[I]:=EPS2/DR+MATEPS;
                 UBOUND[I]:=EPS2/DL+MATEPS
             "END" "ELSE"
             "BEGIN" "FOR" K:=I "STEP" 1 "UNTIL" I01 "DO"
                 LBOUND[K]:=UBOUND[K]:=EPS[K]+MATEPS
             "END";I01:=I01+1
         "END"
     "END" BOUNDS;
     "PROCEDURE" LNGMATVEC(L,U,I,A,B,C,CC,D,DD);"CODE" 34411;
     "PROCEDURE" LNGTAMMAT(L,U,I,J,A,B,C,CC,D,DD);"CODE" 34414;
     "PROCEDURE" LNGMATMAT(L,U,I,J,A,B,C,CC,D,DD);"CODE" 34413;
     "INTEGER" "PROCEDURE" QRISYM(A,N,VAL,EM);"CODE" 34163;
     "REAL""PROCEDURE" VECVEC(L,U,SHIFT,A,B);"CODE" 34010;
     "REAL""PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;
     "REAL""PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014;
     "REAL" "PROCEDURE" INFNRMMAT(LC,UC,LR,UR,K,A);"CODE" 31064;
     "BOOLEAN" STOP;STOP:="FALSE";NORMA:=INFNRMMAT(1,N,1,N,S,A);
     RELERRA:=AUX[0];RELTOLR:=AUX[2];MAXITP1:=AUX[4]+1;
     MATEPS:=RELERRA*NORMA;TOL:=RELTOLR*NORMA;
     "FOR" ITER:=1 "STEP" 1 "UNTIL" MAXITP1 "DO"
     "BEGIN" STOP:="TRUE";MAX:=0;
                                                               "COMMENT"
1SECTION : 3.3.1.1.3.3        (NOVEMBER 1976)                     PAGE 6
                                                                       ;
 
 
         "FOR" J:=1 "STEP" 1 "UNTIL" N "DO"
         "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN"
             LNGMATVEC(J,J,I,VEC,VAL1,0,0,HEAD,TAIL);
             LNGMATMAT(1,N,I,J,A,VEC,-HEAD,-TAIL,R[I,J],TAIL);
             "IF"ABS(R[I,J])>MAX "THEN" MAX:=ABS(R[I,J])
         "END";"IF" MAX > TOL "THEN" STOP:="FALSE";
         "IF" "NOT" STOP "AND" ITER<MAXITP1 "THEN"
         "BEGIN"
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
             LNGTAMMAT(1,N,I,I,VEC,R,VAL1[I],0,RQ[I],RQT[I]);
             "FOR" J:=1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
                 ETA[I]:=R[I,J]-(RQ[J]-VAL1[J])*VEC[I,J];
                 Z[J]:=SQRT(VECVEC(1,N,0,ETA,ETA))
             "END";
             MERGESORT(RQ,PERM,1,N);VECPERM(PERM,1,N,RQ);
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" EPS[I]:=Z[PERM[I]];VAL3[I]:=VAL1[PERM[I]];
                 ROWPERM(PERM,1,N,I,VEC);ROWPERM(PERM,1,N,I,R)
             "END";
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
             "FOR" J:=I "STEP" 1 "UNTIL" N "DO"
             P[I,J]:=P[J,I]:=TAMMAT(1,N,I,J,VEC,R);
         "END";
         "FOR"I0:=1,I1+1 "WHILE" I0<=N "DO"
         "BEGIN" J:=I1:=I0;
             "FOR" J:=J+1 "WHILE" "IF" J>N "THEN" "FALSE" "ELSE"
             RQ[J]-RQ[J-1]<=SQRT((EPS[J]+EPS[J-1])*NORMA) "DO" I1:=J;
             "IF" STOP "OR" ITER=MAXITP1 "THEN"
             BOUNDS(I0,I1,N,LBOUND,UBOUND) "ELSE"
             "BEGIN"
                 "IF" I0=I1 "THEN"
                 "BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"
                     "IF" K=I0 "THEN" Y[K,I0]:=1 "ELSE"
                     R[K,I0]:=P[K,I0];
                     VAL1[I0]:=RQ[I0];VAL2[I0]:=RQT[PERM[I0]]
                 "END" "ELSE"
                 "BEGIN""INTEGER" N1,I0M1,I1P1;"REAL" M1;"ARRAY"EM[0:5];
                     N1:=I1-I0+1;EM[0]:=EM[2]:="-14;EM[4]:=10*N1;
                     "BEGIN" "ARRAY" PP[1:N1,1:N1],VAL4[1:N1];M1:=0;
                         "FOR" K:=I0 "STEP" 1 "UNTIL" I1 "DO"
                         M1:=M1+VAL3[K];M1:=M1/N1;
                                                               "COMMENT"
1SECTION : 3.3.1.1.3.3        (DECEMBER 1979)                     PAGE 7
                                                                       ;
 
 
                         "FOR" I:=1 "STEP" 1 "UNTIL" N1 "DO"
                         "FOR" J:=1 "STEP" 1 "UNTIL" N1 "DO"
                         "BEGIN" PP[I,J]:=P[I+I0-1,J+I0-1];
                             "IF" I=J "THEN"
                             PP[I,J]:=PP[I,J]+VAL3[J+I0-1]-M1
                         "END";"FOR" I:=I0 "STEP" 1 "UNTIL" I1 "DO"
                         "BEGIN" VAL3[I]:=M1;VAL1[I]:=RQ[I];
                             VAL2[I]:=RQT[PERM[I]]
                         "END";
                         QRISYM(PP,N1,VAL4,EM);
                         MERGESORT(VAL4,PERM,1,N1);
                         "FOR" I:=1 "STEP" 1 "UNTIL" N1 "DO"
                         "FOR" J:=1 "STEP" 1 "UNTIL" N1 "DO"
                         P[I+I0-1,J+I0-1]:=PP[I,PERM[J]];
                         I0M1:=I0-1;I1P1:=I1+1;
                         "FOR" J:=I0 "STEP" 1 "UNTIL" I1 "DO"
                         "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" I0M1,
                             I1P1 "STEP" 1 "UNTIL" N "DO"
                             "BEGIN" S:=0;
                                 "FOR" K:=I0 "STEP" 1 "UNTIL" I1 "DO"
                                 S:=S+P[I,K]*P[K,J];
                                 R[I,J]:=S
                             "END";"FOR" I:=I0 "STEP" 1 "UNTIL" I1 "DO"
                             Y[I,J]:=P[I,J]
                         "END" FOR J
                     "END" INNERBLOCK
                 "END" I1>I0
             "END" NOT STOP
         "END" FOR I0;
         "IF" "NOT" STOP "AND" ITER<MAXITP1 "THEN"
         "BEGIN"
             "FOR" J:=1 "STEP" 1 "UNTIL" N "DO"
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
             "IF" VAL3[I]^=VAL3[J] "THEN"
             Y[I,J]:=R[I,J]/(VAL3[J]-VAL3[I]);
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" N "DO"
                 Z[J]:=MATMAT(1,N,I,J,VEC,Y);
                 "FOR" J:=1 "STEP" 1 "UNTIL" N "DO" VEC[I,J]:=Z[J]
             "END";ORTHOG(N,1,N,VEC)
         "END" "ELSE"
         "BEGIN" AUX[5]:=ITER-1;"GOTO" EXIT "END"
     "END" FOR ITER;
     EXIT: AUX[1]:=NORMA;AUX[3]:=MAX
 "END" SYMEIGIMP;
         "EOP"
1SECTION : 3.1.1.1.1.3.1     ( DECEMBER 1978 )                  PAGE 1
        
        
        
 AUTHORS : J.R.BUNCH,L,KAUFMAN,B.N.PARLETT.     
        
        
 CONTRIBUTOR : C.H.CONVENT. 
        
        
 INSTITUTE : UNIVERSITY OF AMSTERDAM. 
        
        
 RECEIVED : 770712. 
        
        
 BRIEF DESCRIPTION :
        
    DECSYM2 CALCULATES THE LDL' DECOMPOSITION OF A SYMMETRIC MATRIX.
    THE MATRIX MAY BE INDEFINITE AND/OR SINGULAR; 
        
        
 KEYWORDS :       
        
    GENERAL SYMMETRIC MATRIX, 
    LDL' DECOMPOSITION,     
    BLOCK DIAGONAL PIVOTING;
        
        
 CALLING SEQUENCE : 
        
    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" DECSYM2(A,N,TOL,AUX,P,DETAUX);  
    "VALUE" N;"INTEGER" N;"REAL" TOL; 
    "ARRAY" A,DETAUX;"INTEGER" "ARRAY" AUX,P;   
    "CODE" 34291; 
        
    THE MEANING OF THE FORMAL PARAMETERS IS :   
    A :     <ARRAY IDENTIFIER>;       
            "ARRAY" A[1:N,1:N];       
            ENTRY : THE SYMMETRIC COEFFICIENT MATRIX;     
            EXIT  : THE ELEMENTS OF THE LDL' DECOMPOSITION OF A ARE 
                    STORED IN THE UPPER TRIANGULAR PART OF A. HERE D
                    IS A BLOCK DIAGONAL MATRIX WITH BLOCKS OF ORDER 1 
                    OR 2. FOR A BLOCK OF ORDER 2 WE ALWAYS HAVE     
                    D[I,I+1]^=0 AND L[I+1,I]=0,SO THAT D AND L' FIT 
                    IN THE UPPER TRIANGULAR PART OF A.    
                    THE STRICTLY LOWER TRIANGULAR PART OF A IS LEFT 
                    UNDISTURBED.      
    N :     <ARITHMETIC EXPRESSION>;  
            THE ORDER OF THE MATRIX;  
1SECTION : 3.1.1.1.1.3.1     ( DECEMBER 1978 )                  PAGE 2
        
        
        
    TOL :  <ARITHMETIC EXPRESSION>;   
           A RELATIVE TOLERANCE, USED TO CONTROL THE      
           CALCULATION OF THE BLOCK DIAGONAL ELEMENTS;    
    AUX :   <ARRAY IDENTIFIER>;       
            "INTEGER" "ARRAY" AUX[2:5]; 
            EXIT :
            AUX[2] : IF THE MATRIX IS SYMMETRIC THEN 1, OTHER-      
                     WISE 0;IN THE LAST CASE NO DECOMPOSITION       
                     IS PERFORMED;    
            AUX[3] : IF THE MATRIX IS SYMMETRIC THEN THE  
                     NUMBER OF ITS POSITIVE EIGENVALUES,  
                     OTHERWISE 0. IF AUX[3]=N THEN THE    
                     MATRIX IS POSITIVE DEFINITE; 
            AUX[4] : IF THE MATRIX IS SYMMETRIC THEN THE  
                     NUMBER OF ITS NEGATIVE EIGENVALUES,  
                     OTHERWISE 0. IF AUX[4]=N THEN THE    
                     MATRIX IS NEGATIVE DEFINITE; 
            AUX[5] : IF THE MATRIX IS SYMMETRIC THEN THE  
                     NUMBER OF ITS ZERO EIGENVALUES, OTHERWISE      
                     N; SO, IF AUX[5]=0 THEN THE MATRIX IS
                     SYMMETRIC AND NON-SINGULAR;
    P :     <ARRAY IDENTIFIER>;       
            "INTEGER" "ARRAY" P[1:N]; 
            EXIT :
            A VECTOR RECORDING
            1) THE INTERCHANGES PERFORMED ON A DURING THE 
               COMPUTATION OF THE DECOMPOSITION AND       
            2) THE BLOCK STRUCTURE OF D.
            IF P[I]>0 AND P[I+1]=0 A 2*2 BLOCK HAS BEEN   
            FOUND I.E. D[I,I+1]^=0 AND L[I+1,I]=0;
    DETAUX : <ARRAY IDENTIFIER>;      
             "ARRAY" DETAUX[1:N];     
             EXIT : 
             IF P[I]>0 AND P[I+1]>0 THEN DETAUX[I] EQUALS 
             THE EXIT-VALUE OF A[I,I].
             IF P[I]>0 AND P[I+1]=0 THEN DETAUX[I]=1 AND  
             DETAUX[I+1] EQUALS THE VALUE OF THE DETERMINANT
             OF THE CORRESPONDING 2*2 DIAGONAL BLOCK AS   
             DETERMINED BY DECSYM2;   
        
        
 PROCEDURES USED :
        
    ELMROW=CP34030. 
    ICHROW=CP34032. 
    ICHROWCOL=CP34033.      
1SECTION : 3.1.1.1.1.3.1     ( DECEMBER 1978 )                  PAGE 3
        
        
        
 RUNNING TIME : ROUGHLY PROPORTIONAL TO N**3.   
        
        
 METHOD AND PERFORMANCE :   
        
    THE PROCEDURE DECSYM2 COMPUTES THE LDL' DECOMPOSITION OF A      
    SYMMETRIC MATRIX,ACCORDING TO A METHOD DUE TO BUNCH,KAUFMAN AND 
    PARLETT (SEE [1],[2]). IT USES BLOCK DIAGONAL PIVOTING. 
    THE BLOCK DIAGONAL MATRIX D IS DELIVERED IN THE BLOCK DIAGONAL  
    OF A. IF P[I]>0 AND P[I+1]=0 A 2*2 BLOCK HAS BEEN FOUND AND     
    FURTHERMORE : L[I+1,I]=0 WHEN D[I,I+1]^=0.  
    THE STRICTLY UPPER TRIANGULAR PART OF L' IS DELIVERED IN THE    
    STRICTLY UPPER TRIANGULAR PART OF A. FOR THE INERTIA PROBLEM    
    IT IS IMPORTANT THAT DECSYM2 CAN ACCEPT SINGULAR MATRICES.      
    NOTE, HOWEVER, THAT IN ORDER TO FIND THE NUMBER OF ZERO 
    EIGENVALUES OF SINGULAR MATRICES, THE SINGULAR VALUE  
    DECOMPOSITION MIGHT BE PREFERRED. 
    BEFORE THE DECOMPOSITION IS PERFORMED A CHECK IS MADE TO SEE    
    WHETHER THE MATRIX IS SYMMETRIC. IF THE MATRIX IS ASYMMETRIC    
    THEN NO DECOMPOSITION IS PERFORMED; 
        
        
 REFERENCES :     
        
    1) J.R.BUNCH,L.KAUFMAN. 
       SOME STABLE METHODS FOR CALCULATING INERTIA AND SOLVING      
       SYMMETRIC LINEAR SYSTEMS.      
       MATHEMATICS OF COMPUTATION 31,P 163-180,1977.      
    2) J.R.BUNCH,L.KAUFMAN,B.N.PARLETT. 
       DECOMPOSITION OF A SYMMETRIC MATRIX.     
       NUMERISCHE MATHEMATIK 27,P 95-109,1976.  
1SECTION : 3.1.1.1.1.3.1     ( DECEMBER 1978 )                  PAGE 4
        
        
        
 SOURCE TEXT :    
 "CODE" 34291;    
 "PROCEDURE" DECSYM2(A,N,TOL,AUX,P,DETAUX);     
 "VALUE" N;"INTEGER" N;"REAL" TOL;    
 "ARRAY" A,DETAUX;"INTEGER" "ARRAY" P,AUX;      
 "BEGIN" "INTEGER" I,J,K,M,IP1,IP2,DUMMY;"BOOLEAN" ONEBYONE,SYM;    
    "REAL" DET,S,T,ALPHA,LAMBDA,SIGMA,AII,AIP1,AIP1I;     
    "PROCEDURE" ELMROW(L,U,I,J,A,B,X);"CODE" 34024;       
    "PROCEDURE" ICHROW(L,U,I,J,A);"CODE" 34032; 
    "PROCEDURE" ICHROWCOL(L,U,I,J,A);"CODE" 34033;
    AUX[3]:=AUX[4]:=0;SYM:="TRUE";I:=0; 
    "FOR" DUMMY:=0 "WHILE" SYM "AND" (I<N) "DO" 
    "BEGIN" I:=I+1;J:=I;    
       "FOR" M:=0 "WHILE" SYM "AND" (J<N) "DO"  
       "BEGIN" J:=J+1;SYM:=SYM "AND" (A[I,J]=A[J,I]) "END"; 
    "END";
    "IF" SYM "THEN" AUX[2]:=1 
         "ELSE" "BEGIN" AUX[2]:=0;"GOTO" ENDDEC "END";    
    ALPHA:=(1+SQRT(17))/8;P[N]:=N;I:=1; 
    "FOR" DUMMY:=0 "WHILE" I<N "DO"   
    "BEGIN" IP1:=I+1;IP2:=I+2;AII:=ABS(A[I,I]);P[I]:=I;   
       LAMBDA:=ABS(A[I,IP1]);J:=IP1;  
       "FOR" M:=IP2 "STEP" 1 "UNTIL" N "DO"     
       "IF" ABS(A[I,M])>LAMBDA "THEN" 
       "BEGIN" J:=M;LAMBDA:=ABS(A[I,M]) "END";  
       T:=ALPHA*LAMBDA;ONEBYONE:="TRUE";
       "IF" AII<T "THEN"    
       "BEGIN" SIGMA:=LAMBDA; 
          "FOR" M:=IP1 "STEP" 1 "UNTIL" J-1 "DO"
             "IF" ABS(A[M,J])>SIGMA "THEN" SIGMA:=ABS(A[M,J]);      
          "FOR" M:=J+1 "STEP" 1 "UNTIL" N "DO"  
             "IF" ABS(A[J,M])>SIGMA "THEN" SIGMA:=ABS(A[J,M]);      
          "IF" SIGMA*AII<LAMBDA "THEN"
          "BEGIN" "IF" ALPHA*SIGMA<ABS(A[J,J]) "THEN"     
             "BEGIN" ICHROW(J+1,N,I,J,A);ICHROWCOL(IP1,J-1,I,J,A);  
                T:=A[I,I];A[I,I]:=A[J,J];A[J,J]:=T;P[I]:=J
             "END"
             "ELSE" 
             "BEGIN" "IF" J>IP1 "THEN"
               "BEGIN" ICHROW(J+1,N,IP1,J,A);ICHROWCOL(IP2,J-1,IP1,J,A);      
                   T:=A[I,I];A[I,I]:=A[J,J];A[J,J]:=T;    
                   T:=A[I,J];A[I,J]:=A[I,IP1];A[I,IP1]:=T 
                "END";      
                DET:=A[I,I]*A[IP1,IP1]-A[I,IP1]**2;AIP1I:=A[I,IP1]/DET;       
                AII:=A[I,I]/DET;AIP1:=A[IP1,IP1]/DET;P[I]:=J;P[IP1]:=0;       
                DETAUX[I]:=1;DETAUX[IP1]:=DET;                 "COMMENT"      
1SECTION : 3.1.1.1.1.3.1     ( DECEMBER 1978 )                  PAGE 5
                                                                  ; 
        
        
                "FOR" J:=IP2 "STEP" 1 "UNTIL" N "DO"      
                "BEGIN" S:=AIP1I*A[IP1,J]-AIP1*A[I,J];    
                   T:=AIP1I*A[I,J]-AII*A[IP1,J];ELMROW(J,N,J,I,A,A,S);
                   ELMROW(J,N,J,IP1,A,A,T);A[I,J]:=S;A[IP1,J]:=T    
                "END";      
                AUX[3]:=AUX[3]+1;AUX[4]:=AUX[4]+1;I:=IP2; 
                ONEBYONE:="FALSE"     
             "END"
          "END"   
       "END";     
       "IF" ONEBYONE "THEN" 
       "BEGIN" "IF" TOL<ABS(A[I,I]) "THEN"      
          "BEGIN" AII:=A[I,I];DETAUX[I]:=A[I,I];
             "IF" AII>0 "THEN" AUX[3]:=AUX[3]+1 "ELSE" AUX[4]:=AUX[4]+1;      
             "FOR" J:=IP1 "STEP" 1 "UNTIL" N "DO" 
            "BEGIN" S:=-A[I,J]/AII;ELMROW(J,N,J,I,A,A,S);A[I,J]:=S "END"      
          "END";I:=IP1      
       "END"      
    "END" WHILE I;
    "IF" I=N "THEN" 
    "BEGIN" "IF" TOL<ABS(A[N,N]) "THEN" 
       "BEGIN" "IF" A[N,N]>0 "THEN" AUX[3]:=AUX[3]+1      
               "ELSE" AUX[4]:=AUX[4]+1
       "END";DETAUX[N]:=A[N,N]
    "END";
   ENDDEC:
    AUX[5]:=N-AUX[3]-AUX[4] 
 "END" DECSYM2;   
         "EOP"    
        
        
1SECTION : 3.1.1.1.1.3.2     ( DECEMBER 1978 )                  PAGE 1
        
        
        
 AUTHORS : J.R.BUNCH,L.KAUFMAN,B.N.PARLETT.     
        
        
 CONTRIBUTOR : C.H.CONVENT. 
        
        
 INSTITUTE : UNIVERSITY OF AMSTERDAM. 
        
        
 RECEIVED : 770712. 
        
        
 BRIEF DESCRIPTION :
        
    DETERMSYM2 CALCULATES THE DETERMINANT OF A SYMMETRIC MATRIX.    
    THE LDL' DECOMPOSITION OF THE MATRIX, AS PRODUCED BY DECSYM2,   
    SHOULD BE AVAILABLE.    
        
        
 KEYWORDS :       
        
    GENERAL SYMMETRIC MATRIX, 
    LDL' DECOMPOSITION,     
    BLOCK DIAGONAL PIVOTING;
        
        
 CALLING SEQUENCE : 
        
     THE HEADING OF THE PROCEDURE READS :       
     "REAL" "PROCEDURE" DETERMSYM2(DETAUX,N,AUX); 
     "VALUE" N;"INTEGER" N; 
     "ARRAY" DETAUX;"INTEGER" "ARRAY" AUX;      
     "CODE" 34294;
        
     DETERMSYM2 : DELIVERS THE CALCULATED VALUE OF THE DETERMINANT OF 
                  THE MATRIX; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     DETAUX :    <ARRAY IDENTIFIER>;  
                 "ARRAY" DETAUX[1:N]; 
                 ENTRY : THE ARRAY DETAUX AS PRODUCED BY DECSYM2;   
     N      :    <ARITHMETIC EXPRESSION>;       
                 THE ORDER OF THE ARRAY DETAUX  
                     ( = THE ORDER OF THE MATRIX );       
     AUX    :    <ARRAY IDENTIFIER>;  
                 "INTEGER" "ARRAY" AUX[2:5];    
                 ENTRY : THE ARRAY AUX AS PRODUCED BY DECSYM2;      
        
        
 PROCEDURES USED : NONE.    
        
1SECTION : 3.1.1.1.1.3.2     ( DECEMBER 1978 )                  PAGE 2
        
        
        
 RUNNING TIME : PROPORTIONAL TO N;    
        
        
        
 METHOD AND PERFORMANCE :   
        
    FIRST OF ALL DECSYM2 SHOULD BE CALLED TO PERFORM THE LDL'       
    DECOMPOSITION OF THE SYMMETRIC MATRIX,ACCORDING TO A METHOD DUE 
    TO BUNCH,KAUFMAN AND PARLETT (SEE [1],[2]). 
    IF A 1*1 BLOCK HAS BEEN COMPUTED FOR D THEN DETAUX[I] CONTAINS  
    THE VALUE OF D[I]. IF A 2*2 BLOCK HAS BEEN COMPUTED FOR D THEN  
    DETAUX[I]=1 AND DETAUX[I+1] CONTAIN THE VALUE OF THE DETERMINANT
    OF THE CORRESPONDING 2*2 BLOCK.   
    THE COMPUTATION OF THE DETERMINANT IS DONE BY CALCULATING THE   
    PRODUCT OF THE ELEMENTS OF DETAUX.
        
        
 REFERENCES :     
        
    1) J.R.BUNCH,L.KAUFMAN. 
       SOME STABLE METHODS FOR CALCULATING INERTIA AND SOLVING      
       SYMMETRIC LINEAR SYSTEMS.      
       MATHEMATICS OF COMPUTATION 31,P 163-180,1977.      
    2) J.R.BUNCH,L.KAUFMAN,B.N.PARLETT. 
       DECOMPOSITION OF A SYMMETRIC MATRIX.     
       NUMERISCHE MATHEMATIK 27,P 95-109,1976.  
        
        
        
        
 EXAMPLE OF USE  :
        
        
 "BEGIN" "COMMENT" EXAMPLE OF USE OF THE PROCEDURE DETERMSYM2;      
         "INTEGER" I,J;"REAL" TOL,DETERMINANT;  
         "REAL" "ARRAY" A[1:5,1:5],DETAUX[1:5]; 
         "INTEGER" "ARRAY" AUX[2:5],P[1:5];     
         "PROCEDURE" DECSYM2(A,N,TOL,AUX,P,DETAUX);"CODE" 34291;    
         "REAL" "PROCEDURE" DETERMSYM2(DETAUX,N,AUX);"CODE" 34294;  
        
         A[1,1]:=A[1,2]:=-3;A[1,3]:=-18;A[1,4]:=-30;A[1,5]:=18;     
         A[2,2]:=-1;A[2,3]:=-4;A[2,4]:=-48;A[2,5]:=8;     
         A[3,3]:=-6;A[3,4]:=-274;A[3,5]:=6;     
         A[4,4]:=119;A[4,5]:=19;      
         A[5,5]:=216;       
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         "FOR" J:=I+1 "STEP" 1 "UNTIL" 5 "DO" A[J,I]:=A[I,J];       
                                                               "COMMENT"      
1SECTION : 3.1.1.1.1.3.2     ( DECEMBER 1978 )                  PAGE 3
                                                                  ; 
        
        
         OUTPUT(61,"(""("THE COEFFICIENT MATRIX :")",/")"); 
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" 5 "DO"       
                 OUTPUT(61,"("-2ZD,4B")",A[I,J]); 
                 OUTPUT(61,"("/")")   
         "END";   
         TOL:="-14; 
         DECSYM2(A,5,TOL,AUX,P,DETAUX); 
         "IF" AUX[2]=1 "THEN" 
              OUTPUT(61,"("/,"("THE MATRIX IS SYMMETRIC")"")")      
         "ELSE" OUTPUT(61,"("/,"("THE MATRIX IS ASYMMETRIC.THE ")", 
                            "("RESULTS ARE MEANINGLESS")"")");      
         DETERMINANT:=DETERMSYM2(DETAUX,5,AUX); 
         OUTPUT(61,"("/,"("THE DETERMINANT OF THE MATRIX : ")",     
                     3Z3D.2D")",DETERMINANT)    
 "END"; 
        
        
 THIS DELIVERS AS RESULT :  
        
        
    THE COEFFICIENT MATRIX :
      -3      -3     -18    -30     18
      -3      -1      -4    -48      8
     -18      -4      -6   -274      6
     -30     -48    -274    119     19
      18       8       6     19    216
        
    THE MATRIX IS SYMMETRIC.
    THE DETERMINANT OF THE MATRIX :    168.00   
        
        
        
        
 SOURCE TEXT :    
        
 "CODE" 34294;    
 "REAL" "PROCEDURE" DETERMSYM2(DETAUX,N,AUX);   
 "VALUE" N;"INTEGER" N;     
 "ARRAY" DETAUX;"INTEGER" "ARRAY" AUX;
 "BEGIN" "INTEGER" I;"REAL" DET;      
         "IF" AUX[5]>0 "THEN" DET:=0 "ELSE"     
         "BEGIN" DET:=1;    
                 "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" DET:=DET*DETAUX[I]
         "END";   
         DETERMSYM2:=DET    
 "END" DETERMSYM2;
        
        
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 1
        
        
        
 AUTHORS : J.R.BUNCH,L.KAUFMAN,B.N.PARLETT.     
        
        
 CONTRIBUTOR : C.H.CONVENT. 
        
        
 INSTITUTE : UNIVERSITY OF AMSTERDAM. 
        
        
 RECEIVED : 770712. 
        
        
 BRIEF DESCRIPTION :
        
 THIS SECTION CONTAINS TWO PROCEDURES : 
 A) SOLSYM2 SOLVES A SYMMETRIC SYSTEM OF LINEAR EQUATIONS, ASSUMING 
    THAT THE MATRIX HAS BEEN DECOMPOSED INTO LDL' FORM BY A CALL OF 
    DECSYM2;      
 B) DECSOLSYM2 CALCULATES THE LDL' DECOMPOSITION OF A SYMMETRIC     
    MATRIX; MOREOVER, IF THIS MATRIX IS NON-SINGULAR, THEN IT SOLVES A
    CORRESPONDING SYSTEM OF LINEAR EQUATIONS;   
        
        
 KEYWORDS :       
        
    GENERAL SYMMETRIC MATRIX, 
    LDL' DECOMPOSITION,     
    BLOCK DIAGONAL PIVOTING;
        
        
 SUBSECTION : SOLSYM2.      
        
 CALLING SEQUENCE : 
        
    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" SOLSYM2(A,N,B,P,DETAUX);
    "VALUE" N;"INTEGER" N;  
    "ARRAY" A,B,DETAUX;"INTEGER" "ARRAY" P;     
    "CODE" 34292; 
        
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 2
        
        
        
    THE MEANING OF THE FORMAL PARAMETERS IS :   
    A :     <ARRAY IDENTIFIER>;       
            "ARRAY" A[1:N,1:N];       
            ENTRY : THE LDL' DECOMPOSITION OF A AS PRODUCED BY      
                    DECSYM2;
    N :     <ARITHMETIC EXPRESSION>;  
            THE ORDER OF THE MATRIX;  
    B :     <ARRAY IDENTIFIER>;       
            "ARRAY" B[1:N]; 
            ENTRY : THE RIGHT-HANDSIDE OF A SYSTEM OF LINEAR
                    EQUATIONS;
            EXIT  : THE CALCULATED SOLUTION VECTOR;       
    P :     <ARRAY IDENTIFIER>;       
            "INTEGER" "ARRAY" P[1:N]; 
            ENTRY : A VECTOR RECORDING THE INTERCHANGES PERFORMED ON
                    THE MATRIX BY THE PROCEDURE DECSYM2. P ALSO CONTAINS      
                    INFORMATION ON THE BLOCKSTRUCTURE OF THE MATRIX AS
                    DECOMPOSED BY DECSYM2;      
    DETAUX : <ARRAY IDENTIFIER>;      
             "ARRAY" DETAUX[1:N];     
             ENTRY : THE ARRAY DETAUX AS PRODUCED BY DECSYM2;       
        
        
 PROCEDURES USED :
        
    MATVEC=CP34011. 
    ELMVECROW=CP34026.      
        
        
 RUNNING TIME : ROUGHLY PROPORTIONAL TO N**2.   
        
        
 METHOD AND PERFORMANCE :   
        
    THE PROCEDURE SOLSYM2 COMPUTES THE SOLUTION OF A SYMMETRIC      
    SYSTEM OF LINEAR EQUATIONS,ASSUMING THAT THE MATRIX HAS BEEN    
    DECOMPOSED INTO LDL' FORM BY A CALL OF DECSYM2. B IS OVERWRITTEN
    WITH THE SOLUTION VECTOR. 
        
        
 REFERENCES :     
        
    1) J.R.BUNCH,L.KAUFMAN. 
       SOME STABLE METHODS FOR CALCULATING INERTIA AND SOLVING      
       SYMMETRIC LINEAR SYSTEMS.      
       MATHEMATICS OF COMPUTATION 31,P 163-180,1977.      
    2) J.R.BUNCH,L.KAUFMAN,B.N.PARLETT. 
       DECOMPOSITION OF A SYMMETRIC MATRIX.     
       NUMERISCHE MATHEMATIK 27,P 95-109,1976.  
        
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 3
        
        
        
 SUBSECTION : DECSOLSYM2.   
        
 CALLING SEQUENCE : 
        
    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" DECSOLSYM2(A,N,B,TOL,AUX);      
    "VALUE" N;"INTEGER" N;"REAL" TOL; 
    "ARRAY" A,B;"INTEGER" "ARRAY" AUX;
    "CODE" 34293; 
        
    THE MEANING OF THE FORMAL PARAMETERS IS :   
    A :     <ARRAY IDENTIFIER>;       
            "ARRAY" A[1:N,1:N];       
            ENTRY : SEE DECSYM2;      
            EXIT  : SEE DECSYM2;      
    N :     <ARITHMETIC EXPRESSION>;  
            THE ORDER OF THE MATRIX;  
    B :     <ARRAY IDENTIFIER>;       
            "ARRAY" B[1:N]; 
            ENTRY : SEE SOLSYM2;      
            EXIT  : THE CALCULATED SOLUTION VECTOR,WHEN A WAS FOUND 
                    TO BE NON-SINGULAR. 
                    B IS LEFT UNDISTURBED OTHERWISE;      
    TOL :   <ARITHMETIC EXPRESSION>;  
            ENTRY : SEE DECSYM2;      
    AUX :   <ARRAY IDENTIFIER>;       
            "INTEGER" "ARRAY" AUX[2:5]; 
            EXIT  : SEE DECSYM2;      
        
 PROCEDURES USED :
        
    DECSYM2=CP34291.
    SOLSYM2=CP34292.
        
        
 RUNNING TIME : ROUGHLY PROPORTIONAL TO N**3.   
        
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 4
        
        
        
 METHOD AND PERFORMANCE :   
        
    THE PROCEDURE DECSOLSYM2 COMPUTES THE SOLUTION OF A SYMMETRIC   
    SYSTEM OF LINEAR EQUATIONS. IT DOES SO BY FIRST CALLING THE     
    PROCEDURE DECSYM2 TO COMPUTE THE LDL' DECOMPOSITION OF THE      
    SYMMETRIC MATRIX,ACCORDING TO A METHOD DUE TO BUNCH,KAUFMAN AND 
    PARLETT (SEE [1],[2]).  
    WHEN THE MATRIX IS FOUND TO BE NON-SINGULAR THE PROCEDURE SOLSYM2 
    IS CALLED TO COMPUTE THE SOLUTION VECTOR,AND THE LATTER OVERWRITES
    B. WHEN THE MATRIX IS FOUND TO BE SINGULAR THE PROCEDURE SOLSYM2
    IS NOT CALLED AND B IS LEFT UNDISTURBED.    
        
        
 REFERENCES :     
        
    1) J.R.BUNCH,L.KAUFMAN. 
       SOME STABLE METHODS FOR CALCULATING INERTIA AND SOLVING      
       SYMMETRIC LINEAR SYSTEMS.      
       MATHEMATICS OF COMPUTATION 31,P 163-180,1977.      
    2) J.R.BUNCH,L.KAUFMAN,B.N.PARLETT. 
       DECOMPOSITION OF A SYMMETRIC MATRIX.     
       NUMERISCHE MATHEMATIK 27,P 95-109,1976.  
        
        
        
        
 EXAMPLES OF USE :
        
        
 "BEGIN" "COMMENT" EXAMPLE OF USE OF THE PROCEDURE DECSOLSYM2;      
         "INTEGER" I,J;"REAL" TOL;    
         "REAL" "ARRAY" A[1:5,1:5],B[1:5];      
         "INTEGER" "ARRAY" AUX[2:5];  
         "PROCEDURE" DECSOLSYM2(A,N,B,TOL,AUX);"CODE" 34293;
        
         A[1,1]:=A[1,2]:=-3;A[1,3]:=-18;A[1,4]:=-30;A[1,5]:=18;     
         A[2,2]:=-1;A[2,3]:=-4;A[2,4]:=-48;A[2,5]:=8;     
         A[3,3]:=-6;A[3,4]:=-274;A[3,5]:=6;     
         A[4,4]:=119;A[4,5]:=19;      
         A[5,5]:=216;       
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         "FOR" J:=I+1 "STEP" 1 "UNTIL" 5 "DO" A[J,I]:=A[I,J];       
        
         OUTPUT(61,"(""("THE COEFFICIENTMATRIX :")",/")");
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" 5 "DO"       
                 OUTPUT(61,"("-2ZD,4B")",A[I,J]); 
                 OUTPUT(61,"("/")")   
         "END";                                                "COMMENT"      
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 5
                                                                  ; 
        
        
         OUTPUT(61,"("/,"("THE RHS-VECTOR :")",/")");     
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         "BEGIN" INPUT(60,"("")",B[I]); 
                 OUTPUT(61,"("-Z3D,4B")",B[I])  
         "END";   
         TOL:="-14; 
         DECSOLSYM2(A,5,B,TOL,AUX);   
         OUTPUT(61,"("2/")"); 
         "IF" AUX[2]=1 "THEN" 
              OUTPUT(61,"(""("THE MATRIX IS SYMMETRIC")",/")")      
         "ELSE" OUTPUT(61,"(""("THE MATRIX IS ASYMMETRIC.THE ")",   
                            "("RESULTS ARE MEANINGLESS")",/")");    
         OUTPUT(61,"("/,"("INERTIA : <")",D,"(",")",D,"(",")",D,    
                     "(">")",/")",AUX[3],AUX[4],AUX[5]);  
         OUTPUT(61,"("/,"("THE COMPUTED SOLUTION :")",/")");
         "FOR" I:=1 "STEP" 1 "UNTIL" 5 "DO"     
         OUTPUT(61,"("-D.5D,4B")",B[I]) 
 "END"; 
        
        
 THIS DELIVERS AS RESULT :  
        
        
    THE COEFFICIENT MATRIX :
      -3      -3     -18    -30      18 
      -3      -1      -4    -48       8 
     -18      -4      -6   -274       6 
     -30     -48    -274    119      19 
      18       8       6     19     216 
        
    THE RHS-VECTOR :
     327     291   1290     275    1720 
        
    THE MATRIX IS SYMMETRIC.
    INERTIA : <3,2,0>       
        
    THE COMPUTED SOLUTION : 
    -7.00000   -2.00000   -1.00000    -4.00000     9.00000
1SECTION : 3.1.1.1.1.3.3     ( DECEMBER 1978 )                  PAGE 6
        
        
        
 SOURCE TEXT(S) : 
        
 "CODE" 34292;    
 "PROCEDURE" SOLSYM2(A,N,B,P,DETAUX); 
 "VALUE" N;"INTEGER" N;     
 "ARRAY" A,B,DETAUX;"INTEGER" "ARRAY" P;
 "BEGIN" "INTEGER" I,II,J,K,IP1,PI,PII,DUMMY;   
    "REAL" DET,TEMP,SAVE;   
    "REAL" "PROCEDURE" MATVEC(L,U,I,A,B);"CODE" 34011;    
    "PROCEDURE" ELMVECROW(L,U,I,A,B,X);"CODE" 34026;      
    I:=1; 
    "FOR" DUMMY:=0 "WHILE" I<N "DO"   
    "BEGIN" IP1:=I+1;PI:=P[I];SAVE:=B[PI];      
       "IF" P[IP1]>0 "THEN" 
       "BEGIN" B[PI]:=B[I];B[I]:=SAVE/A[I,I];   
           ELMVECROW(IP1,N,I,B,A,SAVE);I:=IP1   
       "END"      
       "ELSE"     
       "BEGIN" TEMP:=B[I];B[PI]:=B[IP1];DET:=DETAUX[IP1]; 
           B[I]:=(TEMP*A[IP1,IP1]-SAVE*A[I,IP1])/DET;     
           B[IP1]:=(SAVE*A[I,I]-TEMP*A[I,IP1])/DET;       
           ELMVECROW(I+2,N,I,B,A,TEMP);ELMVECROW(I+2,N,IP1,B,A,SAVE); 
           I:=I+2 
       "END"      
    "END" WHILE I;
    "IF" I=N "THEN" "BEGIN" B[I]:=B[I]/A[I,I];I:=N-1 "END"
       "ELSE" I:=N-2;       
    "FOR" DUMMY:=0 "WHILE" I>0 "DO"   
    "BEGIN" "IF" P[I]=0 "THEN" II:=I-1 "ELSE" II:=I;      
       "FOR" K:=II "STEP" 1 "UNTIL" I "DO"      
       "BEGIN" SAVE:=B[K];SAVE:=SAVE+MATVEC(I+1,N,K,A,B); 
          B[K]:=SAVE
       "END";     
       PII:=P[II];B[I]:=B[PII];B[PII]:=SAVE;I:=II-1       
    "END" WHILE I 
 "END" SOLSYM2;   
         "EOP"    
        
        
 "CODE" 34293;    
 "PROCEDURE" DECSOLSYM2(A,N,B,TOL,AUX); 
 "VALUE" N;"INTEGER" N;"REAL" TOL;    
 "ARRAY" A,B;"INTEGER" "ARRAY" AUX;   
 "BEGIN" "REAL" "ARRAY" DETAUX[1:N];"INTEGER" "ARRAY" P[1:N];       
         "PROCEDURE" DECSYM2(A,N,TOL,AUX,P,DETAUX);"CODE" 34291;    
         "PROCEDURE" SOLSYM2(A,N,B,P,DETAUX);"CODE" 34292;
         DECSYM2(A,N,TOL,AUX,P,DETAUX); 
         "IF" AUX[5]=0 "THEN" SOLSYM2(A,N,B,P,DETAUX)     
 "END" DECSOLSYM2;
         "EOP"    
        
        
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 1      
        
        
        
 AUTHORS: A.BJOERCK AND G.H.GOLUB.    
        
        
 CONTRIBUTOR: J.KOOPMAN.    
        
        
 INSTITUTE: UNIVERSITY OF AMSTERDAM.  
        
        
 RECEIVED: 780701 . 
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS TWO PROCEDURES FOR SOLVING A LINEAR LEAST
     SQUARES PROBLEM SUBJECT TO LINEAR CONSTRAINTS:       
     LSQDECOMP , FOR THE QR-DECOMPOSITION OF A LEAST SQUARES MATRIX,
                 WHERE  THIS MATRIX ALSO CONTAINS THE COEFFICIENTS OF 
                 THE LINEAR CONSTRAINTS (LINEAR EQUATIONS); 
     LSQREFSOL , FOR THE SOLUTION OF THIS LEAST SQUARES PROBLEM, IF 
                 THE MATRIX HAS BEEN DECOMPOSED BY LSQDECOMP.       
        
        
 KEYWORDS:
        
     LEAST SQUARES, 
     LINEAR CONSTRAINTS,    
     HOUSEHOLDER TRIANGULARIZATION,   
     ITERATIVE REFINEMENT.  
        
        
 SUBSECTION: LSQDECOMP.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LSQDECOMP(A, N, M, N1, AUX, AID, CI);    
     "VALUE"N,M,N1;"INTEGER"N,M,N1;"ARRAY" A,AUX,AID;     
     "INTEGER""ARRAY" CI;   
     "CODE" 34137;
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 2      
        
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
        
     A     :<ARRAY IDENTIFIER>;       
            "ARRAY" A[1:N,1:M];       
            ENTRY: THE ORIGINAL LEAST SQUARES MATRIX, WHERE THE FIRST 
                   N1 ROWS SHOULD FORM THE CONSTRAINT MATRIX (I.E. THE
                   FIRST N1 EQUATIONS ARE TO BE STRICTLY SATISFIED);
            EXIT : IN THE UPPER TRIANGLE OF A (THE ELEMENTS A[I,J] WITH       
                   I<J) THE SUPERDIAGONAL PART OF THE UPPER TRIANGULAR
                   MATRIX, PRODUCED BY HOUSEHOLDER TRANSFORMATIONS; IN
                   THE OTHER PART OF THE COLUMNS OF A THE SIGNIFICANT 
                   ELEMENTS OF THE GENERATING VECTORS OF THE HOUSEHOLDER      
                   MATRICES USED FOR THE TRIANGULARIZATION; 
     N     :<ARITHMETIC EXPRESSION>;  
            NUMBER OF ROWS OF THE MATRIX;       
     M     :<ARITHMETIC EXPRESSION>;  
            NUMBER OF COLUMNS OF THE MATRIX;    
     N1    :<ARITHMETIC EXPRESSION>;  
            NUMBER OF LINEAR CONSTRAINTS, I.E. THE FIRST N1 ROWS OF A 
            SET UP A SYSTEM OF N1 LINEAR EQUATIONS THAT MUST BE     
            STRICTLY SATISFIED (OF COURSE, IF THERE ARE NO CON-     
            STRAINTS, N1 MUST BE CHOSEN ZERO);  
     AUX   :<ARRAY IDENTIFIER>;       
            "ARRAY" AUX[2:7]; 
            ENTRY: AUX[2] CONTAINS A RELATIVE TOLERANCE FOR CALCULATING       
                   THE DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX;      
            EXIT:  AUX[3] CONTAINS THE NUMBER OF DIAGONAL ELEMENTS WHICH      
                   ARE NOT NEGLIGIBLE, NORMAL EXIT AUX[3]=M;
     AID   :<ARRAY IDENTIFIER>;       
            "ARRAY" AID[1:M]; 
            NORMAL EXIT (AUX[3]=M): THE DIAGONAL ELEMENTS OF THE UPPER
            TRIANGULAR MATRIX PRODUCED BY THE HOUSEHOLDER TRANSFORMATION      
     CI    :<ARRAY IDENTIFIER>;       
            "INTEGER""ARRAY" CI[1:M]; 
            EXIT: THE PIVOTAL INDICES OF THE INTERCHANGES OF THE COLUMNS      
            OF THE GIVEN MATRIX;      
        
        
 PROCEDURES USED: 
        
     MATMAT = CP34013.      
     TAMMAT = CP34014.      
     ELMCOL = CP34023.      
     ICHCOL = CP34031.      
        
        
 RUNNING TIME:    
        
     ROUGHLY PROPORTIONAL TO N*M**2-M**3/3.     
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 3      
        
        
        
 METHOD AND PERFORMANCE:    
        
     LET A DENOTE THE GIVEN MATRIX. LSQDECOMP PRODUCES AN N-TH ORDER
     ORTHOGONAL MATRIX Q AND AN N*M UPPER TRIANGULAR MATRIX R SUCH THAT       
     R EQUALS QA WITH PERMUTED COLUMNS. 
     THE ORTHOGONAL MATRIX Q IS FORMED AS A PRODUCT OF AT MOST M TRANS-       
     FORMATIONS OF THE FORM (I-BETA U U'). THESE HOUSEHOLDER MATRICES 
     REDUCE A TO THE MATRIX R: AT THE K-TH STAGE THE K-TH COLUMN OF THE       
     (ALREADY MODIFIED) MATRIX IS INTERCHANGED WITH THE COLUMN OF   
     MAXIMUM EUCLIDEAN NORM.THESE INTERCHANGES ARE RECORDED IN THE ARRAY      
     CI. PREMATURE TERMINATION OCCURS IF AT SOME STAGE THE EUCLIDEAN
     NORM OF THE PIVOTAL COLUMN IS LESS THAN SOME TOLERANCE (AUX[2])
     TIMES THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE COLUMNS OF THE 
     MATRIX.      
     LSQDECOMP DELIVERS THE UPPER TRIANGULAR MATRIX, WHERE THE DIAGONAL       
     ELEMENTS ARE STORED IN THE ARRAY AID AND THE OFF-DIAGONAL ELEMENTS       
     IN THE UPPER TRIANGULAR PART OF A. THE SIGNIFICANT ELEMENTS OF THE       
     GENERATING VECTORS OF THE HOUSEHOLDER TRANSFORMATIONS ARE STORED 
     IN THE COLUMNS OF A, I.E. ON AND BELOW THE LEADING DIAGONAL OF 
     A. IT SHOULD BE NOTED THAT FOR THE SOLUTION OF LEAST SQUARES   
     PROBLEMS, ONLY CALLS WITH N>=M ARE USEFUL. 
        
        
 REFERENCES:      
        
     A.BJOERCK AND G.H.GOLUB: 
     ITERATIVE REFINEMENT OF LEAST SQUARES SOLUTIONS BY HOUSEHOLDER 
     TRANSFORMATION,  BIT 7 (1967), PP. 322-337 
        
        
 SUBSECTION: LSQREFSOL.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LSQREFSOL(A, QR, N, M, N1, AUX,AID,CI,B,LDX,X,RES);
     "VALUE"N,M,N1;"INTEGER"N,M,N1;"INTEGER""ARRAY"CI;"REAL" LDX;   
     "ARRAY"A,QR,AUX,AID,B,X,RES;     
     "CODE" 34138;
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 4      
        
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A     :<ARRAY IDENTIFIER>;       
            "ARRAY" A[1:N,1:M];       
            THE ORIGINAL LEAST SQUARES MATRIX, WHERE THE FIRST N1   
            ROWS SET UP A SYSTEM OF LINEAR EQUATIONS THAT MUST BE   
            STRICTLY SATISFIED;       
     QR    :<ARRAY IDENTIFIER>;       
            "ARRAY" QR[1:N,1:M];      
            THE QR-DECOMPOSITION OF THE ORIGINAL LEAST SQUARES      
            MATRIX AS DELIVERED BY A SUCCESSFUL CALL OF LSQDECOMP;  
     N     :<ARITHMETIC EXPRESSION>;  
            NUMBER OF ROWS OF THE MATRICES A AND QR;      
     M     :<ARITHMETIC EXPRESSION>;  
            NUMBER OF COLUMNS OF THE MATRICES A AND QR;   
     N1    :<ARITHMETIC EXPRESSION>;  
            NUMBER OF LINEAR CONSTRAINTS;       
     AUX   :<ARRAY IDENTIFIER>;       
            "ARRAY" AUX[2:7]; 
            ENTRY: AUX[2] CONTAINS A RELATIVE TOLERANCE AS A
                   CRITERION TO STOP ITERATIVE REFINING: IF THE     
                   EUCLIDEAN NORM OF THE CORRECTION IS SMALLER      
                   THAN AUX[2] TIMES THE CURRENT APPROXIMATION OF   
                   THE SOLUTION, THE ITERATIVE REFINING IS STOPPED; 
                   AUX[6]: MAXIMUM NUMBER OF ITERATIONS ALLOWED     
                   (USUALLY AUX[6]=5 WILL BE SUFFICIENT); 
            EXIT : AUX[7]: THE NUMBER OF ITERATIONS PERFORMED;      
    AID    :<ARRAY IDENTIFIER>;       
            "ARRAY" AID[1:M]; 
            THE DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX    
            AS DELIVERED BY A SUCCESSFUL CALL OF LSQDECOMP; 
     CI    :<ARRAY IDENTIFIER>;       
            "INTEGER""ARRAY" CI[1:M]; 
            THE PIVOTAL INDICES AS PRODUCED BY LSQDECOMP; 
     B     :<ARRAY IDENTIFIER>;       
            "ARRAY" B[1:N]; 
            THE RIGHT-HAND SIDE OF THE LEAST SQUARES PROBLEM;       
            FIRST N1 ELEMENTS FORM THE RIGHT HAND SIDES OF THE      
            CONSTRAINTS;    
     LDX   :<REAL VARIABLE>;
            THE EUCLIDEAN NORM OF THE LAST CORRECTION OF THE
            SOLUTION;       
     X     :<ARRAY IDENTIFIER>;       
            "ARRAY" X[1:M]; 
            EXIT: THE SOLUTION VECTOR;
     RES   :<ARRAY IDENTIFIER>;       
            "ARRAY" RES[1:N]; 
            EXIT: THE RESIDUAL VECTOR CORRESPONDING TO THE SOLUTION;
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 5      
        
        
        
 PROCEDURES USED: 
        
     VECVEC = CP34010.      
     MATVEC = CP34011.      
     TAMVEC = CP34012.      
     ELMVECCOL = CP34021.   
     ICHCOL = CP34031.      
     LNG SUB = CP31106.     
     LNGMATVEC = CP34411.   
     LNGTAMVEC = CP34412.   
        
        
 RUNNING TIME:    
        
     ROUGHLY PROPORTIONAL TO C1*M*N-C2*M**2, WHERE C1 AND C2
     ARE CONSTANTS. 
        
        
 METHOD AND PERFORMANCE:    
        
     LSQREFSOL SHOULD BE CALLED AFTER A SUCCESSFUL CALL OF LSQDECOMP
     (I.E. AUX[3]=M). LSQREFSOL YIELDS THE LEAST SQUARES SOLUTION OF
     THE OVERDETERMINED SYSTEM WITH THE DECOMPOSED COEFFICIENT MATRIX 
     IN THE ARRAY A AND THE RIGHT-HAND SIDE IN ARRAY B. THE ORIGINAL
     LEAST SQUARES MATRIX ALSO CONTAINS THE LINEAR CONSTRAINTS (THE 
     FIRST N1 ROWS OF THIS MATRIX SET UP A SYSTEM THAT MUST BE STRICT-
     LY SATISFIED). FIRST, THE ORTHOGONAL TRANSFORMATION WITH THE   
     HOUSEHOLDER MATRICES IS PERFORMED ON THE RIGHT-HAND SIDE. NEXT 
     THE SYSTEM IS SOLVED BY MEANS OF A BACK SUBSTITUTION. IN THIS  
     WAY THE FIRST APPROXIMATION OF THE SOLUTION (WITH PERMUTED     
     COLUMNS) IS OBTAINED.  
     AFTER THIS AN ITERATIVE PROCESS REFINES THE APPROXIMATION      
     UNTIL THE EUCLIDEAN NORM OF THE CORRECTION VECTOR IS NEGLIGIBLY
     SMALL COMPARED TO THE APPROXIMATION OR UNTIL THE MAXIMUM NUMBER
     OF ITERATIONS IS REACHED. FOR A MORE DETAILED DESCRIPTION OF   
     THE ITERATIVE IMPROVEMENT, SEE REF[1].     
     AFTER THE ITERATIVE PROCESS AN APPROXIMATION TO THE SOLUTION IS
     FOUND. HOWEVER, THE ORDER OF THE COMPONENTS POSSIBLY IS NOT    
     CORRECT. THEREFORE THIS ORDER IS RESTORED AT THE END OF THE    
     PROCEDURE (SEE ALSO METHOD AND PERFORMANCE OF LSQDECOMP).      
     THE LEAST SQUARES SOLUTIONS OF SEVERAL OVERDETERMINED SYSTEMS  
     WITH THE SAME CONSTRAINTS AND COEFFICIENT MATRIX CAN BE SOLVED 
     BY SUCCESSIVE CALLS OF LSQREFSOL WITH DIFFERENT RIGHT-HAND     
     SIDES.       
        
        
 REFERENCES:      
        
     [1] A.BJOERCK AND G.H.GOLUB:     
     ITERATIVE REFINEMENT OF LEAST SQUARES SOLUTIONS BY HOUSEHOLDER 
     TRANSFORMATION, BIT 7 (1967), PP. 322-337. 
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
        
     THE FOLLOWING PROGRAM SOLVES THE PROBLEM:  
        
     MINIMIZE  //B - A*X//  
        
     UNDER    X1 + 1000*X2 + 5*X3 = 2016 ,WHERE 
     A IS THE MATRIX
        
        
      1    0    8 
      0    3    2 
      1    2    "-5 
      0    0    0 
        
     AND THE VECTOR B = (25, 12, 5.00003, 1)' , X = (X1, X2, X3)' ; 
        
        
 "BEGIN" "INTEGER" N,M,N1,I,J;
         N := 5; M := 3; N1 := 1;     
         "BEGIN""INTEGER""ARRAY" CI[1:M];       
             "ARRAY" AUX[2:7],QR,A[1:N,1:M],B,RES[1:N],AID,X[1:M];  
             "REAL" LDX;    
             "PROCEDURE" LSQDECOMP(A,N,M,N1,AUX,AID,CI);"CODE"34137;
             "PROCEDURE" LSQREFSOL(A,QR,N,M,N1,AUX,AID,CI,B,LDX,X,RES);       
             "CODE" 34138;  
             A[1,1] := 1; A[1,2] := 1000; A[1,3] := 5;    
             A[2,1] := 1; A[2,2] := 0; A[2,3] := 8;       
             A[3,1] := 0; A[3,2] := 3; A[3,3] := 2;       
             A[4,1] := 1; A[4,2] := 2; A[4,3] := "-5;     
             A[5,1]:=A[5,2]:=A[5,3]:=0; 
             B[1] := 2016; B[2] :=25; B[3] :=12; B[4] := 5.00003;   
             B[5] := 1;  AUX[2] := "-14; AUX[6] := 5;     
             "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" 
             "FOR" J:=1 "STEP" 1 "UNTIL" M "DO" QR[I,J] := A[I,J];  
        
             LSQDECOMP(QR,N,M,N1,AUX,AID,CI);   
             LSQREFSOL(A,QR,N,M,N1,AUX,AID,CI,B,LDX,X,RES); 
        
 OUTPUT(61,"(""("THE SOLUTION VECTOR: ")",//")"); 
 "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" OUTPUT(61,"("/")",X[I]);
 OUTPUT(61,"("//,"("THE RESIDUAL VECTOR: ")",//")");      
 "FOR" J:=N1+1 "STEP" 1 "UNTIL" N "DO" OUTPUT(61,"("/")",RES[J]);   
 OUTPUT(61,"("///,"("NUMBER OF ITERATIONS: ")",D,//,      
        "("NORM LAST CORRECTION OF X: ")",N")",AUX[7],LDX)
        "END"     
"END"   
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 7      
        
        
        
 DELIVERS:
        
        
 THE SOLUTION VECTOR:       
        
 +1.0000000000000"+000      
 +2.0000000000000"+000      
 +3.0000000000000"+000      
        
        
 THE RESIDUAL VECTOR:       
        
 -5.2734444477081"-016      
 +2.1280091641666"-015      
 +5.3479806840033"-016      
 +1.0000000000000"+000      
        
 NUMBER OF ITERATIONS: 2    
        
 NORM LAST CORRECTION OF X:  +2.1657844626990"-015
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34137;    
 "PROCEDURE" LSQDECOMP( A, N ,M ,N1 ,AUX ,AID ,CI );      
 "VALUE" N , M ,N1;"INTEGER" N,M,N1;"ARRAY" A,AUX, AID;   
 "INTEGER""ARRAY" CI;       
 "BEGIN""INTEGER"J,K,KPIV,NR,S;"BOOLEAN" FSUM;  
     "REAL" BETA,SIGMA,NORM,AIDK,AKK,W,EPS;     
     "ARRAY" SUM[1:M];      
     "REAL""PROCEDURE" MATMAT(L ,U ,I ,J ,A ,B );"CODE" 34013;      
     "REAL""PROCEDURE" TAMMAT(L ,U ,I ,J ,A ,B );"CODE" 34014;      
     "PROCEDURE" ELMCOL(L ,U ,I ,J ,A ,B ,X );"CODE" 34023; 
     "PROCEDURE" ICHCOL(L ,U ,I ,J ,A );"CODE" 34031;     
        
     NORM:=0 ; AUX[3]:=M;NR:=N1;FSUM:="TRUE";   
     "FOR" K:=1 "STEP" 1 "UNTIL" M "DO" 
         "BEGIN" "IF" K=N1+1 "THEN" "BEGIN" FSUM:="TRUE"; NR:=N "END";
         "IF" FSUM "THEN"   
         "FOR" J:=K "STEP" 1 "UNTIL" M "DO"     
         "BEGIN" W:=SUM[J]:= TAMMAT(K ,NR ,J ,J ,A ,A);   
                 "IF" W>NORM "THEN" NORM:=W     
         "END"; FSUM:="FALSE";EPS:=AUX[2]*SQRT(NORM);     
         SIGMA:=SUM[K]; KPIV:=K;      
         "FOR" J:=K+1 "STEP" 1 "UNTIL" M "DO"   
         "IF" SUM[J]>SIGMA "THEN"     
         "BEGIN" SIGMA:=SUM[J]; KPIV:=J "END";  
                                                               "COMMENT"      
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 8      
                                                                  ; 
        
        
         "IF" KPIV^=K "THEN"
         "BEGIN" SUM[KPIV]:=SUM[K]; ICHCOL( 1 ,N ,K ,KPIV ,A) "END";
         CI[K]:=KPIV; AKK:=A[K,K];    
         SIGMA:=TAMMAT(K ,NR ,K ,K ,A ,A); W:=SQRT(SIGMA);
         AIDK:=AID[K]:="IF" AKK<0 "THEN" W "ELSE" -W;     
         "IF" W<EPS "THEN"  
         "BEGIN" AUX[3]:=K-1;"GOTO" ENDDEC "END"; 
         BETA:= 1/(SIGMA-AKK*AIDK); A[K,K]:=AKK-AIDK;     
         "FOR" J:=K+1 "STEP" 1 "UNTIL" M "DO"   
         "BEGIN" ELMCOL(K ,NR ,J ,K ,A ,A ,-BETA*TAMMAT(K ,NR ,     
         K ,J ,A ,A)); SUM[J]:=SUM[J]-A[K,J]**2 
         "END";   
         "IF" K=N1 "THEN"   
         "FOR" J:=N1+1 "STEP" 1 "UNTIL" N "DO"  
         "FOR" S:=1 "STEP" 1 "UNTIL" M "DO"     
         "BEGIN" NR:="IF" S>N1 "THEN" N1 "ELSE" S-1;      
         W:=A[J,S]-MATMAT(1 ,NR , J ,S ,A ,A);  
         A[J,S]:="IF" S>N1 "THEN" W "ELSE" W/AID[S]       
         "END"    
     "END" FOR K; 
   ENDDEC:
 "END" LSQDECOMP; 
         "EOP"    
 "CODE" 34138;    
 "PROCEDURE" LSQREFSOL(A, QR, N, M, N1, AUX, AID, CI, B,LDX,X,RES); 
 "VALUE"N,M,N1;"INTEGER"N,M,N1;"INTEGER""ARRAY"CI;"REAL"LDX;
 "ARRAY"QR,A,AID,AUX,X,RES,B; 
 "BEGIN""INTEGER"I,J,K,S;   
        "REAL"C1,NEXVE,NDX,NDR,D,DD,OP,OPL,CORRNORM;      
        "ARRAY"F[1:N],G[1:M]; 
        "REAL""PROCEDURE"VECVEC(L,U,S,A,B);"CODE"34010;   
        "REAL""PROCEDURE"MATVEC(L,U,S,A,B);"CODE"34011;   
        "REAL""PROCEDURE"TAMVEC(L,U,S,A,B);"CODE"34012;   
        "PROCEDURE" ICHCOL(L,U,I,J,A);"CODE" 34031;       
        "PROCEDURE"LNG SUB(A,AA,B,BB,C,CC);"CODE"31106;   
        "PROCEDURE"ELMVECCOL(L,U,I,A,B,X);"CODE"34021;    
        "PROCEDURE"LNGTAMVEC(L,U,I,A,B,C,CC,D,DD);"CODE"34412;      
        "PROCEDURE"LNGMATVEC(L,U,I,A,B,C,CC,D,DD);"CODE"34411;      
        "PROCEDURE" HOUSEHOLDER(P, Q, R, E);    
        "VALUE" P,Q,R,E;"INTEGER" P,Q,R,E;      
        "BEGIN" "FOR" S:=P "STEP" Q "UNTIL" R "DO"
                ELMVECCOL(S, E, S, F, QR,TAMVEC(S, E, S,QR, F)/(QR[S,S]*      
                AID[S]))    
        "END";    
        "FOR" J:=1 "STEP" 1 "UNTIL" M "DO"      
        "BEGIN" S:=CI[J];"IF" S^=J"THEN" ICHCOL(1,N,J,S,A) "END";   
        "FOR"J:=1"STEP"1"UNTIL"M"DO"X[J]:=G[J]:=0;
        "FOR"I:=1"STEP"1"UNTIL"N"DO"  
        "BEGIN"RES[I]:=0;F[I]:=B[I]"END";       
                                                               "COMMENT"      
1SECTION : 3.1.1.2.1.4       (DECEMBER 1978)                      PAGE 9      
                                                                  ; 
        
        
        "FOR"K:=0,1,K+1     
           "WHILE" (CORRNORM>AUX[2]*NEXVE & K<=AUX[6])    
        "DO"      
        "BEGIN"NDX:=NDR:=0; 
           "IF"K^=0"THEN"   
           "BEGIN""FOR"I:=1"STEP"1"UNTIL"N"DO"RES[I]:=RES[I]+F[I];  
                  "FOR"S:=1"STEP"1"UNTIL"M"DO"  
                  "BEGIN" X[S]:=X[S]+G[S];      
                         LNGTAMVEC(1,N,S,A,RES,0,0,D,DD); 
                         G[S]:=(-D-TAMVEC(1,S-1,S,QR,G))/AID[S]     
                  "END";    
                  "FOR"I:=1"STEP"1"UNTIL"N"DO"  
                  "BEGIN" LNGMATVEC(1, M, I, A, X,
                          "IF" I>N1 "THEN" RES[I] "ELSE" 0, 0, D, DD);
                         LNG SUB(B[I],0,D,DD,OP,OPL);     
                         F[I]:=OP     
                  "END"     
           "END"; 
           NEXVE:=SQRT(VECVEC(1,M,0,X,X)+VECVEC(1,N,0,RES,RES));    
           HOUSEHOLDER(1, 1, N1, N1); 
           "FOR" I:=N1+1 "STEP" 1 "UNTIL" N "DO"
           F[I]:=F[I]-MATVEC(1,N1, I, QR, F);   
           HOUSEHOLDER(N1+1, 1, M, N);
           "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"   
           "BEGIN" C1:=F[I];F[I]:=G[I]; 
               G[I]:="IF" I>N1 "THEN" C1-G[I] "ELSE" C1   
           "END"; 
           "FOR"S:=M"STEP"-1"UNTIL"1"DO"
           "BEGIN"G[S]:=(G[S]-MATVEC(S+1,M,S,QR,G))/AID[S]; 
                  NDX:=NDX+G[S]**2    
           "END"; 
           HOUSEHOLDER(M, -1, N1+1, N); 
           "FOR" S:=1 "STEP" 1 "UNTIL" N1 "DO"  
           F[S]:=F[S]-TAMVEC(N1+1, N, S, QR, F);
           HOUSEHOLDER(N1, -1, 1, N1);
           AUX[7]:=K;       
           "FOR"I:=1"STEP"1"UNTIL"N"DO"NDR:=NDR+F[I]**2;  
           CORRNORM:=SQRT(NDX+NDR)    
        "END"FOR K; 
        LDX:=SQRT(NDX);     
        "FOR" S:=M "STEP" -1 "UNTIL" 1 "DO"     
        "BEGIN" J:=CI[S];"IF" J^=S "THEN"       
           "BEGIN" C1:=X[J];X[J]:=X[S];X[S]:=C1;ICHCOL(1,N,J,S,A)"END"
        "END"     
 "END" LSQREFSOL; 
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 1
 
 
 
 AUTHORS:  TH.J. DEKKER AND TH.H.P. REYMER
 
 
 CONTRIBUTOR:  TH.H.P. REYMER
 
 
 INSTITUTE:  UNIVERSITY OF AMSTERDAM
 
 
 RECEIVED:  770427;
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS TWO PROCEDURES:
     A) ZERPOL   CALCULATES  ALL  ROOTS   OF  A   POLYNOMIAL  WITH  REAL
        COEFFICIENTS BY MEANS OF LAGUERRE'S METHOD;
     B) BOUNDS   CALCULATES  UPPERBOUNDS FOR THE ABSOLUTE ERROR IN GIVEN
        APPROXIMATED ZEROS OF A POLYNOMIAL WITH REAL COEFFICIENTS;
 
 
 KEYWORDS:
 
     ZEROS;
     REAL POLYNOMIAL;
     LAGUERRE'S METHOD;
     ERROR BOUNDS;
 
 
 SUBSECTION: ZERPOL
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "INTEGER" "PROCEDURE" ZERPOL(N, A, EM, RE, IM, D);
     "VALUE" N; "INTEGER" N; "ARRAY" A, EM, RE, IM, D;
     "CODE" 34501 ;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     N:      <ARITHMETIC EXPRESSION>;
             ENTRY:  THE DEGREE OF THE POLYNOMIAL;
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[0 : N];
             ENTRY:  THE COEFFICIENTS OF THE POLYNOMIAL, IN  SUCH  A WAY
                     THAT
                     P(Z) = (..(A[N]*Z + A[N-1])*Z +..+ A[1])*Z + A[0];
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 2
 
 
 
     EM:     <ARRAY IDENTIFIER>;
             "ARRAY" EM[0 : 4];
             ENTRY:  EM[0]: MACHINE PRECISION;
                     EM[1]: THE MAXIMAL NUMBER OF ITERATIONS ALLOWED FOR
                            EACH ZERO;
             EXIT:   EM[2]: FAIL INDICATION;
                            0  SUCCESSFUL CALL;
                            1  UPON ENTRY DEGREE N <= 0;
                            2  UPON ENTRY LEADING COEFFICIENT A[N] = 0;
                            3  NUMBER OF ITERATIONS EXCEEDED EM[1];
                     EM[3]: NUMBER OF NEW STARTS IN THE LAST ITERATION;
                     EM[4]: TOTAL NUMBER OF ITERATIONS PERFORMED;
             FOR THE CD CYBER 70 SYSTEM SUITABLE VALUES ARE:
             EM[0]:= "-14;
             EM[1]:= 40; IF, UPON EXIT, EM[2] = 3  AND EM[3] < 5 THEN IT
                         MAY BE  USEFUL TO  START  AGAIN  WITH  A HIGHER
                         VALUE OF EM[1];
     RE, IM: <ARRAY IDENTIFIERS>;
             "ARRAY" RE, IM[1 : N];
             EXIT:   THE  REAL  AND  IMAGINARY PARTS OF THE ZEROS OF THE
                     POLYNOMIAL; THE  MEMBERS  OF  EACH  NONREAL COMPLEX
                     CONJUGATE PAIR ARE CONSECUTIVE;
     D:      <ARRAY IDENTIFIER>;
             "ARRAY" D[0 : N];
             EXIT:   IF THE CALL IS UNSUCCESSFUL AND ONLY N-K ZEROS HAVE
                     BEEN   FOUND,  THEN   D[0 :  K]   CONTAINS   THE
                     COEFFICIENTS   OF   THE   (DEFLATED)    POLYNOMIAL;
                     MOREOVER,  THEN THE ZEROS  FOUND  ARE  DELIVERED IN
                     RE, IM[ K + 1 : N], WHEREAS THE REMAINING  PARTS
                     OF RE AND IM CONTAIN NO INFORMATION;
 
     ZERPOL:= THE NUMBER, K, OF ZEROS NOT FOUND;
 
 
 PROCEDURES USED:
 
     DWARF   = CP30003;
     GIANT   = CP30004;
     COMABS  = CP34340;
     COMSQRT = CP34343;
 
 
 REQUIRED CENTRAL MEMORY:
 
     TOTAL SIZE OF LOCAL ARRAYS IS N + 16 REAL LOCATIONS;
 
 
 RUNNING TIME:
 
     ROUGHLY PROPORTIONAL TO N**2;
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 3
 
 
 
 METHOD AND PERFORMANCE:
 
     THE  PROCEDURE  USES  LAGUERRE'S  METHOD TO FIND ZEROS OF THE GIVEN
     POLYNOMIAL (SEE [2]);  WHEN A ZERO HAS BEEN FOUND, A COMPOSITE
     DEFLATION TECHNIQUE IS USED TO  OBTAIN A  NEW POLYNOMIAL OF LOWER
     DEGREE (SEE [1], [3], [4]); IF CONVERGENCE IS NOT APPARENT, SEVERAL
     RESTARTS, THE NUMBER OF WHICH DEPENDS ON  EM[1]  BUT HAS A MAXIMUM
     OF  5, ARE MADE IN THE NEIGHBOURHOOD OF THE ABSOLUTE LARGEST ZERO;
     THE  ACCURACY OF THE  CALCULATED  ZEROS  STRONGLY  DEPENDS  ON  THE
     POLYNOMIAL; A ROUGH INDICATION FOR THE ERROR IN A CALCULATED ZERO Z
     FOLLOWS FROM P(Z) / DP(Z), WHERE P DENOTES THE GIVEN POLYNOMIAL AND
     DP ITS  FIRST DERIVATIVE (SEE E.G. [5]); TO FIND A  TRUE UPPERBOUND
     FOR THESE ERRORS, ONE CAN USE PROCEDURE BOUNDS (SEE NEXT
     SUBSECTION); FOR  A MORE  DETAILED DESCRIPTION OF THE PROCEDURE AND
     TEST RESULTS SEE [4];
 
 
 REFERENCES:
 
     [1] D.A. ADAMS, A STOPPING CRITERION FOR POLYNOMIAL ROOT FINDING,
         CACM 10, NO 10, PP. 655-658, OCTOBER 1967;
     [2] T.J. DEKKER, NEWTON-LAGUERRE ITERATION,
         MATHEMATISCH CENTRUM MR82, 1966;
     [3] G. PETERS AND J.H. WILKINSON, PRACTICAL PROBLEMS ARISING IN THE
         SOLUTION OF POLYNOMIAL EQUATIONS,
         J. INST. MATHS APPLICS 1971, NO. 8, PP. 16-35;
     [4] TH.H.P. REYMER, BEREKENING VAN NULPUNTEN VAN REELE POLYNOMEN EN
         FOUTGRENZEN VOOR DEZE NULPUNTEN,
         DOCTORAAL SCRIPTIE UVA, APRIL 1977;
     [5] J.H. WILKINSON, ROUNDING ERRORS IN ALGEBRAIC PROCESSES,
         PRENTICE HALL, 1963;
 
 
 EXAMPLE OF USE:
 
     FOR AN EXAMPLE OF USE SEE PROCEDURE BOUNDS (NEXT SUBSECTION);
 
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 4
 
 
 
 SUBSECTION: BOUNDS
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" BOUNDS(N,A,RE,IM,RELE,ABSE,RECENTRE,IMCENTRE,BOUND);
     "VALUE" N, RELE, ABSE; "INTEGER" N; "REAL" RELE, ABSE;
     "ARRAY" A, RE, IM, RECENTRE, IMCENTRE, BOUND;
     "CODE" 34502 ;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     N:        <ARITHMETIC EXPRESSION>;
               ENTRY:  DEGREE OF THE POLYNOMIAL;
     A:        <ARRAY IDENTIFIER>;
               "ARRAY" A[0 : N];
               ENTRY:  THE  COEFFICIENTS  OF  THE  POLYNOMIAL  OF  WHICH
                       RE[J] + I * IM[J]  ARE  THE  APPROXIMATED  ZEROS,
                       IN SUCH A WAY THAT
                       P(Z) = (..(A[N]*Z + A[N-1])*Z +..+A[1])*Z + A[0]
     RE, IM:   <ARRAY IDENTIFIERS>;
               "ARRAY" RE, IM[1 : N];
               ENTRY:  REAL AND IMAGINARY PARTS OF APPROXIMATED ZEROS OF
                       A  POLYNOMIAL,  SUCH  THAT  THE  MEMBERS  OF EACH
                       NONREAL COMPLEX CONJUGATE PAIR ARE CONSECUTIVE;
               EXIT:   A PERMUTATION OF THE INPUT DATA;
     RELE:     <ARITHMETIC EXPRESSION>;
               ENTRY:  RELATIVE ERROR IN THE  NON-VANISHING COEFFICIENTS
                       A[J] OF THE GIVEN POLYNOMIAL;
     ABSE:     <ARITHMETIC EXPRESSION>;
               ENTRY:  ABSOLUTE ERROR IN THE VANISHING COEFFICIENTS A[J]
                       GIVEN  POLYNOMIAL;  IF  THERE  ARE  NO  VANISHING
                       COEFFICIENTS, ABSE SHOULD BE ZERO;
     RECENTRE, IMCENTRE: <ARRAY IDENTIFIERS>;
               "ARRAY" RECENTRE, IMCENTRE[1 : N];
               EXIT:   REAL  AND IMAGINARY PARTS OF THE CENTERS OF DISKS
                       IN WHICH  SOME  NUMBER OF ZEROS OF THE POLYNOMIAL
                       GIVEN BY A ARE SITUATED;
                       THE NUMBER OF IDENTICAL CENTERS  DENOTES
                       THE NUMBER OF ZEROS IN THAT DISK;
     BOUND:    <ARRAY IDENTIFIER>;
               "ARRAY" BOUND[1 : N];
               EXIT:   RADIUS  OF THE  DISKS  WHOSE  CENTERS  ARE  GIVEN
                       CORRESPONDINGLY IN RECENTRE AND IMCENTRE;
 
 
 PROCEDURES USED:
 
     ARREB = CP30002;
     GIANT = CP30004;
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 5
 
 
 
 REQUIRED CENTRAL MEMORY:
 
     TOTAL SIZE OF LOCAL ARRAYS IS AT MOST 7 * N REAL LOCATIONS;
 
 
 RUNNING TIME:
 
     APPROXIMATELY OF ORDER N**2;
 
 METHOD AND PERFORMANCE:
 
     FROM  THE  APPROXIMATED  ZEROS  A  POLYNOMIAL IS  RECONSTRUCTED AND
     COMPARED  WITH  THE  GIVEN  POLYNOMIAL; SUBSEQUENTLY, THE PROCEDURE
     CALCULATES  DISKS  SUCH THAT THE NUMBER OF GIVEN APPROXIMATED ZEROS
     WITHIN EACH DISK EQUALS THE NUMBER OF ZEROS OF THE GIVEN POLYNOMIAL
     WITHIN  THAT  DISK;  UPON  EXIT  EVERY  TWO NON-IDENTICAL DISKS ARE
     DISJOINT;
     FOR A MORE DETAILED DESCRIPTION SEE [1], [2];
 
 
 REFERENCES:
 
     [1] G. PETERS AND J.H. WILKINSON, PRACTICAL PROBLEMS ARISING IN THE
         SOLUTION OF POLYNOMIAL EQUATIONS, J.INST.MATHS APPLICS 1971,
         NO 8, PP. 16-35;
     [2] TH.H.P. REYMER, BEREKENING VAN NULPUNTEN VAN REELE POLYNOMEN EN
         FOUTGRENZEN VOOR DEZE NULPUNTEN,
         DOCTORAAL SCRIPTIE UVA, APRIL 1977;
 
 
 EXAMPLE OF USE:
 
 "BEGIN" "INTEGER" I, J;
   "ARRAY" A, D[0:7], RE, IM[1:7], EM[0:4];
   "INTEGER""PROCEDURE" ZERPOL(N,A,EM,RE,IM,D); "CODE"34501;
 
   A[7]:= 1; A[6]:= -3; A[5]:= -3; A[4]:= 25; A[3]:= -46;
   A[2]:= 38; A[1]:= -12; A[0]:= 0;
   EM[0]:= "-14; EM[1]:= 40;
   I:= ZERPOL(7, A, EM, RE, IM, D);
   OUTPUT(61,"(""("COEFFICIENTS OF POLYNOMIAL:")",//")");
   "FOR" J:=7 "STEP" -1 "UNTIL" 0 "DO"
   OUTPUT(61,"("-2ZD3B")",A[J]);
   OUTPUT(61,"("//,"("NUMBER NOT FOUND ZEROS ")",3ZD,/,
      "("FAIL INDICATION ")",3ZD,/,"("NUMBER NEW STARTS ")",3ZD,/,
      "("NUMBER OF ITERATIONS ")",3ZD,/")",I,EM[2],EM[3],EM[4]);
   OUTPUT(61,"("/,"("ZEROS: ")",/")");
   "FOR" J:= I+1 "STEP" 1 "UNTIL" 7 "DO" "IF" IM[J] = 0
   "THEN" OUTPUT(61,"("/,N")",RE[J])
   "ELSE" OUTPUT(61,"("/,2(N)")",RE[J],IM[J]);
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 6
 
 
 
   "IF" I = 0 "THEN"
   "BEGIN" "ARRAY" RECENTRE, IMCENTRE, BOUND[1:7];
     "PROCEDURE" BOUNDS(N,A,RE,IM,RELE,ABSE,RECENTRE,IMCENTRE,BOUND);
     "CODE" 34502;
 
     BOUNDS(7, A, RE, IM, 0, 0, RECENTRE, IMCENTRE, BOUND);
     OUTPUT(61,"("2/,"("REAL AND IMAG. PART OF CENTRE + RADIUS")",/")");
     "FOR" J:= 1 "STEP" 1 "UNTIL" 7 "DO"
     OUTPUT(61,"("/,3(N)")",RECENTRE[J],IMCENTRE[J],BOUND[J])
   "END"
 "END"
 
 
 RESULTS :
 
 
 COEFFICIENTS OF POLYNOMIAL:
 
    1     -3     -3     25    -46     38    -12      0
 
 NUMBER NOT FOUND ZEROS    0
 FAIL INDICATION    0
 NUMBER NEW STARTS    0
 NUMBER OF ITERATIONS   11
 
 ZEROS:
 
 +2.0000000000000"+000
 -3.0000000000000"+000
 +1.0000000000000"+000  -1.0000000000000"+000
 +1.0000000000000"+000  +1.0000000000000"+000
 +1.0000000083024"+000
 +9.9999999169752"-001
 +0.0000000000000"+000
 
 REAL AND IMAG. PART OF CENTRE + RADIUS
 
 +2.0000000000000"+000  +0.0000000000000"+000  +1.3238111117716"-011
 -3.0000000000000"+000  +0.0000000000000"+000  +3.8857510604494"-013
 +1.0000000000000"+000  -1.0000000000000"+000  +4.0912729775463"-012
 +1.0000000000000"+000  +1.0000000000000"+000  +4.0912729775463"-012
 +9.9999999999998"-001  +0.0000000000000"+000  +2.2533888428865"-006
 +9.9999999999998"-001  +0.0000000000000"+000  +2.2533888428865"-006
 +0.0000000000000"+000  +0.0000000000000"+000  +0.0000000000000"+000
 
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 7
 
 
 
 SOURCE TEXT(S) :
 
 
 "CODE" 34501;
 "INTEGER""PROCEDURE" ZERPOL(N, A, EM, RE, IM, D);
   "VALUE" N; "INTEGER" N; "ARRAY" A, EM, RE, IM, D;
 "BEGIN" "INTEGER" I, TOTIT, IT, FAIL, START, UP, MAX, GIEX, ITMAX;
    "REAL" X, Y, NEWF, OLDF, MAXRAD, AE, TOL, H1, H2, LN2;
    "ARRAY" F[0 : 5], TRIES[1 : 10];
 
   "REAL""PROCEDURE" DWARF; "CODE" 30003;
   "REAL""PROCEDURE" GIANT; "CODE" 30004;
   "REAL""PROCEDURE" COMABS(XR, XI); "CODE" 34340;
   "PROCEDURE" COMSQRT(AR, AI, PR, PI); "CODE" 34343;
 
   "BOOLEAN""PROCEDURE" FUNCTION;
   "BEGIN" "INTEGER" K, M1, M2;
     "REAL" P, Q, QSQRT, F01, F02, F03, F11, F12, F13,
            F21, F22, F23, STOP;
     IT:= IT + 1;
     P:= 2 * X; Q:= -(X * X + Y * Y); QSQRT:= SQRT(-Q);
     F01:= F11:= F21:= D[0]; F02:= F12:= F22:= 0;
     M1:= N - 4; M2:= N - 2;
     STOP:= ABS(F01) * 0.8;
     "FOR" K:= 1 "STEP" 1 "UNTIL" M1 "DO"
     "BEGIN" F03:= F02; F02:= F01; F01:= D[K] + P * F02 + Q * F03;
       F13:= F12; F12:= F11; F11:= F01 + P * F12 + Q * F13;
       F23:= F22; F22:= F21; F21:= F11 + P * F22 + Q * F23;
       STOP:= QSQRT * STOP + ABS(F01)
     "END";
     "IF" M1 < 0 "THEN" M1:= 0;
     "FOR" K:= M1 + 1 "STEP" 1 "UNTIL" M2 "DO"
     "BEGIN" F03:= F02; F02:= F01; F01:= D[K] + P * F02 + Q * F03;
       F13:= F12; F12:= F11; F11:= F01 + P * F12 + Q * F13;
       STOP:= QSQRT * STOP + ABS(F01)
     "END";
     "IF" N = 3 "THEN" F21:= 0;
     F03:= F02; F02:= F01; F01:= D[N - 1] + P * F02 + Q * F03;
     F[0]:= D[N] + X * F01 + Q * F02;
     F[1]:= Y * F01;
     F[2]:= F01 - 2 * F12 * Y * Y;
     F[3]:= 2 * Y * (- X * F12 + F11);
     F[4]:= 2 * (- X * F12 + F11) - 8 * Y * Y * (- X * F22 + F21);
     F[5]:= Y * (6 * F12 - 8 * Y * Y * F22);
     STOP:= QSQRT * (QSQRT * STOP + ABS(F01)) + ABS(F[0]);
     NEWF:= F02:= COMABS(F[0], F[1]);
     FUNCTION:= F02 < (2 * ABS(X * F01) - 8 * (ABS(F[0]) + ABS(F01)
        * QSQRT) + 10 * STOP) * TOL * (1 + TOL) ** (4 * N + 3)
   "END" OF FUNCTION;                                          "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 8
                                                                  ;
 
 
   "BOOLEAN""PROCEDURE" CONTROL;
   "IF" IT > ITMAX "THEN"
   "BEGIN" TOTIT:= TOTIT + IT; FAIL:= 3; "GOTO" EXIT "END"
   "ELSE" "IF" IT = 0 "THEN"
   "BEGIN" "INTEGER" I, H; "REAL" H1, SIDE;
     MAXRAD:= 0; MAX:= (GIEX - LN(ABS(D[0])) / LN2) / N;
     "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
     "BEGIN" H1:= "IF" D[I] = 0 "THEN" 0
       "ELSE" EXP(LN(ABS(D[I] / D[0])) / I);
       "IF" H1 > MAXRAD "THEN" MAXRAD:= H1
     "END";
     "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
     "IF" D[I] ^= 0 "THEN"
     "BEGIN" H:= (GIEX - LN(ABS(D[I])) / LN2) / (N - I);
       "IF" H < MAX "THEN" MAX:= H
     "END";
     MAX:= MAX * LN2 / LN(N);
     SIDE:= - D[1] / D[0];
     SIDE:= "IF" ABS(SIDE) < TOL "THEN" 0 "ELSE" SIGN(SIDE);
     "IF" SIDE = 0 "THEN"
     "BEGIN" TRIES[7]:= TRIES[2]:= MAXRAD; TRIES[9]:= -MAXRAD;
       TRIES[6]:= TRIES[4]:= TRIES[3]:= MAXRAD / SQRT(2);
       TRIES[5]:= -TRIES[3]; TRIES[10]:= TRIES[8]:= TRIES[1]:= 0
     "END" "ELSE"
     "BEGIN" TRIES[8]:= TRIES[4]:= MAXRAD/ SQRT(2);
        TRIES[1]:= SIDE * MAXRAD; TRIES[3]:= TRIES[4] * SIDE;
        TRIES[6]:= MAXRAD; TRIES[7]:= -TRIES[3];
        TRIES[9]:= -TRIES[1]; TRIES[2]:= TRIES[5]:= TRIES[10]:= 0
     "END";
     "IF" COMABS(X, Y) > 2 * MAXRAD "THEN" X:= Y:= 0;
     CONTROL:= "FALSE"
   "END" "ELSE"
   "BEGIN" "IF" IT > 1 & NEWF >= OLDF "THEN"
     "BEGIN" UP:= UP+ 1;
       "IF" UP = 5 & START < 5 "THEN"
       "BEGIN" START:= START + 1; UP:= 0; X:= TRIES[2 * START - 1];
         Y:= TRIES[2 * START]; CONTROL:= "FALSE"
       "END" "ELSE" CONTROL:= "TRUE"
     "END" "ELSE" CONTROL:= "TRUE"
   "END" OF CONTROL;                                           "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 9
                                                                  ;
 
 
   "PROCEDURE" DEFLATION;
   "IF" X = 0 & Y = 0 "THEN" N:= N - 1 "ELSE"
   "BEGIN" "INTEGER" I, SPLIT; "REAL" H1, H2;
     "ARRAY" B[0 : N - 1];
     "IF" Y = 0 "THEN"
     "BEGIN" N:= N - 1; B[N]:= -D[N + 1] / X;
       "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
       B[N - I]:= (B[N - I + 1] - D[N - I + 1]) / X;
       "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
       D[I]:= D[I] + D[I - 1] * X
     "END" "ELSE"
     "BEGIN" H1:= - 2 * X; H2:= X * X + Y * Y;
       N:= N - 2;
       B[N]:= D[N + 2] / H2; B[N - 1]:= (D[N + 1] - H1 * B[N]) / H2;
       "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
       B[N - I]:=  (D[N - I + 2] - H1 * B[N - I + 1] - B[N - I + 2])/H2;
       D[1]:= D[1] - H1 * D[0];
       "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
       D[I]:= D[I] - H1 * D[I-1] - H2 * D[I-2]
     "END";
     SPLIT:= N;
     H2:= ABS(D[N] - B[N]) / (ABS(D[N]) + ABS(B[N]));
     "FOR" I:= N - 1 "STEP" -1 "UNTIL" 0 "DO"
     "BEGIN" H1:= ABS(D[I]) + ABS(B[I]);
       "IF" H1 > TOL "THEN"
       "BEGIN" H1:= ABS(D[I] - B[I]) / H1;
         "IF" H1 < H2 "THEN" "BEGIN" H2:= H1; SPLIT:= I "END"
       "END"
     "END";
     "FOR" I := SPLIT + 1 "STEP" 1 "UNTIL" N "DO" D[I]:= B[I];
     D[SPLIT]:= (D[SPLIT] + B[SPLIT]) / 2
   "END" OF DEFLATION;
 
   "PROCEDURE" LAGUERRE;
   "BEGIN" "INTEGER" M;
     "REAL" S1RE, S1IM, S2RE, S2IM, DX, DY, H1, H2, H3, H4, H5, H6;
     "IF" ABS(F[0]) > ABS(F[1]) "THEN"
     "BEGIN" H1:= F[0]; H6:= F[1] / H1; H2:= F[2] + H6 * F[3];
       H3:= F[3] - H6 * F[2]; H4:= F[4] + H6 * F[5];
       H5:= F[5] - H6 * F[4]; H6:= H6 * F[1] + H1
     "END" "ELSE"
     "BEGIN" H1:= F[1]; H6:= F[0] / H1; H2:= H6 * F[2] + F[3];
       H3:= H6 * F[3] - F[2]; H4:= H6 * F[4] + F[5];
       H5:= H6 * F[5] - F[4]; H6:= H6 * F[0] + F[1]
     "END";                                                    "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 10
                                                                  ;
 
 
     S1RE:= H2 / H6; S1IM:= H3 / H6;
     H2:= S1RE * S1RE - S1IM * S1IM; H3:= 2 * S1RE * S1IM;
     S2RE:= H2 - H4 / H6; S2IM:= H3 - H5 / H6;
     H1:= S2RE * S2RE + S2IM * S2IM;
     H1:= "IF" H1 ^= 0 "THEN" (S2RE * H2 + S2IM * H3) / H1 "ELSE" 1;
     M:= "IF" H1 >= N - 1 "THEN" ("IF" N > 1 "THEN" N - 1 "ELSE" 1)
         "ELSE" "IF" H1 > 1 "THEN" H1 "ELSE" 1;
     H1:= (N - M) / M;
     COMSQRT(H1 * (N * S2RE - H2), H1 * (N * S2IM - H3), H2, H3);
     "IF" S1RE * H2 + S1IM * H3 < 0 "THEN"
     "BEGIN" H2:= - H2; H3:= - H3 "END";
     H2:= S1RE + H2; H3:= S1IM + H3;
     H1:= H2 * H2 + H3 * H3;
     "IF" H1 = 0 "THEN" "BEGIN" DX:= -N; DY:= N "END" "ELSE"
     "BEGIN" DX:= - N * H2 / H1; DY:= N * H3 / H1 "END";
     H1:= ABS(X) * TOL + AE; H2:= ABS(Y) * TOL+ AE;
     "IF" ABS(DX) < H1 & ABS(DY) < H2 "THEN"
     "BEGIN" DX:= "IF" DX = 0 "THEN" H1 "ELSE" SIGN(DX) * H1;
       DY:= "IF" DY = 0 "THEN" H2 "ELSE" SIGN(DY) * H2
     "END";
     X:= X + DX; Y:= Y + DY;
     "IF" COMABS(X, Y) > 2 * MAXRAD "THEN"
     "BEGIN" H1:= "IF" ABS(X) > ABS(Y) "THEN" ABS(X) "ELSE" ABS(Y);
       H2:= LN(H1) / LN2 + 1 - MAX;
       "IF" H2 > 0 "THEN"
       "BEGIN" H2:= 2 ** H2; X:= X / H2; Y:= Y / H2 "END"
     "END"
   "END" OF LAGUERRE;
 
   TOTIT:= IT:= FAIL:= UP:= START:= 0; LN2:= LN(2);
   NEWF:= GIANT; AE:= DWARF; GIEX:= LN(NEWF) / LN2 - 40;
   TOL:= EM[0]; ITMAX:= EM[1];
   "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" D[I]:= A[N-I];
   "IF" N <= 0 "THEN"
   "BEGIN" FAIL:= 1; "GOTO" EXIT "END"
   "ELSE" "IF" D[0] = 0 "THEN"
   "BEGIN" FAIL:= 2; "GOTO" EXIT "END";
   "FOR" I:= 1 "WHILE" D[N] = 0 & N > 0 "DO"
   "BEGIN" RE[N]:= IM[N]:= 0; N:= N - 1 "END";
   X:= Y:= 0;                                                  "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 11
                                                                  ;
 
 
   "FOR" I:= 1 "WHILE" N > 2 "DO"
   "BEGIN" "IF" CONTROL "THEN" LAGUERRE;
     OLDF:= NEWF;
     "IF" FUNCTION "THEN"
     "BEGIN" "IF" Y ^= 0 & ABS(Y) < .1 "THEN"
       "BEGIN" "REAL" H; H:= Y; Y:= 0;
          "IF" ^ FUNCTION "THEN" Y:= H
       "END";
       RE[N]:= X; IM[N]:= Y;
       "IF" Y ^= 0 "THEN" "BEGIN" RE[N - 1]:= X; IM[N - 1]:= -Y "END";
       DEFLATION; TOTIT:= TOTIT + IT; UP:= START:= IT:= 0
     "END"
   "END";
   "IF" N = 1 "THEN" "BEGIN" RE[1]:= - D[1] / D[0]; IM[1]:= 0 "END"
   "ELSE"
   "BEGIN" "REAL" H1, H2;
     H1:= - 0.5 * D[1] / D[0]; H2:= H1 * H1 - D[2] / D[0];
     "IF" H2 >= 0 "THEN"
     "BEGIN" RE[2]:= "IF" H1 < 0 "THEN" H1 - SQRT(H2)
             "ELSE" H1 + SQRT(H2);
       RE[1]:= D[2] / (D[0] * RE[2]);
       IM[2]:= IM[1]:= 0
     "END" "ELSE"
     "BEGIN" RE[2]:= RE[1]:= H1;
       IM[2]:= SQRT(-H2); IM[1]:= -IM[2]
     "END"
   "END"; N:= 0;
 EXIT:  EM[2]:= FAIL; EM[3]:= START; EM[4]:= TOTIT;
   "FOR" I:= (N-1) "DIV" 2 "STEP" -1 "UNTIL" 0 "DO"
   "BEGIN" TOL := D[I]; D[I]:= D[N-I]; D[N-I]:= TOL
   "END";
   ZERPOL:= N
 "END" OF ZERPOL;
         "EOP"
 
 
 "CODE" 34502;
 "PROCEDURE" BOUNDS(N,A,RE,IM,RELE,ABSE,RECENTRE,IMCENTRE,BOUND);
   "VALUE" N, RELE, ABSE; "INTEGER" N; "REAL" RELE, ABSE;
   "ARRAY" RE, IM, A, RECENTRE, IMCENTRE, BOUND;
 "BEGIN" "INTEGER" I, J, K, L, INDEX1, INDEX2; "BOOLEAN" GOON;
   "REAL" H, MIN, RECENT, IMCENT, GIA, XK, YK, ZK, CORR;
   "ARRAY" RC, C, RCE[0:N], CLUST[1:N];
 
   "REAL""PROCEDURE" GIANT; "CODE" 30004;
   "REAL""PROCEDURE" ARREB; "CODE" 30002;                      "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 12
                                                                  ;
 
 
 
   "REAL""PROCEDURE" G(RAD, RECENT, IMCENT, K, M);
     "VALUE" RAD, RECENT, IMCENT, K, M; "REAL" RAD, RECENT, IMCENT;
     "INTEGER" K, M;
   "BEGIN" "REAL" S, H1, H2; "INTEGER" I;
     S:= SQRT(RECENT * RECENT + IMCENT * IMCENT) + RAD;
     H1:= RC[1]; H2:= RC[0];
     "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" H1:= H1*S + RC[I];
     "FOR" I:= 1 "STEP" 1 "UNTIL" M-1, M+K "STEP" 1 "UNTIL" N "DO"
     H2:= H2 * ABS(SQRT((RE[I]-RECENT)**2 + (IM[I]-IMCENT)**2) - RAD);
     G:= "IF" H1=0 "THEN" 0 "ELSE" "IF" H2=0 "THEN" -10 "ELSE" H1 / H2
   "END";
 
   "PROCEDURE" KCLUSTER(K, M);
     "VALUE" K, M; "INTEGER" K, M;
   "BEGIN" "INTEGER" I, J, STOP, L; "BOOLEAN" NONZERO;
     "REAL" RECENT, IMCENT, D, PROD, RAD, GR, R;
     "ARRAY" DIST[M: M+K-1];
     RECENT:= RE[M]; IMCENT:= IM[M]; STOP:= M+K-1;
     L:= SIGN(IMCENT); NONZERO:= L ^= 0;
     "FOR" I:= M+1 "STEP" 1 "UNTIL" STOP "DO"
     "BEGIN" RECENT:= RECENT+RE[I];
       "IF" NONZERO "THEN"
       "BEGIN" NONZERO:= L = SIGN(IM[I]); IMCENT:= IMCENT+IM[I] "END"
     "END";
     RECENT:= RECENT/K; IMCENT:= "IF" NONZERO "THEN" IMCENT/K "ELSE" 0;
     D:= 0; RAD:= 0;
     "FOR" I:= M "STEP" 1 "UNTIL" STOP "DO"
     "BEGIN" RECENTRE[I]:= RECENT; IMCENTRE[I]:= IMCENT;
       DIST[I]:= SQRT((RE[I] -RECENT)**2 + (IM[I]-IMCENT)**2);
       "IF" D < DIST[I] "THEN" D:= DIST[I]
     "END";
     GR:= ABS(G(0, RECENT, IMCENT, K, M));
     "IF" GR > 0 "THEN"
     "BEGIN" "FOR" J:= 1, 1 "WHILE" PROD <= GR "DO"
       "BEGIN" R:= RAD; RAD:= D + EXP(LN(1.1*GR)/K);
         "IF" RAD = R "THEN" RAD:= EXP(LN(1.1)/K) * RAD;
         GR:= G(RAD, RECENT, IMCENT, K, M);
         PROD:= 1;
         "FOR" I:= M "STEP" 1 "UNTIL" STOP "DO"
         PROD:= PROD*(RAD-DIST[I])
       "END"
     "END";
     "FOR" I:= M "STEP" 1 "UNTIL" STOP "DO"
     "BEGIN" BOUND[I]:= RAD; CLUST[I]:= K "END";
   "END";                                                      "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1979)                      PAGE 13
                                                                  ;
 
 
   "PROCEDURE" SHIFT(INDEX, NEW);
     "VALUE" INDEX, NEW; "INTEGER" INDEX, NEW;
   "BEGIN" "INTEGER" J, PLACE, CLUSTIN;
     "REAL" BOUNDIN, IMCENT, RECENT;
     "REAL""ARRAY" WA1, WA2[1:CLUST[INDEX]];
     CLUSTIN:= CLUST[INDEX]; BOUNDIN:= BOUND[INDEX];
     IMCENT:= IMCENTRE[INDEX]; RECENT:= RECENTRE[INDEX];
     "FOR" J:= 1 "STEP" 1 "UNTIL" CLUSTIN "DO"
     "BEGIN" PLACE:=INDEX+J-1; WA1[J]:= RE[PLACE]; WA2[J]:= IM[PLACE];
     "END";
     "FOR" J:= INDEX-1 "STEP" -1 "UNTIL" NEW "DO"
     "BEGIN" PLACE:= J+CLUSTIN;
       RE[PLACE]:= RE[J]; IM[PLACE]:= IM[J]; CLUST[PLACE]:= CLUST[J];
       BOUND[PLACE]:= BOUND[J]; RECENTRE[PLACE]:= RECENTRE[J];
       IMCENTRE[PLACE]:= IMCENTRE[J]
     "END";
     "FOR" J:= NEW+CLUSTIN-1 "STEP" -1 "UNTIL" NEW "DO"
     "BEGIN" PLACE:= J+1-NEW;
       RE[J]:= WA1[PLACE]; IM[J]:= WA2[PLACE];
       BOUND[J]:= BOUNDIN; CLUST[J]:= CLUSTIN;
       RECENTRE[J]:= RECENT; IMCENTRE[J]:= IMCENT
     "END"
   "END";
 
   GIA:= GIANT;
   RC[0]:= C[0]:= A[N]; RCE[0]:= ABS(C[0]); K:= 0;
   "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
   "BEGIN" RC[I]:= RCE[I]:= 0 ; C[I]:= A[N-I] "END";
   "FOR" I:= 0 "WHILE" K < N "DO"
   "BEGIN" K:= K+1; XK:= RE[K]; YK:= IM[K]; ZK:= XK*XK+YK*YK;
     "FOR" J:= K "STEP" -1 "UNTIL" 1 "DO"
     RCE[J]:= RCE[J]+RCE[J-1]*SQRT(ZK);
     "IF" YK = 0 "THEN"
     "BEGIN" "FOR" J:= K "STEP" -1 "UNTIL" 1 "DO"
        RC[J]:= RC[J]-XK*RC[J-1]
     "END" "ELSE"
     "BEGIN" K:= K+1;
       "IF" K <= N & XK = RE[K] & YK = -IM[K] "THEN"
       "BEGIN" XK:= -2*XK;
         "FOR" J:= K "STEP" -1 "UNTIL" 1 "DO"
         RCE[J]:= RCE[J]+RCE[J-1]*SQRT(ZK);
         "FOR" J:= K "STEP" -1 "UNTIL" 2 "DO"
         RC[J]:= RC[J]+XK*RC[J-1]+ZK*RC[J-2];
         RC[1]:= RC[1]+XK*RC[0]
       "END"
     "END"
   "END";
   RC[0]:= RCE[0]; CORR:= 1.06*ARREB;
   "FOR" I:= 1 "STEP" 1 "UNTIL" N-1 "DO"
   RC[I]:= ABS(RC[I]-C[I])+RCE[I]*CORR*(N+I-2)+RELE*ABS(C[I])+ABSE;
   RC[N]:= ABS(RC[N]-C[N])+RCE[N]*CORR*(N-1)+RELE*ABS(C[N])+ABSE;
   "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" KCLUSTER(1, I);         "COMMENT"
1SECTION : 3.6.1             (DECEMBER 1978)                      PAGE 14
                                                                  ;
 
 
   GOON:= "TRUE";
   "FOR" L:= 1 "WHILE" GOON "DO"
   "BEGIN" INDEX1:= INDEX2:= 0; MIN:= GIANT; I:= N-CLUST[N]+1;
     "FOR" I:= I "WHILE" I >= 2 "DO"
     "BEGIN" J:= I; RECENT:= RECENTRE[I]; IMCENT:= IMCENTRE[I];
       "FOR" J:= J "WHILE" J >= 2 "DO"
       "BEGIN" J:= J-CLUST[J-1];
         H:= SQRT((RECENT-RECENTRE[J])**2 + (IMCENT-IMCENTRE[J])**2);
         "IF" H < BOUND[I] + BOUND[J] & H <= MIN "THEN"
         "BEGIN" INDEX1:= J; INDEX2:= I; MIN:= H "END"
       "END"; I:= I-CLUST[I-1]
     "END";
     "IF" INDEX1 = 0 "THEN" GOON:= "FALSE" "ELSE"
     "BEGIN" "IF" IMCENTRE[INDEX1] = 0 "THEN"
       "BEGIN" "IF" IMCENTRE[INDEX2] ^= 0 "THEN"
         CLUST[INDEX2]:= 2*CLUST[INDEX2]
       "END" "ELSE" "IF" IMCENTRE[INDEX2] = 0 "THEN"
       CLUST[INDEX1]:= 2*CLUST[INDEX1];
       K:= INDEX1+CLUST[INDEX1];
       "IF" K ^= INDEX2 "THEN" SHIFT(INDEX2, K);
       K:= CLUST[INDEX1]+CLUST[K];
       KCLUSTER(K, INDEX1)
     "END"
   "END"
 "END";
         "EOP"
1SECTION: 4.2.3.1            (NOVEMBER 1978)                      PAGE 1      
        
        
        
 AUTHOR:         C.G. VAN DER LAAN.   
        
        
 CONTRIBUTORS:   C.G. VAN DER LAAN AND M. VOORINTHOLT.    
        
        
 INSTITUTE:      REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN.
        
        
 RECEIVED:       780601.    
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS THE PROCEDURES GSSWTS, GSSWTSSYM AND RECCOF.       
     RECCOF CALCULATES FROM A GIVEN WEIGHT FUNCTION ON [-1,1] THE   
     RECURRENCE COEFFICIENTS OF THE CORRESPONDING ORTHOGONAL
     POLYNOMIALS; GSSWTS AND GSSWTSSYM CALCULATE FROM THE RECURRENCE
     COEFFICIENTS THE GAUSSIAN WEIGHTS OF THE CORRESPONDING WEIGHT  
     FUNCTION.    
        
        
 KEYWORDS:
        
     RECURRENCE COEFFICIENTS ORTHOGONAL POLYNOMIALS,      
     GAUSSIAN WEIGHTS,      
     GAUSSIAN QUADRATURE.   
        
        
 SUBSECTION: RECCOF.
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" RECCOF(N,M,X,WX,B,C,L,SYM);    
     "VALUE" N,M,SYM; "INTEGER" N,M; "BOOLEAN" SYM;       
     "REAL" X,WX; "ARRAY" B,C,L;      
     "CODE" 31254;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:       <ARITHMETIC EXPRESSION>;
              ENTRY: UPPER BOUND FOR THE INDICES OF THE ARRAYS B, C, L
                     (N>=0);
     M:       <ARITHMETIC EXPRESSION>;
              ENTRY: THE NUMBER OF POINTS USED IN THE GAUSS-CHEBYSHEV 
                     QUADRATURE RULE FOR CALCULATING THE APPROXIMATION
                     OF THE INTEGRAL REPRESENTATIONS OF B[K],C[K]   
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 2      
        
        
        
     SYM:     <BOOLEAN EXPRESSION>;   
              ENTRY: "IF" SYM 
                         "THEN" WEIGHT FUNCTION ON [-1,1] IS EVEN   
                         "ELSE" WEIGHT FUNCTION ON [-1,1] IS NOT EVEN;
     X,WX:    <ARITHMETIC EXPRESSION>;
              ENTRY: JENSEN VARIABLES WITH WX AN EXPRESSION OF X    
                     DENOTING THE WEIGHT FUNCTION ON [-1,1];
     B,C,L:   <ARRAY IDENTIFIER>;     
              "ARRAY" B,C,L[0:N];     
              EXIT:  THE APPROXIMATE RECURRENCE COEFFICIENTS FOR    
                     P[K+1](X) = (X-B[K])*P[K](X) - C[K]*P[K-1](X), 
                                                      K=0,1,2,...N, 
                     AND THE APPROXIMATE SQUARE LENGTHS   
                             X = +1   
                     L[K] = INTEGRAL ( W(X) * P[K](X) **  2 ) DX    
                             X = -1   
        
        
 PROCEDURES USED:    ORTPOL = CP31044.
        
        
 RUNNING TIME:       PROPORTIONAL TO M*N**2.    
        
        
 METHOD AND PERFORMANCE:    
        
     THE RECURRENCE COEFFICIENTS ARE REPRESENTED BY       
                        X = +1
              B[K] = ( INTEGRAL ( W(X) * X * P[K](X) ** 2 ) DX ) / L[K],      
                        X = -1
        
              C[K] = L[K] / L[K-1]],  
     WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL.     
     THE INTEGRALS ARE APPROXIMATED BY THE M-POINTS GAUSS-CHEBYSHEV 
     RULE AS      
           X = +1                       M       
          INTEGRAL(F(X))DX := PI / M * SUM SIN(THETA[J])*F(COS(THETA[J]))     
           X = -1                      J=1      
     WITH THETA[J] = (2*J-1) * PI / (2*M)      (SEE GAUTSCHI, 1968A). 
     THE VALUE OF M IS TO BE SUPPLIED BY THE USER.
        
        
 REFERENCES:      
        
     GAUTSCHI, W. (1968A):  
     CONTRUCTION OF GAUSS-CHRISTOFFEL FORMULAS. 
     MATH. COMP., 22,P.251-270.       
        
     GAUTSCHI, W. (1968B):  
     GAUSSIAN QUADRATURE FORMULAS.    
     COMM. ACM. CALGO 331.  
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 3      
        
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING PROGRAM DELIVERS AN APPROXIMATION      
     FOR THE RECURSION COEFFICIENTS C[1] AND C[2], OF THE CHEBYSHEV 
     POLYNOMIALS OF THE SECOND KIND;  
        
     "BEGIN"      
         "PROCEDURE" RECCOF(N,M,X,WX,B,C,L,SYM);
         "VALUE" N,M,SYM; "INTEGER" N,M; "BOOLEAN" SYM;   
         "REAL" X,WX; "ARRAY" B,C,L;  
         "CODE" 31254;      
         "REAL" X; "ARRAY" B,C,L[0:2];
         RECCOF(2,200,X,SQRT(1 - X**2),B,C,L,"TRUE");     
         OUTPUT(61,"("2/,2(3B,-ZD.3D)")",C[1],C[2]);      
     "END";       
        
     RESULTS:     
        
          0.250     0.250   
        
        
 SUBSECTION: GSSWTS.
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" GSSWTS(N,ZER,B,C,W); 
     "VALUE" N; "INTEGER" N; "ARRAY" ZER,B,C,W; 
     "CODE" 31253;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:       <ARITHMETIC EXPRESSION>;
              ENTRY: THE NUMBER OF WEIGHTS TO BE COMPUTED; UPPER    
                     BOUND FOR THE ARRAYS Z AND W (N>=1); 
     ZER:     <ARRAY IDENTIFIER>;     
              "ARRAY" ZER[1:N];       
              ENTRY: THE ZEROS OF THE N-TH DEGREE ORTHOGONAL POLYNOMIAL;      
     B,C:     <ARRAY IDENTIFIER>;     
              "ARRAY" B[0:N-1], C[1:N-1];       
              ENTRY: THE RECURRENCE COEFFICIENTS; 
     W:       <ARRAY IDENTIFIER>;     
              "ARRAY" W[1:N]; 
              EXIT : THE GAUSSIAN WEIGHTS DIVIDED BY THE  
                     INTEGRAL OVER THE WEIGHT FUNCTION.   
        
        
 PROCEDURES USED:    ALLORTPOL = CP 31045.      
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 4      
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE GAUSSIAN WEIGHTS OF AN N-POINTS RULE DIVIDED BY THE INTEGRAL 
     OF THE WEIGHT FUNCTION MAY BE REPRESENTED AS 
              W[K] = 1/(...((P[N-1](Z)**2/C[N-1]+P[N-2](Z)**2)/C[N-2]+...     
                                  ...+P[1](Z)**2)/C[1]+1)  , K=1,2,...N       
     WITH     Z= K-TH ZERO OF P[N](X).  (GAUTSCHI, 1970). 
        
     ALLZERORTPOL AND GSSWTS MAY BE USED TO GENERATE GAUSSIAN       
     QUADRATURE RULES PROVIDED THE RECURRENCE COEFFICIENTS AND THE  
     INTEGRAL OF THE WEIGHT FUNCTION ARE KNOWN. 
     FOR EXAMPLE THE GAUSS-LAGUERRE QUADRATURE RULE APPLIED TO F    
     MAY BE OBTAINED BY THE CALLS     
        
              "FOR" K:=1 "STEP" 1 "UNTIL" N-1 "DO"
              "BEGIN"       
                 B[K]:=2*K+ALPHA+1;   
                 C[K]:=K*(K+ALPHA)    
              "END";
              B[0]:=ALPHA+1;
              ALLZERORTPOL(N,ZER,B,C);
              GSSWTS(N,ZER,B,C,W);    
              GAUSSRULE:=0; 
              "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"
              GAUSSRULE:=GAUSSRULE+W[K]*F(ZER[K]);
              GAUSSRULE:=GAUSSRULE*GAMMA(ALPHA+1) 
        
     GAUSSRULE CONTAINS THE VALUE OF THE APPOXIMATING GAUSS 
     QUADRATURE RULE AND ZER[1:N],W[1:N] CONTAIN THE GAUSSIAN       
     ABSCISSA AND WEIGHTS.  
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 5      
        
        
        
     IN THE FOLLOWING TABEL WE SUMMARIZE CLASSICAL QUADRATURE RULES 
        
            :  WEIGHT   :    RECURRENCE COEFFICIENTS    :   INTEGRAL
  GAUSSIAN  : FUNCTION  :---------------:---------------:      OF   
 QUADRATURE :   W(X)    :      B[K]     :     C[K]      :WEIGHT FUNCTION      
 -----------:-----------:---------------:---------------:---------------      
            :           :               :               : 
 LEGENDRE   :     1     :       0       :K**2*(4*K**2-1):      2    
            :           :               :               : 
  CHEBYSHEV : 1/SQRT(1- :       0       :  1/2  ,  K=1  :      PI   
 (1-ST KIND): X**2)     :               :  1/4  ,  K>1  : 
            :           :               :               : 
  CHEBYSHEV : SQRT(1-   :       0       :      1/4      :    PI/2   
 (2-ND KIND): X**2)     :               :               : 
            :           :               :               : 
 JACOBI     : (1-X)**   : -(ALPHA**2-   : 4*(1+ALPHA)*  : 2**(ALPHA+
            : ALPHA*(1+ : BETA**2)/((2* : (1+BETA)/     : BETA+1)*  
            : X)**BETA  : K+ALPHA+BETA)*: ((ALPHA+BETA+ : GAMMA(ALPHA+
            :           : (2*K+ALPHA+   : 2)**2*(ALPHA+ : 1)*GAMMA(BETA+      
            :           : BETA+2))      : BETA+3)) ,K=1 : 1)/GAMMA(ALPHA+     
            :           :               :               : BETA+2)   
            :           :               : 4*K*(K+ALPHA)*: 
            :           :               : (K+BETA)*(K+  : 
            :           :               : ALPHA+BETA)/  : 
            :           :               : ((2*K+ALPHA+  : 
            :           :               : BETA)**2*     : 
            :           :               : ((2*K+ALPHA+  : 
            :           :               : BETA)**2-1))  : 
            :           :               :          ,K>1 : 
            :           :               :               : 
            :           :               :(ALPHA,BETA>-1): 
            :           :               :               : 
 LAGUERRE   : EXP(-X)*  : 2*K+ALPHA+1   : K*(K+ALPHA)   : GAMMA(ALPHA+1)      
            : X**ALPHA  :               :               : 
            :           :               :               : 
 HERMITE    : EXP(-X**2):       0       :      K/2      :   SQRT(PI)
        
 (THE INTEGRATION INTERVALS ARE:  [-INFINITY,INFINITY]  FOR HERMITE;
                                          [0,INFINITY]  FOR LAGUERRE; 
                                         [-1,1]         FOR THE OTHERS.)      
 FOR NON-CLASSICAL WEIGHT FUNCTIONS ON A FINITE INTERVAL THE RECURSION
 COEFFICIENTS (AND THE SQUARE LENGTHS OF THE CORRESPONDING ORTHOGONAL-
 POLYNOMIALS) CAN BE OBTAINED BY THE PROCEDURE RECCOF (THIS SECTION). 
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 6      
        
        
        
 REFERENCES:      
        
     GAUTSCHI, W. (1970):   
     GENERATION OF GAUSSIAN QUADRATURE RULES AND
     ORTHOGONAL POLYNOMIALS.
     IN:  COLLOQUIUM APPROXIMATIETHEORIE,       
          MC SYLLABUS 14.   
        
        
 EXAMPLE OF USE: SEE SUBSECTION GSSWTSSYM.      
        
        
 SUBSECTION: GSSWTSSYM.     
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" GSSWTSSYM(N,ZER,C,W);
     "VALUE" N; "INTEGER" N; "ARRAY" ZER,C,W;   
     "CODE" 31252;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:       <ARITHMETIC EXPRESSION>;
              ENTRY: THE WEIGHTS OF AN N-POINTS GAUSS RULE ARE      
                     TO BE CALCULATED (BECAUSE OF SYMMETRY ONLY     
                     (N+1)//2 OF THE VALUES ARE DELIVERED); 
     ZER:     <ARRAY IDENTIFIER>;     
              "ARRAY" ZER[1:N//2]     
              ENTRY: THE NEGATIVE ZEROS OF THE N-TH DEGREE ORTHOGONAL 
                     POLYNOMIAL (ZER[I] < ZER[I+1], I=1,2,...,N//2-1);
              (IF N IS ODD THEN 0 IS ALSO A ZERO.)
     C:       <ARRAY IDENTIFIER>;     
              "ARRAY" C[1:N-1];       
              ENTRY: THE RECURRENCE COEFFICIENTS; 
     W:       <ARRAY IDENTIFIER>;     
              "ARRAY" W[1:(N+1)//2];  
              EXIT : PART OF THE GAUSSIAN WEIGHTS DIVIDED BY THE    
                     INTEGRAL OF THE WEIGHT FUNCTION.     
              (NOTE THAT W[N+1-K]=W[K]  , K=1,2,...,(N+1)//2.)      
        
        
 PROCEDURES USED:    ALLORTPOLSYM = CP 31049.   
        
        
 METHOD AND PERFORMANCE: SEE SUBSECTION GSSWTS; THIS PROCEDURE IS   
                         SUPPLIED FOR STORAGE ECONOMICAL REASONS.   
        
        
 REFERENCES:  SEE SUBSECTION GSSWTS.  
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 7      
        
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING PROGRAM DELIVERS THE GAUSSIAN WEIGHTS  
     FOR THE 5-POINTS GAUSS-CHEBYSHEV QUADRATURE RULE BY MEANS OF   
     THE PROCEDURE GSSWTSSYM (C[1]=0.5; C[K]=0.25, K=2,3,...;       
     ZER[I] = COS((2*(N-I) - 1) / (2*N) * PI), I=1,2,...,N//2.      
        
     "BEGIN"      
         "PROCEDURE" GSSWTSSYM(N,ZER,C,W);      
         "VALUE" N; "INTEGER" N; "ARRAY" ZER,C,W; 
         "CODE" 31252;      
         "REAL" PI; "INTEGER" I;      
         "ARRAY" ZER[1:2], W[1:3], C[1:4];      
         PI:=4*ARCTAN(1);   
         C[1]:=.5;
         "FOR" I:=2 "STEP" 1 "UNTIL" 4 "DO"     
         C[I]:=.25; 
         ZER[1]:=COS(.9*PI);
         ZER[2]:=COS(.7*PI);
         GSSWTSSYM(5,ZER,C,W);
         OUTPUT(61,"("2/,5(3B,-ZD.3D)")",W[1]*PI,W[2]*PI,W[3]*PI,   
                                                 W[2]*PI,W[1]*PI);  
     "END";       
        
     RESULTS:     
          0.628     0.628     0.628     0.628     0.628   
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 8      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 31254;    
 "PROCEDURE" RECCOF(N,M,X,WX,B,C,L,SYM);
 "VALUE" N,M,SYM;"INTEGER" N,M; "BOOLEAN" SYM;  
 "REAL" X,WX;"ARRAY" B,C,L; 
 "BEGIN" "INTEGER" I,J,UP;"REAL" R,S,PIM,H,HH,ARG,SA;     
 "REAL""PROCEDURE" ORTPOL(N,X,B,C);   
 "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C; 
 "CODE" 31044;    
 PIM:=4*ARCTAN(1)/M;
 "IF" SYM "THEN" "BEGIN"    
 "FOR" J:=0 "STEP" 1 "UNTIL" N"DO"    
 "BEGIN" R:=B[J]:=0;UP:=M "DIV" 2;    
    "FOR" I:=1 "STEP" 1 "UNTIL" UP"DO"
    "BEGIN" ARG:=(I-.5)*PIM;X:=COS(ARG);
       R:=R+SIN(ARG)*WX*ORTPOL(J,X,B,C)**2;     
    "END";"IF" UP*2=M "THEN" L[J]:=2*R*PIM "ELSE" 
    "BEGIN" X:=0;L[J]:=(2*R+WX*ORTPOL(J,0,B,C)**2)*PIM "END";       
    C[J]:="IF" J=0 "THEN" 0 "ELSE" L[J]/L[J-1]; 
 "END" "END" "ELSE" 
 "FOR" J:=0 "STEP" 1 "UNTIL" N "DO"   
 "BEGIN" R:=S:=0;UP:=M"DIV" 2;
    "FOR" I:=1 "STEP" 1 "UNTIL" UP "DO" 
    "BEGIN" ARG:=(I-.5)*PIM;SA:=SIN(ARG);X:=COS(ARG);     
       H:=WX*ORTPOL(J,X,B,C)**2;X:=-X;HH:=WX*ORTPOL(J,X,B,C)**2;    
       R:=R+(H+HH)*SA;S:=S+(HH-H)*X*SA; 
    "END";B[J]:=S*PIM;      
    "IF" UP*2=M "THEN" L[J]:=R*PIM"ELSE"
    "BEGIN" X:=0;L[J]:=(R+WX*ORTPOL(J,0,B,C)**2)*PIM "END"; 
    C[J]:="IF" J=0 "THEN" 0 "ELSE" L[J]/L[J-1]; 
 "END"; 
 "END" RECCOF;    
         "EOP"    
1SECTION 4.2.3.1             (NOVEMBER 1978)                      PAGE 9      
        
        
        
 "CODE" 31253;    
 "PROCEDURE" GSSWTS(N,ZER,B,C)RESULTS:(W);      
 "VALUE" N; "INTEGER" N;    
 "ARRAY" ZER,B,C,W; 
 "BEGIN"
     "INTEGER" J,K; "REAL" S; "ARRAY" P[0:N-1]; 
     "PROCEDURE" ALLORTPOL(N,X,B,C,P);
     "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C,P;   
     "CODE" 31045;
     "FOR" J:=1 "STEP" 1 "UNTIL" N "DO" 
     "BEGIN"      
         ALLORTPOL(N-1,ZER[J],B,C,P); 
         S:=0.0;  
         "FOR" K:=N-1 "STEP" -1 "UNTIL" 1 "DO"  
         S:=(S+P[K]**2)/C[K]; 
         W[J]:=1/(1+S);     
     "END"
 "END" GSSWTS;    
         "EOP"    
        
        
 "CODE" 31252;    
 "PROCEDURE" GSSWTSSYM(N,ZER,C)RESULTS:(W);     
 "VALUE" N; "INTEGER" N;    
 "ARRAY" ZER,C,W; 
 "BEGIN"
    "PROCEDURE" ALLORTPOLSYM(N,X,C,P);
    "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C,P;      
    "CODE" 31049; 
    "INTEGER" LOW,UP,DUMMY; 
    "ARRAY" P[0:N-1];       
    LOW:=1; UP:=N;
       "FOR" DUMMY:=1       
    "WHILE" LOW < UP "DO"   
       "BEGIN" "INTEGER" I; "REAL" S; 
       ALLORTPOLSYM( N-1,ZER[LOW],C,P );
       S:=P[N-1]**2;
       "FOR" I:=N-1 "STEP" -1 "UNTIL" 1 "DO"    
          S:=S/C[I] + (P[I-1])**2;    
       W[LOW]:=1/S; LOW:=LOW+1; UP:=UP-1;       
       "END";     
    "IF" LOW = UP 
    "THEN" "BEGIN" "INTEGER" TWOI; "REAL" S; S:=1.0;      
       "FOR" TWOI:=N-1 "STEP" -2 "UNTIL" 2 "DO" 
          S:=S*C[TWOI-1]/C[TWOI]+1;   
       W[LOW]:=1/S; 
       "END";     
 "END" GSSWTSSYM; 
         "EOP"    
1SECTION : 6.4.3             (DECEMBER 1978)                     PAGE 1       
        
        
        
 AUTHOR : N.M. TEMME.       
        
        
 CONTRIBUTOR : R. MONTIJN.  
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED : 780801. 
        
        
 BRIEF DESCRIPTION :
        
     THIS SECTION CONTAINS THE PROCEDURE  LOG ONE PLUS X  FOR       
     THE COMPUTATION OF  LN(1+X) FOR X > -1.    
        
        
 KEYWORDS : LOGARITHMIC FUNCTION.     
        
        
 CALLING SEQUENCE : 
        
     THE HEADING OF THE PROCEDURE READS :       
     "REAL" "PROCEDURE" LOG ONE PLUS X(X); "VALUE" X; "REAL" X;     
     "CODE" 35130;
     LOG ONE PLUS X : DELIVERS THE VALUE OF LN(1+X);      
        
     THE MEANING OF THE FORMAL PARAMETER IS :   
     X:      <ARITHMETIC  EXPRESSION>;
             ENTRY :  THE ARGUMENT OF LN(1+X), X > -1.    
        
        
 PROCEDURES USED : NONE.    
        
        
 RUNNING TIME : THE ALGORITHM NEEDS 9 MULTIPLICATIONS.    
        
        
 METHOD AND PERFORMANCE :   
        
     FOR  X < -0.2928 OR X > 0.4142  THE PROCEDURE USES   
     THE STANDARD FUNCTION LN, FOR -0.2928 <= X <= 0.4142 A 
     POLYNOMIAL APPROXIMATION IS USED.
     WE USE AN APPROXIMATION BASED ON THE BEST APPROXIMATON FOR     
     THE INTERVAL 1/SQRT(2)-1 <= X <= SQRT(2)-1, OF WHICH THE       
     COEFFICIENTS ARE GIVEN IN HART (1968); CF. P. 111, INDEX 2665. 
     THE PROCEDURE LOG ONE PLUS X COMPUTES LN(1+X) WITH RELATIVE    
     ACCURACY COMPARABLE WITH THE MACHINE ACCURACY.       
1SECTION : 6.4.3             (DECEMBER 1978)                     PAGE 2       
        
        
        
     AS IS WELL KNOWN, FOR SMALL ABS(X) RELATIVE ACCURACY IS LOST   
     WHEN COMPUTING LN(1+X) BY USING THE STANDARD FUNCTION LN.      
     THE PROCEDURE IS USED IN THE PROCEDURES ARCSINH AND ARCTANH,   
     SECTION 6.4.2. 
        
        
 REFERENCES : HART, J.F. CS. (1968), COMPUTER APPROXIMATIONS,       
              WILEY, NEW YORK.
        
        
 EXAMPLE OF USE : 
        
     WE COMPUTE LN(EXP(X)) FOR SMALL POSITIVE X. IN ORDER TO
     PRESERVE RELATIVE ACCURACY WE WRITE
        
     LN ( EXP(X) )   =    LN ( 1+ EXP(X)-1 )    
                     =    LN ( 1+ 2* EXP(X/2)* SINH(X/2) ). 
        
     THE FOLOWING PROGRAM   
        
     "BEGIN" "REAL" X,Y;    
         "REAL" "PROCEDURE" SINH(X)          ; "CODE" 35111;
         "REAL" "PROCEDURE" LOG ONE PLUS X(X); "CODE" 35130;
         "FOR" X:= "-1, "-10, "-50, "-100, "-250 "DO"     
         "BEGIN" Y:= LOG ONE PLUS X( 2*EXP(X/2)*SINH(X/2) );
             OUTPUT(61,"("N,/")",Y)   
         "END"    
     "END";       
        
     PRINTS THE FOLOWING RESULTS :    
        
     +1.0000000000000"-001  
     +1.0000000000000"-010  
     +1.0000000000000"-050  
     +1.0000000000000"-100  
     +1.0000000000000"-250  
        
        
 SOURCE TEXT :    
 "CODE" 35130;    
 "REAL" "PROCEDURE" LOG ONE PLUS X(X); "VALUE" X; "REAL" X; 
 "COMMENT" COMPUTES LN(1+X) FOR X > -1; 
 "IF" X = 0 "THEN" LOG ONE PLUS X:= 0 "ELSE"    
 "IF" X < -0.2928 "OR" X > 0.4142 "THEN" LOG ONE PLUS X:= LN(1+X) "ELSE"      
 "BEGIN" "REAL" Y,Z;
     Z:= X/(X+2); Y:= Z*Z;  
     LOG ONE PLUS X:= Z*(2+ Y*
     ( .66666 66666 63366 + Y*
     ( .40000 00012 06045 + Y*
     ( .28571 40915 90488 + Y*
     ( .22223 82333 2791  + Y*
     ( .18111 36267 967   + Y*
       .16948 21248 8)))))) 
 "END" LOG ONE PLUS X;      
        "EOP"     
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR: M. BAKKER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
 
 
 RECEIVED: 791231.
 
 
 BRIEF DESCRIPTION:
 
       THE PROCEDURE NONLIN FEMLAGSKEW SOLVES A NONLINEAR TWO POINT
       BOUNDARY VALUE PROBLEM WITH SPHERICAL COORDINATES.
       IT SOLVES THE DIFFERENTIAL EQUATION
 
           (X**NC*Y')'/X**NC = F(X, Y, Y'), 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,
       BOUNDARY VALUE PROBLEMS,
       RITZ-GALERKIN METHOD,
       SPHERICAL COORDINATES,
       GLOBAL METHODS.
 
 
 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 27, 1976.
 
       [3]  BABUSKA, I.,
            NUMERICAL STABILITY IN PROBLEMS OF LINEAR ALGEBRA,
            S.I.A.M.  J. NUM. ANAL., VOL.9, P. 53-77 (1972).
 
       [4]  BAKKER, M.,
            GALERKIN METHODS IN SPHERICAL REGIONS, TO APPEAR.
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 2
 
 
 
 SUBSECTION: NONLIN FEM LAG SKEW.
 
 
 CALLING SEQUENCE:
 
       THE HEADING OF THE PROCEDURE READS:
 
       "PROCEDURE" NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
       "INTEGER" N, NC;
       "REAL" "PROCEDURE" F, FY, FZ;
       "ARRAY" X, Y, E;
       "CODE" 33314;
 
       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 SEGMENT [A,B];
 
       Y:  <ARRAY IDENTIFIER>;
           "ARRAY" Y[0:N];
           ENTRY: Y[I] (I = 0, 1, ... , N) IS AN INITIAL APPROXIMATE
           SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
 
           (3)  (Y'*X**NC)'/X**NC = F(X,Y,Y')    , 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];
 
           EXIT: Y[I] (I = 0, 1, ... , N) IS THE GALERKIN
           SOLUTION AT X[I] OF THE (3)-(4);
 
       F:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF F READS:
           "REAL" "PROCEDURE" F(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           F(X,Y,Z) IS THE RIGHT HAND SIDE OF (3);
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 3
 
 
 
       FY:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF FY READS:
           "REAL" "PROCEDURE" FY(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           FY(X,Y,Z) IS THE DERIVATIVE OF F WITH RESPECT TO Y;
 
       FZ:  <PROCEDURE IDENTIFIER>;
           THE HEADING OF FZ READS:
           "REAL" "PROCEDURE" FZ(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           FZ(X,Y,Z) IS THE DERIVATIVE OF F WITH RESPECT TO Z;
 
       E:  <ARRAY IDENTIFIER>;
           "ARRAY" E[1:6];
           E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (4);
           E[1] AND E[4] ARE NOT ALLOWED TO VANISH BOTH.
 
 
 PROCEDURES USED: DUPVEC CP 31030.
 
 
 REQUIRED CENTRAL MEMORY:
 
       FIVE AUXILIARY ARRAYS OF N REALS ARE USED.
 
 
 RUNNING TIME:
 
       LET IT BE THE NUMBER OF NEWTON ITERATIONS; THEN
       IT*N  EVALUATIONS OF F, FY, FZ ARE NEEDED;
 
 
 DATA AND RESULTS:
 
      THE FUNCTIONS F, FY AND FZ ARE REQUIRED TO BE SUFFICIENTLY
      SMOOTH IN THEIR VARIABLES ON THE INTERIOR OF EVERY SEGMENT
      <X[I],X[I+1]> (I = 0, ..., N - 1);
 
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 4
 
 
 
 METHOD AND PERFORMANCE:
 
      LET Y[0](X) BE SOME INITIAL APPROXIMATION OF Y(X); THEN
      THE NONLINEAR PROBLEM IS SOLVED BY SUCCESIVELY SOLVING
 
          - (D[K]'*X**NC)'/X**NC
 
          + FY(X,Y[K](X),Y[K]'(X))*D[K](X)
 
          + FZ(X,Y[K](X),Y[K]'(X))*D[K]'(X)
 
        =   (Y[K]'*X**NC)'/X**NC
 
          - F(X,Y[K],Y[K]'(X)), X[0] < X < X[N],
 
 
          E[1]*D[K](X[0]) + E[2]*D[K]'(X[0]) = 0;
 
          E[4]*D[K](X[N]) + E[5]*D[K]'(X[N]) = 0;
 
      WITH GALERKIN'S METHOD (SEE PREVIOUS SECTION) AND PUTTING
 
 
          Y[K+1](X) = Y[K](X) + D[K](X), K = 0,1,...
 
 
      THIS IS THE SO-CALLED NEWTON-KANTOROWITCH METHOD;
 
 EXAMPLE OF USE:
 
       WE SOLVE THE BOUNDARY VALUE PROBLEM
 
 
           (Y'*X**2)'/X**2 = EXP(Y)+EXP(Y')-EXP(1-X**2)-EXP(2*X)-6;
           0 < X < 1, Y'(0) = Y(1) = 0;
 
       FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
 
           E[2] = E[4] = 1; E[1] = E[3] = E[5] = E[6] = 0;
 
       THE ANALYTIC SOLUTION IS Y(X) = 1 - X**2; WE APPROXIMATE
       THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I/N,
       I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE THE MAXIMUM ERROR;
       THE PROGRAM READS AS FOLLOWS:
 
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 5
 
 
 
       "BEGIN" "INTEGER" N, NC;
       "FOR" NC:= 0,1,2 "DO"
       "FOR" N:= 25, 50 "DO"
       "BEGIN" "INTEGER" I;"ARRAY" X, Y[0:N], E[1:6]; "REAL" RHO, D;
 
           "REAL" "PROCEDURE" F(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           F:= EXP(Y)+EXP(Z)-EXP(1-X**2)-EXP(-2*X)-2-2*NC;
 
           "REAL" "PROCEDURE" FY(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           FY:= EXP(Y);
 
           "REAL" "PROCEDURE" FZ(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
           FZ:= EXP(Z);
 
           "PROCEDURE" NONLIN FEM LAG SKEW(X,Y,N,F,FY,FZ,NC,E);
           "CODE" 33314;
           E[2]:= E[4]:= 1; E[1]:= E[3]:= E[5]:= E[6]:= 0;
           "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
           "BEGIN" X[I]:= I/N; Y[I]: = 0 "END";
           OUTPUT(61,"("//,4B"("N = ")"ZD,4B"("NC = ")"ZD")",N,NC);
           NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
           RHO:= 0;
           "FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
           "BEGIN" D:= ABS(Y[I] - 1 + X[I]**2);
               "IF" RHO < D "THEN" RHO:= D
           "END";
           OUTPUT(61,"("24B"("MAX.ERROR= ")",D.DD"+ZD")",RHO)
       "END"
       "END"
 
       RESULTS:
 
     N = 25    NC =  0                        MAX.ERROR= 2.47" -4
 
     N = 50    NC =  0                        MAX.ERROR= 6.19" -5
 
     N = 25    NC =  1                        MAX.ERROR= 1.41" -3
 
     N = 50    NC =  1                        MAX.ERROR= 3.99" -4
 
     N = 25    NC =  2                        MAX.ERROR= 2.44" -3
 
     N = 50    NC =  2                        MAX.ERROR= 7.02" -4
 
       ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
       0.25 WHEN THE MESH SIZE IS HALVED.
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 6
 
 
 
 SOURCE TEXT(S):
0"CODE" 33314;
  "PROCEDURE" NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
  "INTEGER" N, NC;
  "REAL" "PROCEDURE" F, FY, FZ;
  "ARRAY" X, Y, E;
  "BEGIN" "INTEGER" L, L1, IT;
    "REAL" XL1, XL, H, A12, A21, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
      PLM, PRM, PL1, PL3, PL1PL2, PL1PL3, PL2PL2, PL2PL3,
      PR1PR2, PR1PR3, PR2PR3, PL1QL2, PL1QL3, PL2QL1, PL2QL2, PL2QL3,
      PL3QL1, PL3QL2, PR1QR2, PR1QR3, PR2QR1, PR2QR2, PR2QR3, PR3QR1,
      PR3QR2, H2RM, ZL1, ZL, E1, E2, E3, E4, E5, E6, EPS, RHO;
    "ARRAY" T, SUPER, SUB, CHI, GI[0:N-1], Z[0:N];
 
    "PROCEDURE" DUPVEC(L,U,S,A,B); "CODE" 31030;
 
    "PROCEDURE" ELEMENT MAT VEC EVALUATION 1;
    "BEGIN" "REAL"  XM,VL,VR,WL,WR,PR,QM,RM,FM,XL12,XL1XL,XL2,ZM,ZACCM;
      "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" XL12:= XL1*XL1/12; XL1XL:=XL1*XL/6; XL2:=XL*XL/12;
        VL:= 3*XL12 + XL1XL + XL2;
        VR:= 3*XL2 + XL1XL + XL12
      "END";
      WL:= H*VL; WR:=H*VR; PR:= VR/(VL +VR);
      XM:= XL1 + H*PR; ZM:= PR*ZL + (1 - PR)*ZL1;
      ZACCM:= (ZL - ZL1)/H ; QM:= FZ(XM,ZM,ZACCM);
      RM:= FY(XM, ZM, ZACCM); FM:= F(XM,ZM,ZACCM);
      TAU1:= WL*RM; TAU2:=WR*RM;
      B1:= WL*FM - ZACCM*(VL +VR); B2:= WR*FM + ZACCM*(VL + VR);
      A12:= - (VL + VR)/H + VL*QM + (1 - PR)*PR*RM*(WL + WR);
      A21:= - (VL + VR)/H - VR*QM + (1 - PR)*PR*RM*(WL + WR);
    "END" ELEM. M.V. EV.;
 
 
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 7
                                                                       ;
 
 
    "PROCEDURE" BOUNDARY CONDITIONS;
    "IF" L=1 "AND" E2 = 0 "THEN"
    "BEGIN" TAU1:= 1; B1:= A12:= 0 "END"
    "ELSE" "IF" L=1 "AND" E2 ^= 0 "THEN"
    "BEGIN" TAU1:= TAU1 - E1/E2
    "END" "ELSE" "IF" L=N "AND" E5 = 0 "THEN"
    "BEGIN" TAU2:= 1; B2:= A21:= 0
    "END" "ELSE" "IF" L=N "AND" E5 ^= 0 "THEN"
    "BEGIN" TAU2:= TAU2 + E4/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;
                                                               "COMMENT"
1SECTION : 5.2.1.2.1.3        (DECEMBER 1979)                     PAGE 8
                                                                       ;
 
 
    DUPVEC(0,N,0,Z,Y);
      E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
    "FOR" IT:= 1, IT + 1 "WHILE" EPS > RHO "DO"
    "BEGIN" L:= 0;XL:= X[0]; ZL:= Z[0];
      "FOR" L:= L + 1 "WHILE" L <= N "DO"
      "BEGIN" XL1:= XL; L1:= L - 1; XL:= X[L]; H:= XL - XL1;
         ZL1:= ZL; ZL:= Z[L];
         ELEMENT MAT VEC EVALUATION 1;
         "IF" L=1 "OR" L=N "THEN" BOUNDARY CONDITIONS;
         FORWARD BABUSKA
      "END";
      BACKWARD BABUSKA;
      EPS:= 0; RHO:= 1;
      "FOR" L:= 0 "STEP" 1 "UNTIL" N "DO"
      "BEGIN" RHO:= RHO + ABS(Z[L]);
        EPS:= EPS + ABS(Y[L]); Z[L]:= Z[L] - Y[L]
      "END";
      RHO:= "-14*RHO
    "END";
      DUPVEC(0,N,0,Y,Z)
  "END" NONLIN FEM LAG SKEW;
         "EOP"