麻省理工ComputationalMethod

				function [Q,H]=arnoldi(A,b,k)
				
				m=size(A,1);
				Q=zeros(m,k);
				H=zeros(k+1,k);
				Q(:,1)=b/norm(b);
				for n=1:k
				  v=A*Q(:,n);
				  for j=1:n
				    H(j,n)=Q(:,j)'*v;
				    v=v-H(j,n)*Q(:,j);
				  end
				  H(n+1,n)=norm(v);
				  Q(:,n+1)=v/H(n+1,n);
				end