PROGRAM KRYLOW(INPUT,OUTPUT);                     
(*$R+                           *)                
(* EW BESTIMMUNG: NUMERIK ZEHN3 *)                
CONST N=4;EPS=5.0E-4;                             
      LEN=16; (* FORMATIERUNG   *)                
                                                  
TYPE VEKTOR=ARRAY[1..N] OF REAL;                  
     ITERATION=0..N;                              
     MATRIX=ARRAY[1..N] OF VEKTOR;                
     ITVEKTOR=ARRAY[ITERATION] OF VEKTOR;         
     ITSKALAR=ARRAY[ITERATION] OF REAL;           
                                                  
VAR I,J,K,M,IV,IW:INTEGER;                        
    S,MX,P,X,D:REAL;                              
    LINABH:BOOLEAN;                               
    B:MATRIX;         (* AUSGANGSMATRIX     *)    
    A,V,W:ITVEKTOR;   (* ITERIERTE VEKTOREN *)    
    VN:ITSKALAR;      (* VEKTORNORM         *)    
    AL,A1,A2,G:                                   
             ITSKALAR;(* KOEFFIZIENTEN      *)    
    Z:VEKTOR;                                     
(*$"KRYLOW1.INC"*)                                
BEGIN                                             
 (* MATRIX VORBELEGEN *)                          
 FOR I:= 1 TO N DO                                
  BEGIN                                           
   FOR J:= 1 TO N DO READ(B[I,J]);                
   READLN                                         
  END;                                            
 (* MAXIMALE ZEILENSUMME VON B *)                 
 MX:=0.0;                                         
 FOR J:= 1 TO N DO                                
  BEGIN S:=0.0;                                   
   FOR K:= 1 TO N DO                              
    S:=S+ABS(B[J,K]);                             
   MX:=MAX(MX,S)                                  
  END;                                            
 FOR IV:=1 TO N DO                                
  BEGIN (* FUER ALLE BASISVEKTOREN *)             
   FOR J:=1 TO N DO                               
    A[0,J]:=ORD(IV=J);                            
                                                  
   V[0]:=A[0];W[0]:=A[0];VN[0]:=1.0;              
   (* STARTVEKTOREN SIND BEREITS NORMIERT *)      
                                                  
   J:=1; (* ITERATION FUER A[J] *)                
   REPEAT                                         
    FOR I:=1 TO N DO                              
     BEGIN S:=0.0;                                
      FOR K:=1 TO N DO                            
       S:=S+B[I,K]*A[J-1,K];                      
      A[J,I]:=S;                                  
     END;                                         
     (* ORTHOGONALISIERE A --> V *)               
     V[J]:=A[J];                                  
     FOR K:=0 TO J-1 DO                           
      BEGIN S:=0.0;                               
       FOR I:=1 TO N DO                           
        S:=S+A[J,I]*W[K,I];                       
       FOR I:=1 TO N DO                           
        V[J,I]:=V[J,I]-S*W[K,I];                  
      END;                                        
     (* BERECHNE İİVİİ            *)              
     S:=0.0;                                      
     FOR I:=1 TO N DO                             
      S:=S+SQR(V[J,I]);                           
     S:=SQRT(S);VN[J]:=S;                         
     LINABH:= S<EPS;                              
     IF NOT LINABH THEN                           
       (* NORMIERE V --> W          *)            
      BEGIN                                       
       FOR I:=1 TO N DO                           
        W[J,I]:=V[J,I]/S;                         
       J:=J+1                                     
      END;                                        
   UNTIL LINABH;                                  
   (* BERECHNE DIE KOEFFIZIENTEN DES*)            
   (* MINIMALPOLYNOMS DER VEKTOREN  *)            
   M:=J;                                          
   WRITELN("DIE VEKTOREN LAUTEN:");               
   FOR J:=0 TO M DO                               
    BEGIN                                         
     WRITE("A[",J,"]=");                          
     FOR I:=1 TO N DO WRITE(A[J,I]:LEN);          
     WRITELN;WRITE("V[",J,"]=");                  
     FOR I:=1 TO N DO WRITE(V[J,I]:LEN);          
     WRITELN;WAIT                                 
    END;                                          
   AL[M]:=1.0;                                    
   J:=M-1;AL[J]:=0.0;                             
   FOR I:=1 TO N DO                               
    AL[J]:=AL[J]+A[M,I]*V[J,I];                   
   AL[J]:=-(AL[J]/SQR(VN[J]));                    
   (* REKURSIV WEITERRECHNEN: *)                  
   FOR J:= M-2 DOWNTO 0 DO                        
    BEGIN S:=0.0;                                 
     FOR I:=1 TO N DO                             
      S:=S+A[M,I]*V[J,I];                         
     FOR K:=J+1 TO M-1 DO                         
      BEGIN P:=0.0;                               
       FOR I:= 1 TO N DO                          
        P:=P+A[K,I]*V[J,I];                       
       S:=S+AL[K]*P                               
      END;                                        
     AL[J]:=-(S/SQR(VN[J]));                      
    END;                                          
   WRITELN("DIE KOEFFIZIENTEN (0..M) LAUTEN:");   
   FOR I:=0 TO M DO WRITE(AL[I]:LEN);             
   WRITELN;                                       
   A1:=AL;                                        
   (* BERECHNE NULLSTELLEN NACH NEWTON-VERFAHREN*)
   FOR IW:=1 TO M DO                              
    BEGIN                                         
     X:=0.0; (* ANFANGSWERT *)                    
     REPEAT                                       
      A2[M]:=A1[M];                               
      FOR K:=M-1 DOWNTO 0 DO                      
       A2[K]:=A1[K]+X*A2[K+1];                    
      P:=A2[M];(*1.ABLEITUNG MIT HORNER-SCHEMA  *)
      FOR K:=M-1 DOWNTO 1 DO                      
       P:=A2[K]+K*P;                              
       D:=A2[0]/P;X:=X-D;                         
     UNTIL ABS(D)<=X*1E-6;                        
     WRITELN("EIGENWERT=",X:LEN);                 
     (* LINEARFAKTOR ABSPALTEN (HORNER-SCHEMA)  *)
     G[M]:=AL[M];                                 
     FOR K:=M-1 DOWNTO 0 DO                       
      G[K]:=AL[K]+X*G[K+1];                       
     FOR I:= 0 TO M-1 DO                          
      A1[I]:=A2[I+1];                             
     (* G ENTHAELT KOEFFIZIENTEN DES TEILERS    *)
     P:=0.0; (* EIGENVEKTOREN BERECHNEN *)        
     FOR I:= 1 TO N DO                            
      BEGIN S:=0.0;                               
       FOR K:=0 TO M-1 DO                         
        S:=S+A[K,I]*G[K+1];                       
       Z[I]:=S;P:=P+S*S;                          
      END;                                        
     P:=SQRT(P);                                  
     FOR I:= 1 TO N DO                            
      Z[I]:=Z[I]/P;                               
     WRITE("EIGENVEKTOR:");                       
     FOR I:= 1 TO N DO WRITE(Z[I]:LEN);WRITELN;   
     (* DEFEKT PRUEFEN *)                         
     WRITE("DEFEKT:");                            
     D:=0.0;                                      
     FOR I:=1 TO N DO                             
      BEGIN S:=-X*Z[I];                           
       FOR K:=1 TO N DO                           
        S:=S+B[I,K]*Z[K];                         
       D:=D+S*S                                   
      END;                                        
     D:=SQRT(D);WRITELN(D);                       
    END;                                          
   WRITELN;WRITELN;WAIT                           
  END;(* FUER ALLE BASISVEKTOREN (IV)*)           
END.