Jeremy Knight Math 639 Algorithm 9.4 Wielandt Deflation
<Text-field style="Heading 1" layout="Heading 1">Wielandt(n,A,lambda,v,x,TOL,N,display)</Text-field> Wielandt:=proc(n,A,lambda,v,x,TOL,N,display) #Wieland Deflation algorithm to approximate the second most dominant eigenvalue and an associated eigenvector of the nxn matrix A given an approximation lambda to the dominant eigenvalue, an approximation v to a corresponding eigenvector, and a vector x in R^(n-1), tolerance TOL, max iterations N local B,i,j,k,u,w,w2,xb,max,index,mu,muhat, temp, PowerMethod; option trace; ##subroutine section PowerMethod:=proc(n,A,x0,TOL,N,display) #n=dimension, A=Array, x=vector, TOL=tolerance, N=max iteraions #This version has been modified to use Aitkens delta^2 procedure local i,j,k,p,x,xp,y,mu,yp, w,EuclidNorm,ERR,temp,mu0,mu1,muhat: #option trace; k:=1: mu0:=0: mu1:=0: x:=x0: if display=1 then print("iterations, eigenvalue, eigenvector, muhat"); end if; #subroutine for finding the Euclidean Norm and #returning the norm and the index of the norm. EuclidNorm:=proc(n,x) #n=vector size, x= 1 by n array local euclid_norm,j,p: euclid_norm:=abs(x[1]): p:=1; for j from 2 to n do if abs(x[j])>euclid_norm then euclid_norm:=abs(x[j]): p:=j: end if: end do: return(euclid_norm,p); end proc: #find the smallest integer p with 1<=p<=n and abs(x[p])=norm(X) xp,p:=EuclidNorm(n,x); #normalize x x:=x/x[p]: while k<=N do #Calculate next eigenvector approximate y:=Array(1..n): for i from 1 to n do for j from 1 to n do y[i]:=y[i]+A[i,j]*x[j]: end do: end do: mu:=y[p]: #use Aitken's delta^2 procedure to speed convergence muhat:=mu0-((mu1-mu0)^2/(mu-2*mu1+mu0)): #find the smallest integer p with 1<=p<=n and abs(y[p])=norm(y) yp,p:=EuclidNorm(n,y): if y[p]=0 then print("Eigenvector",evalf(x)); print("A has the eigenvalue 0, select new vector x and restart"); break; end if: #Modify with abs(x)-abs(y/y[p] to find an accurate error measurement when signs alternate in eigenvector. ERR,temp:=EuclidNorm(n,abs(x)-abs(y/y[p])): x:=y/y[p]; #check for tolerance reached if ERR<TOL and k>=4 then print("The procedure was successful"); print("eigenvalue approximate (muhat)=",evalf(muhat)); print("eigenvector approximate=",evalf(x)); break; end if; if display=1 then print(k, evalf(mu), evalf(x), evalf(muhat)); end if; k:=k+1: mu0:=mu1: mu1:=mu: end do; if k=N then print("The maximum number of iterations exceeded"); end if: if ERR>TOL then print("tolerance not reached, Method Fails"); muhat:=0: x:=Array(1..n): end if; return(muhat,x): end proc: ##End subroutine section max:=v[1]: i:=1; for index from 2 to n do if v[index]>max then max:=v[index]: i:=index: end if: end do: print("i",i); B:=Array(1..n-1,1..n-1): if i<>1 and i<>n then for k from i to n-1 do for j from 1 to i-1 do B[k,j]:=A[k+1,j]-v[k+1]/v[i]*A[i,j]: B[j,k]:=A[j,k+1]-v[j]/v[i]*A[i,k+1]: end do: end do: end if: if i<>n then for k from i to n-1 do for j from i to n-1 do B[k,j]:=A[k+1,j+1]-v[k+1]/v[i]*A[i,j+1] end do: end do: end if: muhat,xb:=PowerMethod(n-1,B,x,TOL,N,display): if muhat=0 and xb=0 then print("Method Failed"): stop: else w2:=xb: mu:=muhat: end if: w:=Array(1..n): if i<>1 then for k from 1 to i-1 do w[k]:=w2[k]: end do: end if: if i<>n then for k from i+1 to n do w[k]:=w2[k-1]: end do: end if: #Compute the eigenvector using Eq.9.6 temp:=0: for j from 1 to n do temp:=temp+A[i,j]*w[j]: end do: u:=Array(1..n): for k from 1 to n do u[k]:=(mu-lambda)*w[k]+temp*v[k]/v[i]: end do: print("eigenvalue =", mu); print("eigenvector = ",u); end proc:
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= Example: We will use the example matrix from the InversePowerMethod algorithm which gave us the approximation: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbXNHRiQ2I1ErZWlnZW52YWx1ZUYnLUkjbW9HRiQ2LVEiLEYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdRJXRydWVGJy8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1JI21uR0YkNiRRLDYuMDAwMDAxNzAxRidGM0Yz LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbXNHRiQ2I1EsZWlnZW52ZWN0b3JGJy1JI21vR0YkNi1RIixGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGOC8lKnN5bW1ldHJpY0dGOC8lKGxhcmdlb3BHRjgvJS5tb3ZhYmxlbGltaXRzR0Y4LyUnYWNjZW50R0Y4LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictSShtYWN0aW9uR0YkNiUtSShtZmVuY2VkR0YkNigtSSdtdGFibGVHRiQ2NS1JJG10ckdGJDYoLUkkbXRkR0YkNigtSSNtbkdGJDYkUSwxLjAwMDAwMDAwMEYnRjMvJSlyb3dhbGlnbkdRIUYnLyUsY29sdW1uYWxpZ25HRltvLyUrZ3JvdXBhbGlnbkdGW28vJShyb3dzcGFuR1EiMUYnLyUrY29sdW1uc3BhbkdGYm8tRlk2KC1GZm42JFEtMC43MTQyODU4MzUyRidGM0ZpbkZcb0Zeb0Zgb0Zjby1GWTYoLUYjNiUtRjA2LVEqJnVtaW51czA7RidGM0Y2L0Y6RjhGPEY+RkBGQkZEL0ZHUSwwLjIyMjIyMjJlbUYnL0ZKRmNwLUZmbjYkUS0wLjI0OTk5OTU0OTlGJ0YzRjNGaW5GXG9GXm9GYG9GY29GaW5GXG9GXm8vJSZhbGlnbkdRJWF4aXNGJy9Gam5RKWJhc2VsaW5lRicvRl1vUSdjZW50ZXJGJy9GX29RJ3xmcmxlZnR8aHJGJy8lL2FsaWdubWVudHNjb3BlR0Y7LyUsY29sdW1ud2lkdGhHUSVhdXRvRicvJSZ3aWR0aEdGZXEvJStyb3dzcGFjaW5nR1EmMS4wZXhGJy8lLmNvbHVtbnNwYWNpbmdHUSYwLjhlbUYnLyUpcm93bGluZXNHUSVub25lRicvJSxjb2x1bW5saW5lc0dGYHIvJSZmcmFtZUdGYHIvJS1mcmFtZXNwYWNpbmdHUSwwLjRlbX4wLjVleEYnLyUqZXF1YWxyb3dzR0Y4LyUtZXF1YWxjb2x1bW5zR0Y4LyUtZGlzcGxheXN0eWxlR0Y4LyUlc2lkZUdRJnJpZ2h0RicvJTBtaW5sYWJlbHNwYWNpbmdHRl1yRjMvSSttc2VtYW50aWNzR0YkUSdWZWN0b3JGJy8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Zjcy8lK2FjdGlvbnR5cGVHUS5ydGFibGVhZGRyZXNzRicvJSlydGFibGVpZEdRKjI1MTMwOTgyNEYnLyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GJy8lKXJlYWRvbmx5R0Y4RjM= A:=Array([[-4.0,14,0],[-5.0,13,0],[-1,0,2]]); x0:=Array(1..2,[1.0,1.0]); q:=19.0/3; v:=Array([1.000000000,0.7142858352,-.2499995499]); LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKm8tMl4j LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKndHMl4j JCIrTExMTGohIio= LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKmtkMl4j Wielandt(3,A,6.00001701,v,x0,0.00001,30,1); {--> enter Wielandt, args = 3, Array(1..3, 1..3, {(1, 1) = -4.0, (1, 2) = 14, (1, 3) = 0, (2, 1) = -5.0, (2, 2) = 13, (2, 3) = 0, (3, 1) = -1, (3, 2) = 0, (3, 3) = 2}), 6.00001701, Array(1..3, {(1) = 1.000000000, (2) = .7142858352, (3) = -.2499995499}), Array(1..2, {(1) = 1.0, (2) = 1.0}), 0.1e-4, 30, 1 Zio2KEkibkc2IkkiQUdGJUkjeDBHRiVJJFRPTEdGJUkiTkdGJUkoZGlzcGxheUdGJTYySSJpR0YlSSJqR0YlSSJrR0YlSSJwR0YlSSJ4R0YlSSN4cEdGJUkieUdGJUkjbXVHRiVJI3lwR0YlSSJ3R0YlSStFdWNsaWROb3JtR0YlSSRFUlJHRiVJJXRlbXBHRiVJJG11MEdGJUkkbXUxR0YlSSZtdWhhdEdGJUYlRiVDLj44JiIiIj44MSIiIT44MkZCPjgoOSZAJC85KUY/LUkmcHJpbnRHJSpwcm90ZWN0ZWRHNiNRS2l0ZXJhdGlvbnMsfmVpZ2VudmFsdWUsfmVpZ2VudmVjdG9yLH5tdWhhdEYlPjguZio2JEYkRjA2JUksZXVjbGlkX25vcm1HRiVGLUYvRiVGJUMmPjgkLUkkYWJzR0ZNNiMmOSU2I0Y/Rj0/KDglIiIjRj85JEkldHJ1ZUdGTUAkMkZYLUZaNiMmRmduNiNGam5DJD5GWEZgbz5GPkZqbk82JEZYRj5GJUYlRiU+NiQ4KTgnLUZRNiRGXG9GRj5GRiomRkZGPyZGRjYjRlxwISIiPyhGJUY/Rj9GJTFGPjkoQy8+OCotSSZBcnJheUdGTTYjO0Y/RlxvPyhGWEY/Rj9GXG9GXW8/KEZqbkY/Rj9GXG9GXW8+JkZpcDYjRlgsJkZhcUY/KiYmRmduNiRGWEZqbkY/JkZGRmNvRj9GPz44KyZGaXBGYnA+ODMsJkZBRj8qJiwmRkRGP0ZBRmNwRltvLChGaXFGP0ZEISIjRkFGP0ZjcEZjcD42JDgsRlxwLUZRNiRGXG9GaXBAJC9GanFGQkMlLUZMNiRRLEVpZ2VudmVjdG9yRiUtSSZldmFsZkdGTTYjRkYtRkw2I1FYQX5oYXN+dGhlfmVpZ2VudmFsdWV+MCx+c2VsZWN0fm5ld352ZWN0b3J+eH5hbmR+cmVzdGFydEYlWz42JDgvODAtRlE2JEZcbywmLUZaRl9zRj8tRlo2IyomRmlwRj9GanFGY3BGY3A+RkZGXnRAJDMyRmZzOScxIiIlRj5DJi1GTDYjUT1UaGV+cHJvY2VkdXJlfndhc35zdWNjZXNzZnVsRiUtRkw2JFFAZWlnZW52YWx1ZX5hcHByb3hpbWF0ZX4obXVoYXQpPUYlLUZeczYjRlxyLUZMNiRROWVpZ2VudmVjdG9yfmFwcHJveGltYXRlPUYlRl1zRmNzQCRGSS1GTDYmRj4tRl5zNiNGaXFGXXNGXXU+Rj4sJkY+Rj9GP0Y/PkZBRkQ+RkRGaXFAJC9GPkZmcC1GTDYjUUpUaGV+bWF4aW11bX5udW1iZXJ+b2Z+aXRlcmF0aW9uc35leGNlZWRlZEYlQCQyRmN0RmZzQyUtRkw2I1FEdG9sZXJhbmNlfm5vdH5yZWFjaGVkLH5NZXRob2R+RmFpbHNGJT5GXHJGQj5GRkZqcE82JEZcckZGRiVGJUYl JCIrKysrKzUhIio= IiIi NiRRImk2IiIiIg== LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKk9DM14j JCIqSikqKioqSCEiKQ== JCIiIUYj JCIrKnAkKioqXCQhIio= JCIiIyIiIQ== UUtpdGVyYXRpb25zLH5laWdlbnZhbHVlLH5laWdlbnZlY3Rvcix+bXVoYXQ2Ig== NiYiIiIkIis1JCkqKioqSCEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipDKzVeIyQiIiFGMA== NiYiIiMkIislZSkzNFIhIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK1o2XiMkIitbXk0vVkYm NiYiIiQkIitTWDZsTSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipfRzheIyQiK2lNejVPRiY= NiYiIiUkIisjemElb0shIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqO3k5XiMkIisuWTI3SkYm NiYiIiYkIistJnBVOyQhIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqNzk7XiMkIitQLipvLyRGJg== NiYiIickIitOciNRNSQhIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqISkpeTZEJCIrKjM1Li0kRiY= NiYiIigkIitfPyFwMSQhIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqQ0k+XiMkIitXXSMqM0lGJg== NiYiIikkIitPeWlWSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiozXj9eIyQiK3JoJVIrJEYm NiYiIiokIitMeG1HSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiopKVJBXiMkIismM1w8KyRGJg== NiYiIzUkIitYLiQqPUkhIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqI3pQN0QkIitrY3grSUYm NiYiIzYkIitDMGE3SSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipfPUReIyQiK2lPTStJRiY= NiYiIzckIishKVxLM0khIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqXyNwN0QkIit2PDorSUYm NiYiIzgkIitoU2AwSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIio7VkdeIyQiKyhcbSsrJEYm NiYiIzkkIis7P28uSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiovPkleIyQiK0cnRysrJEYm NiYiIzokIisuNlgtSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiozZ0peIyQiK3M8LCtJRiY= NiYiIzskIityQGosSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipTV0xeIyQiKzFWKytJRiY= NiYiIzwkIitlcDMsSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipnd01eIyQiK2s0KytJRiY= NiYiIz0kIis4UXMrSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipPTU9eIyQiKysmKioqKipIRiY= NiYiIz4kIitpPVsrSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiovblBeIyQiK0spKSoqKipIRiY= NiYiIz8kIitFMUsrSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipvQlJeIyQiK1UmKSoqKipIRiY= NiYiI0AkIitsSkArSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipDalNeIyQiK0AlKSoqKipIRiY= NiYiI0EkIitPOjkrSSEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiprP1VeIyQiK1skKSoqKipIRiY= UT1UaGV+cHJvY2VkdXJlfndhc35zdWNjZXNzZnVsNiI= NiRRQGVpZ2VudmFsdWV+YXBwcm94aW1hdGV+KG11aGF0KT02IiQiK1EkKSoqKipIISIq NiRROWVpZ2VudmVjdG9yfmFwcHJveGltYXRlPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKmspUjlE NiQkIitRJCkqKioqSCEiKi1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIipXdVZeIw== LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKld1Vl4j JCIrUSQpKioqKkghIio= LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKiUpcFdeIw== JCIrS11LZEchIzU= JCIrKioqKioqKioqKiEjNQ== IiIh JCIiIUYj JCIrWF1EK1MhIio= JCIrWF1EK1MhIio= LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKjcyWF4j JCIrWF1EK1MhIio= JCIrckU3Kz8hIio= JCErTDEzK1MhIio= NiRRLWVpZ2VudmFsdWV+PTYiJCIrUSQpKioqKkghIio= NiRRL2VpZ2VudmVjdG9yfj1+NiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEdGJCIqNzJYXiM= <-- exit Wielandt (now at top level) = }
<Text-field style="Heading 1" layout="Heading 1">InversePowerMethod(n,A,q,x0,TOL,N,display)</Text-field> InversePowerMethod:=proc(n,A,q,x0,TOL,N,display) #n=dimension,A=Matrix,x=vector,TOL=tolerance,N=max iterations local i,j,k,p,x,y,yp,ERR,A2,mu,mu0,mu1, muhat, temp,MVProduct,EuclidNorm,IdentityMatrix,euclid_norm,gauss: #option trace; ##Load Subroutines MVProduct:=proc(n,A,x) #subroutine for calculating y=Ax for Matrix A and vector x local i,j,p; p:=Array(1..n); for i from 1 to n do for j from 1 to n do p(i):=p(i)+A(i,j)*x(j): end do: end do: return(p): end proc: EuclidNorm:=proc(n,x) #subroutine to compute euclidean norm and identify the index of the norm #n=vector size, x= 1 by n array local euclid_norm,j,p: euclid_norm:=abs(x[1]): p:=1; for j from 2 to n do if abs(x[j])>euclid_norm then euclid_norm:=abs(x[j]): p:=j: end if: end do: return(euclid_norm,p); end proc: IdentityMatrix:=proc(n) local i,j,M; M:=Array(1..n,1..n): for i from 1 to n do M(i,i):=1: end do: return(M); end proc: gauss:=proc(aa, bb, n) local i, j, k, p, m, tmp, maxval, tmpval, s, x, nrow, a, b; a:=Array(1..n,1..n); b:=Array(1..n); x:=Array(1..n); nrow:=Array(1..n); # copy matrix into a local variable so original # matrix is not changed by the procedure for i from 1 to n do for j from 1 to n do a[i,j]:= aa[i,j]; end do; end do; # copy vector for same reason for i from 1 to n do b[i]:= bb[i]; end do; # initialize nrow for i from 1 to n do nrow[i]:= i; end do; s:=Array(1..n); # initialize scale factors for i from 1 to n do maxval:=abs(a[i,1]); for j from 2 to n do tmpval:=abs(a[i,j]); if ( tmpval > maxval ) then maxval:=tmpval; end if; end do; s[i]:=maxval; end do; for i from 1 to n-1 do # find pivot row maxval:= abs(a[nrow[i],i])/s[nrow[i]]; p:=i; for j from i+1 to n do tmpval:= abs(a[nrow[j],i])/s[nrow[j]]; if (tmpval > maxval) then p:=j; maxval:=tmpval; end if; end do; if ( maxval = 0 ) then return(0); end if; # switch rows if ( p <> i ) then tmp:=nrow[p]; nrow[p]:=nrow[i]; nrow[i]:=tmp; end if; # zero out column i for k from i+1 to n do m:=a[nrow[k],i]/a[nrow[i],i]; a[nrow[k],i]:=0; for j from i+1 to n do a[nrow[k],j]:=a[nrow[k],j]-m*a[nrow[i],j]; end do; b[nrow[k]]:=b[nrow[k]]-m*b[nrow[i]]; end do; if ( a[nrow[n],n] = 0 ) then return(0); end if; end do; # backward substitution x[n]:= b[nrow[n]]/a[nrow[n],n]; for j from n-1 to 1 by -1 do x[j]:=b[nrow[j]]; for k from n to j+1 by -1 do x[j]:=x[j]- a[nrow[j],k]*x[k]; end do; x[j]:=x[j]/a[nrow[j],j]; end do; return(eval(x)); end: ##End subroutine section k:=1: mu0:=0: mu1:=0: if display=1 then print ("k, x, mu, muhat"): end if: #find the smallest integer p with 1<=p<=n and abs(x[p])=norm(x) euclid_norm,p:=EuclidNorm(n,x0): x:=x0/x0[p]: A2:=A-q*IdentityMatrix(n): while k<=N do y:=gauss(A2,x,n): if y=0 then print("q is an eigenvalue", q): break: end if: mu:=y[p]: #use Aitken's delta^2 procedure to speed convergence muhat:=mu0-((mu1-mu0)^2/(mu-2*mu1+mu0)): #find the smallest integer p with 1<=p<=n and abs(y[p])=norm(y); temp,p:=EuclidNorm(n,y): ERR,temp:=EuclidNorm(n,x -(y/y[p])): x:=y/y[p]: if ERR<TOL and k>=4 then mu:=1/mu+q: print("tolerance reached"); print("eigenvalue",mu); print("eigenvector",x); break: end if; if display=1 then print(k-1, x, (1.0/mu+q), (1.0/muhat+q)): end if; k:=k+1; mu0:=mu1: mu1:=mu: end do: if k=N then print("Maximum number of iterations exceeded"): end if: end proc: