function h=arnoldi(A,ns)
%
% Dieser Code demonstriert den Verlust der Orthogonalität, wenn
% der Arnoldi-Prozess ohne Reorthogonalisierung benutzt wird
%
% Input:
%
% A = (volle oder sparse) Matrix der Dimension n mal n
% ns = Anzahl der Arnoldi-Schritte
%
% Output:
%
% h = berechnete Hessenbergmatrix der
%     Dimension ns x ns

n=size(A,1);
  
qhold = randn(n,1);  % Erzeugt zufälligen Anfangsvektor
qhold = qhold/norm(qhold);

for versuch = 1:2
   if versuch == 1
      disp(' ')
      disp('Arnoldi ohne Reorthogonalisierung:')
   else
      disp('Arnoldi mit Reorthogonalisierung:')
   end
   q = zeros(n,ns);   
   q(:,1)=qhold;
   for k = 1:ns
      p = A*q(:,k);
      for j = 1:k
         h(j,k) = q(:,j)'*p;       
	 p = p - q(:,j)*h(j,k);   
      end
      
      if versuch == 2    % Reorthogonalisiere
         for j = 1:k
	    d = q(:,j)'*p;
	    p = p - q(:,j)*d;
	    h(j,k) = h(j,k) + d;
	 end
      end
      
      if k == ns break; end
      q(:,k+1) = p;
      h(k+1,k) = norm(q(:,k+1));
      if h(k+1,k) == 0
         disp('Invarianten Unterraum gefunden.')
	 space_dimension = k
	 return
      else
         q(:,k+1) = q(:,k+1)/h(k+1,k);
      end
   end
   if versuch == 1
      firstcol = q(:,2:ns)'*q(:,1)      % Drucke einige Ergebnisse aus
   end
   orthodefect = norm(eye(ns)-q(:,1:ns)'*q(:,1:ns))
end      	 	    
	 
