## how to avoid floating point error and figure out a way to break out of this loop such that it generates correct answer

### Nitish Reddy Kotkur (view profile)

on 21 Oct 2019 at 23:41
Latest activity Answered by Christine Tobler

### Christine Tobler (view profile)

on 28 Oct 2019 at 16:44
function [] = lanczos(A, m)
this is a function whick takes A as input which is a matrix and reduces it to smaller size. And i have got a for loop
for j=2:inf
w = A*V(:,j) - beta(j)*V(:,j-1);
alpha(j) = w'*V(:,j);
w = w - alpha(j)*V(:,j);
beta(j+1) = norm(w,2);
loopcnt = loopcnt + 1;
if abs(beta(j+1)) = 0
break
end
end
now i need the control to come out of loop when beta(j+1) is equal to zero ,but its looping continuously may be because of some floating point error.so i tried something like this
if abs(beta(j+1)) < 0.0001 should come out of loop which worked fine for smaller matrix sizes.
but when matrix size got bigger even its not working and loop is running continuosly can some one suggest me a way to avoid this problem and run the loop properly and break it when beta(j+1) become zero

Show 1 older comment
Nitish Reddy Kotkur

### Nitish Reddy Kotkur (view profile)

on 22 Oct 2019 at 1:05
function [] = lanczos(A, m)
A = readmatrix('output1.txt','Whitespace',' []'); %output1.txt contains matrix
[n,k] = size(A);
V = zeros(k,k);
V(:,2) = 0;
V(1,2)=1;
V(:,2)=V(:,2)';
V(:,2)=V(:,2)/norm(V(:,2),2);
beta(2)=0;
loopcnt = 0;
for j=2:inf
w = A*V(:,j) - beta(j)*V(:,j-1);
alpha(j) = w'*V(:,j);
w = w - alpha(j)*V(:,j);
beta(j+1) = norm(w,2);
loopcnt = loopcnt + 1;
if abs(beta(j+1)) < 0.0001
break
end
V(:,j+1) = w/beta(j+1);
end
m=loopcnt-1;
T=sparse(m+1,m+1);
for i=2:m+1
T(i-1,i-1)=alpha(i);
T(i-1,i)=beta(i+1);
T(i,i-1)=beta(i+1);
end
disp(['number of loops: ' num2str(loopcnt)]);
T(m+1,m+1)=alpha(m+2);
V = V(:,2:loopcnt+1)
disp(['approximating eigenvalues are: ', num2str(eigs(T)')]);
T
j
end
This is my entire code,and i have highlighted the loop part. could you help me with it
David Hill

### David Hill (view profile)

on 22 Oct 2019 at 1:31
j=2;
while abs(beta(j))>=0.0001
w = A*V(:,j) - beta(j)*V(:,j-1);
alpha(j) = w'*V(:,j);
w = w - alpha(j)*V(:,j);
beta(j+1) = norm(w,2);
j=j+1;
end
If you are sure abs(beta(j)) will approach zero, you could troubleshoot by printing beta(j+1) during each loop by removing the semicolon. You should also preallocate alpha and beta if you have some idea how many cycles the loop will make (it will dramatically increase the speed).
Nitish Reddy Kotkur

### Nitish Reddy Kotkur (view profile)

on 22 Oct 2019 at 2:50
well its not working the problem is figuring out perfect suitable tolerance,which i am unable to find