Mercurial > matrix-functions
comparison mftoolbox/arnoldi.m @ 0:8f23314345f4 draft
Create local repository for matrix toolboxes. Step #0 done.
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Wed, 06 May 2015 14:56:53 +0200 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:8f23314345f4 |
|---|---|
| 1 function [Q,H] = arnoldi(A,q1,m) | |
| 2 %ARNOLDI Arnoldi iteration | |
| 3 % [Q,H] = ARNOLDI(A,q1,M) carries out M iterations of the | |
| 4 % Arnoldi iteration with N-by-N matrix A and starting vector q1 | |
| 5 % (which need not have unit 2-norm). For M < N it produces | |
| 6 % an N-by-(M+1) matrix Q with orthonormal columns and an | |
| 7 % (M+1)-by-M upper Hessenberg matrix H such that | |
| 8 % A*Q(:,1:M) = Q(:,1:M)*H(1:M,1:M) + H(M+1,M)*Q(:,M+1)*E_M', | |
| 9 % where E_M is the M'th column of the M-by-M identity matrix. | |
| 10 | |
| 11 n = length(A); | |
| 12 if nargin < 3, m = n; end | |
| 13 q1 = q1/norm(q1); | |
| 14 Q = zeros(n,m); Q(:,1) = q1; | |
| 15 H = zeros(min(m+1,m),n); | |
| 16 | |
| 17 for k=1:m | |
| 18 z = A*Q(:,k); | |
| 19 for i=1:k | |
| 20 H(i,k) = Q(:,i)'*z; | |
| 21 z = z - H(i,k)*Q(:,i); | |
| 22 end | |
| 23 if k < n | |
| 24 H(k+1,k) = norm(z); | |
| 25 if H(k+1,k) == 0, return, end | |
| 26 Q(:,k+1) = z/H(k+1,k); | |
| 27 end | |
| 28 end |