Jeremy Knight Math 639 2/22/11 Householder's Method To obtain a symmetric tridiagonal matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNictRi82JVEibkYnRjJGNS1JI21vR0YkNi1RIitGJy9GNlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkovJSlzdHJldGNoeUdGSi8lKnN5bW1ldHJpY0dGSi8lKGxhcmdlb3BHRkovJS5tb3ZhYmxlbGltaXRzR0ZKLyUnYWNjZW50R0ZKLyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGWS1JI21uR0YkNiRRIjFGJ0ZGRjJGNUZGRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGRg== similar to the symmetric matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdXBHRiQ2JUYrLUYjNiUtSShtZmVuY2VkR0YkNiQtRiM2JS1JI21uR0YkNiRRIjFGJ0Y5Ri9GMkY5Ri9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictRjY2LVEifkYnRjlGO0Y+RkBGQkZERkZGSC9GS1EmMC4wZW1GJy9GTkZeb0Y5 construct the following matrices LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNiUtSSNtbkdGJDYkUSIyRicvRjZRJ25vcm1hbEYnRjJGNUZDRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSIsRidGQy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0ZOLyUqc3ltbWV0cmljR0ZOLyUobGFyZ2VvcEdGTi8lLm1vdmFibGVsaW1pdHNHRk4vJSdhY2NlbnRHRk4vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlRi4tRiM2JS1GOzYkLUYjNiUtRkA2JFEiM0YnRkNGMkY1RkNGMkY1RkVGSC1GLDYlRi4tRiM2JS1GOzYkLUYjNiUtRkA2JFEiNEYnRkNGMkY1RkNGMkY1RkVGSC1GSTYtUSMuLkYnRkNGTC9GUEZORlFGU0ZVRldGWS9GZm5RLDAuMjIyMjIyMmVtRicvRmluRmduLUZJNi1RIi5GJ0ZDRkxGZHBGUUZTRlVGV0ZZRmVuRmdwRkgtRiw2JUYuLUYjNiUtRjs2JC1GIzYnLUYvNiVRIm5GJ0YyRjUtRkk2LVEoJm1pbnVzO0YnRkNGTEZkcEZRRlNGVUZXRllGZXAvRmluRmZwLUZANiRRIjFGJ0ZDRjJGNUZDRjJGNUZFRkM= , where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNiUtRi82JVEia0YnRjJGNUYyRjUvRjZRJ25vcm1hbEYnRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSI9RidGQi8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGTS8lKXN0cmV0Y2h5R0ZNLyUqc3ltbWV0cmljR0ZNLyUobGFyZ2VvcEdGTS8lLm1vdmFibGVsaW1pdHNHRk0vJSdhY2NlbnRHRk0vJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0Zmbi1GOzYkLUYjNiQtSShtc3Vic3VwR0YkNigtRi82JVEiYUYnRjJGNS1GIzYlLUYvNiVRI2lqRidGMkY1RjJGNUY4RkQvJS9zdWJzY3JpcHRzaGlmdEdGRi9JK21zZW1hbnRpY3NHRiRRLFtub25lLG5vbmVdRidGQkZCRkI= for each LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYvLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjFGJ0Y5LUY2Ni1RIixGJ0Y5RjsvRj9GMUZARkJGREZGRkgvRktRJjAuMGVtRicvRk5RLDAuMzMzMzMzM2VtRictRlA2JFEiMkYnRjlGUy1GNjYtUSMuLkYnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZYLUY2Ni1RIi5GJ0Y5RjtGPkZARkJGREZGRkhGV0Zdb0ZTLUYsNiVRIm5GJ0YvRjItRjY2LVEoJm1pbnVzO0YnRjlGO0Y+RkBGQkZERkZGSEZbby9GTkZcb0ZPRjk=.
<Text-field style="Heading 1" layout="Heading 1">Householder(n,A0,display)</Text-field> Householder:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix local A,i,j,k,l,q,u,v,w,z,alpha,PROD,temp,RSQ: #option trace: A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(1RSQ))*v = 1/2r * v v:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: for j from k+2 to n do v[j]:=A[j,k]: end do: #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A): end proc:
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=
<Text-field style="Heading 1" layout="Heading 1">Householder2(n,A0,display)</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYvLUkjbWlHRiQ2JVElVGhpc0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUSJ+RidGMi8lJmZlbmNlR0YxLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRjEvJSpzeW1tZXRyaWNHRjEvJShsYXJnZW9wR0YxLyUubW92YWJsZWxpbWl0c0dGMS8lJ2FjY2VudEdGMS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRkktRiw2JVEodmVyc2lvbkYnRi9GMkY1LUYsNiVRJXdpbGxGJ0YvRjJGNS1GLDYlUSZwcmludEYnRi9GMkY1LUYsNiVRJHRoZUYnRi9GMkY1LUYsNiVRKW1hdHJpY2VzRidGL0YyRjVGMg==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYrLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiUEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIixGJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRkkvJSpzeW1tZXRyaWNHRkkvJShsYXJnZW9wR0ZJLyUubW92YWJsZWxpbWl0c0dGSS8lJ2FjY2VudEdGSS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVGLi1GIzYlLUY7NiRRIjJGJ0Y+RjJGNUZARkMtRiw2JUYuLUYjNiUtRjs2JFEiM0YnRj5GMkY1RkBGQy1GRDYtUSMuLkYnRj5GRy9GS0ZJRkxGTkZQRlJGVC9GV1EsMC4yMjIyMjIyZW1GJy9GWkZYLUZENi1RIi5GJ0Y+RkdGZ29GTEZORlBGUkZURlZGam9GPg== where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiUEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RJyZzZG90O0YnRj4vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkkvJSlzdHJldGNoeUdGSS8lKnN5bW1ldHJpY0dGSS8lKGxhcmdlb3BHRkkvJS5tb3ZhYmxlbGltaXRzR0ZJLyUnYWNjZW50R0ZJLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGWC1GLDYlRi4tRiM2JS1GOzYkUSIyRidGPkYyRjVGQEZDLUZENi1RIy4uRidGPkZHRkpGTEZORlBGUkZUL0ZXUSwwLjIyMjIyMjJlbUYnRlktRkQ2LVEiLkYnRj5GR0ZKRkxGTkZQRlJGVEZWRllGQy1GLDYlRi4tRiM2JS1GLzYlUSJuRidGMkY1RjJGNUZALUZENi1RIj1GJ0Y+RkdGSkZMRk5GUEZSRlQvRldRLDAuMjc3Nzc3OGVtRicvRlpGX3AtRiw2JS1GLzYlUSJRRidGMkY1RmZvRkBGPg== Householder2:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix #option trace; local A,i,j,k,l,P,q,u,v,w,z,alpha,PROD,temp,RSQ; A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(RSQ))*v = 1/2r * v v:=Array(1..n): w:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: w[k+1]:=v[k+1]/(2*(0.5*RSQ)^(1/2)); for j from k+2 to n do v[j]:=A[j,k]: w[j]:=v[j]/(2*(0.5*RSQ)^(1/2)): end do: #Compute P^(k)=I-2w*w P:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do if i=j then P[i,j]:=1-2*w[i]*w[j]: else P[i,j]:=0-2*w[i]*w[j]: end if: end do: end do: print("v",v); print("w",w); print("P^",k,P); #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A); end proc:
JSFH
<Text-field style="Heading 1" layout="Heading 1">Householder3(n,A0,display)</Text-field> This version will print the values of T and Q such that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiVEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdXBHRiQ2JS1GLDYlUSRRQVFGJ0YvRjItRiM2JS1GLDYlUSJ0RidGL0YyRi9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGOQ== Householder3:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix #option trace; local A,i,j,k,l,P,q,Q,u,v,w,z,alpha,PROD,temp,RSQ,SqrMatrixProduct; SqrMatrixProduct:=proc(n,A,B) #multiplies two square nxn matrices A and B local i,j,k,C: C:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do for k from 1 to n do C[i,k]:=C[i,k]+A[i,j]*B[j,k]; end do: end do: end do: return(C); end proc: Q:=Array(1..n,1..n): A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(RSQ))*v = 1/2r * v v:=Array(1..n): w:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: w[k+1]:=v[k+1]/(2*(0.5*RSQ)^(1/2)); for j from k+2 to n do v[j]:=A[j,k]: w[j]:=v[j]/(2*(0.5*RSQ)^(1/2)): end do: #Compute P^(k)=I-2w*w P:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do if i=j then P[i,j]:=1-2*w[i]*w[j]: else P[i,j]:=0-2*w[i]*w[j]: end if: end do:k end do: # print("v",v); # print("w",w); # print("P^",k,P); if k=1 then Q:=P: else Q:=SqrMatrixProduct(n,Q,P): end if: #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A,Q); end proc:
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= A:=Array([[4.0,1,-2,2],[1.0,2,0,1],[-2.0,0,3,-2],[2.0,1,-2,-1]]); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIzo9RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtSShtYWN0aW9uR0YkNiUtSShtZmVuY2VkR0YkNigtRiM2Ji1JJ210YWJsZUdGJDY4LUkkbXRyR0YkNiktSSRtdGRHRiQ2KC1JI21uR0YkNiRRJDQuMEYnRjkvJSlyb3dhbGlnbkdRIUYnLyUsY29sdW1uYWxpZ25HRmBvLyUrZ3JvdXBhbGlnbkdGYG8vJShyb3dzcGFuR1EiMUYnLyUrY29sdW1uc3BhbkdGZ28tRmhuNigtRltvNiRRIjFGJ0Y5Rl5vRmFvRmNvRmVvRmhvLUZobjYoLUYjNiUtRjY2LVEqJnVtaW51czA7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORmdwLUZbbzYkUSIyRidGOUY5Rl5vRmFvRmNvRmVvRmhvLUZobjYoRmlwRl5vRmFvRmNvRmVvRmhvRl5vRmFvRmNvLUZlbjYpLUZobjYoLUZbbzYkUSQxLjBGJ0Y5Rl5vRmFvRmNvRmVvRmhvRlxxLUZobjYoLUZbbzYkUSIwRidGOUZeb0Zhb0Zjb0Zlb0Zob0Zqb0Zeb0Zhb0Zjby1GZW42KS1GaG42KC1GIzYlRmNwLUZbbzYkUSQyLjBGJ0Y5RjlGXm9GYW9GY29GZW9GaG9GZXEtRmhuNigtRltvNiRRIjNGJ0Y5Rl5vRmFvRmNvRmVvRmhvRl9wRl5vRmFvRmNvLUZlbjYpLUZobjYoRmByRl5vRmFvRmNvRmVvRmhvRmpvRl9wLUZobjYoLUYjNiVGY3BGXHBGOUZeb0Zhb0Zjb0Zlb0Zob0Zeb0Zhb0Zjby8lJmFsaWduR1ElYXhpc0YnL0Zfb1EpYmFzZWxpbmVGJy9GYm9RJ2NlbnRlckYnL0Zkb1EnfGZybGVmdHxockYnLyUvYWxpZ25tZW50c2NvcGVHRjEvJSxjb2x1bW53aWR0aEdRJWF1dG9GJy8lJndpZHRoR0ZddC8lK3Jvd3NwYWNpbmdHUSYxLjBleEYnLyUuY29sdW1uc3BhY2luZ0dRJjAuOGVtRicvJSlyb3dsaW5lc0dRJW5vbmVGJy8lLGNvbHVtbmxpbmVzR0ZodC8lJmZyYW1lR0ZodC8lLWZyYW1lc3BhY2luZ0dRLDAuNGVtfjAuNWV4RicvJSplcXVhbHJvd3NHRj0vJS1lcXVhbGNvbHVtbnNHRj0vJS1kaXNwbGF5c3R5bGVHRj0vJSVzaWRlR1EmcmlnaHRGJy8lMG1pbmxhYmVsc3BhY2luZ0dGZXQvJStmb3JlZ3JvdW5kR1EoWzAsMCwwXUYnLyUpcmVhZG9ubHlHRj1GOUY5L0krbXNlbWFudGljc0dGJFEnTWF0cml4RicvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRidGYHYvJSthY3Rpb250eXBlR1EucnRhYmxlYWRkcmVzc0YnLyUpcnRhYmxlaWRHUSoyNTE0MDk2OTZGJ0ZbdkZedkY5 Householder(4,A,1); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbW5HRiQ2JFEiMUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RIixGJ0YvLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjgvJSpzeW1tZXRyaWNHRjgvJShsYXJnZW9wR0Y4LyUubW92YWJsZWxpbWl0c0dGOC8lJ2FjY2VudEdGOC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUkobWFjdGlvbkdGJDYlLUkobWZlbmNlZEdGJDYoLUYjNiYtSSdtdGFibGVHRiQ2OC1JJG10ckdGJDYpLUkkbXRkR0YkNigtRiw2JFEkNC4wRidGLy8lKXJvd2FsaWduR1EhRicvJSxjb2x1bW5hbGlnbkdGXG8vJStncm91cGFsaWduR0Zcby8lKHJvd3NwYW5HUSIxRicvJStjb2x1bW5zcGFuR0Zjby1GZW42KC1GIzYlLUYzNi1RKiZ1bWludXMwO0YnRi9GNi9GOkY4RjxGPkZARkJGRC9GR1EsMC4yMjIyMjIyZW1GJy9GSkZfcC1GLDYkUSwzLjAwMDAwMDAwMEYnRi9GL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GLDYkUSIwRidGL0ZqbkZdb0Zfb0Zhb0Zkb0ZkcEZqbkZdb0Zfby1GWDYpRmZvLUZlbjYoLUYsNiRRLDMuMzMzMzMzMzM0RidGL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GLDYkUS0wLjk5OTk5OTk5OTZGJ0YvRmpuRl1vRl9vRmFvRmRvLUZlbjYoLUYsNiRRLDEuMzMzMzMzMzMzRidGL0ZqbkZdb0Zfb0Zhb0Zkb0ZqbkZdb0Zfby1GWDYpRmRwRmBxLUZlbjYoLUYsNiRRLDEuNjY2NjY2NjY3RidGL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GIzYlRmpvRmdxRi9Gam5GXW9GX29GYW9GZG9Gam5GXW9GX28tRlg2KUZkcEZlcUZhci1GZW42KC1GIzYlRmpvLUYsNiRRIzEuRidGL0YvRmpuRl1vRl9vRmFvRmRvRmpuRl1vRl9vLyUmYWxpZ25HUSVheGlzRicvRltvUSliYXNlbGluZUYnL0Zeb1EnY2VudGVyRicvRmBvUSd8ZnJsZWZ0fGhyRicvJS9hbGlnbm1lbnRzY29wZUdGOy8lLGNvbHVtbndpZHRoR1ElYXV0b0YnLyUmd2lkdGhHRlt0LyUrcm93c3BhY2luZ0dRJjEuMGV4RicvJS5jb2x1bW5zcGFjaW5nR1EmMC44ZW1GJy8lKXJvd2xpbmVzR1Elbm9uZUYnLyUsY29sdW1ubGluZXNHRmZ0LyUmZnJhbWVHRmZ0LyUtZnJhbWVzcGFjaW5nR1EsMC40ZW1+MC41ZXhGJy8lKmVxdWFscm93c0dGOC8lLWVxdWFsY29sdW1uc0dGOC8lLWRpc3BsYXlzdHlsZUdGOC8lJXNpZGVHUSZyaWdodEYnLyUwbWlubGFiZWxzcGFjaW5nR0ZjdC8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJSlyZWFkb25seUdGOEYvRi8vSSttc2VtYW50aWNzR0YkUSdNYXRyaXhGJy8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Zedi8lK2FjdGlvbnR5cGVHUS5ydGFibGVhZGRyZXNzRicvJSlydGFibGVpZEdRKjI1MTQwOTY5NkYnRml1Rlx2Ri8= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbW5HRiQ2JFEiMkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RIixGJ0YvLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjgvJSpzeW1tZXRyaWNHRjgvJShsYXJnZW9wR0Y4LyUubW92YWJsZWxpbWl0c0dGOC8lJ2FjY2VudEdGOC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUkobWFjdGlvbkdGJDYlLUkobWZlbmNlZEdGJDYoLUYjNiYtSSdtdGFibGVHRiQ2OC1JJG10ckdGJDYpLUkkbXRkR0YkNigtRiw2JFEkNC4wRidGLy8lKXJvd2FsaWduR1EhRicvJSxjb2x1bW5hbGlnbkdGXG8vJStncm91cGFsaWduR0Zcby8lKHJvd3NwYW5HUSIxRicvJStjb2x1bW5zcGFuR0Zjby1GZW42KC1GIzYlLUYzNi1RKiZ1bWludXMwO0YnRi9GNi9GOkY4RjxGPkZARkJGRC9GR1EsMC4yMjIyMjIyZW1GJy9GSkZfcC1GLDYkUSwzLjAwMDAwMDAwMEYnRi9GL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GLDYkUSIwRidGL0ZqbkZdb0Zfb0Zhb0Zkb0ZkcEZqbkZdb0Zfby1GWDYpRmZvLUZlbjYoLUYsNiRRLDMuMzMzMzMzMzM0RidGL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GIzYlRmpvLUYsNiRRLDEuNjY2NjY2NjY2RidGL0YvRmpuRl1vRl9vRmFvRmRvRmRwRmpuRl1vRl9vLUZYNilGZHBGYHEtRmVuNigtRiM2JUZqby1GLDYkUSwxLjMyMDAwMDAwM0YnRi9GL0ZqbkZdb0Zfb0Zhb0Zkby1GZW42KC1GLDYkUS0wLjkwNjY2NjY2ODZGJ0YvRmpuRl1vRl9vRmFvRmRvRmpuRl1vRl9vLUZYNilGZHBGZHBGYHItRmVuNigtRiw2JFEsMS45ODY2NjY2NjhGJ0YvRmpuRl1vRl9vRmFvRmRvRmpuRl1vRl9vLyUmYWxpZ25HUSVheGlzRicvRltvUSliYXNlbGluZUYnL0Zeb1EnY2VudGVyRicvRmBvUSd8ZnJsZWZ0fGhyRicvJS9hbGlnbm1lbnRzY29wZUdGOy8lLGNvbHVtbndpZHRoR1ElYXV0b0YnLyUmd2lkdGhHRmlzLyUrcm93c3BhY2luZ0dRJjEuMGV4RicvJS5jb2x1bW5zcGFjaW5nR1EmMC44ZW1GJy8lKXJvd2xpbmVzR1Elbm9uZUYnLyUsY29sdW1ubGluZXNHRmR0LyUmZnJhbWVHRmR0LyUtZnJhbWVzcGFjaW5nR1EsMC40ZW1+MC41ZXhGJy8lKmVxdWFscm93c0dGOC8lLWVxdWFsY29sdW1uc0dGOC8lLWRpc3BsYXlzdHlsZUdGOC8lJXNpZGVHUSZyaWdodEYnLyUwbWlubGFiZWxzcGFjaW5nR0ZhdC8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJSlyZWFkb25seUdGOEYvRi8vSSttc2VtYW50aWNzR0YkUSdNYXRyaXhGJy8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Zcdi8lK2FjdGlvbnR5cGVHUS5ydGFibGVhZGRyZXNzRicvJSlydGFibGVpZEdRKjI1MTQwOTY5NkYnRmd1Rmp1Ri8= LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKidwNDlE LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQwOTY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIzEiJSIlJCIjUyEiIiQhKysrKytJISIqIiIhRi1GKiQiK01MTExMRiwkISttbW1tO0YsRi1GLUYwJCErLisrPzhGLCQiKydvbW0xKiEjNUYtRi1GNCQiK29tbScpPkYsNiI=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQwOTY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIzEiJSIlJCIjUyEiIiQhKysrKytJISIqIiIhRi1GKiQiK01MTExMRiwkISttbW1tO0YsRi1GLUYwJCErLisrPzhGLCQiKydvbW0xKiEjNUYtRi1GNCQiK29tbScpPkYsNiI=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQwOTY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIzEiJSIlJCIjUyEiIiQhKysrKytJISIqIiIhRi1GKiQiK01MTExMRiwkISttbW1tO0YsRi1GLUYwJCErLisrPzhGLCQiKydvbW0xKiEjNUYtRi1GNCQiK29tbScpPkYsNiI=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQwOTY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIzEiJSIlJCIjUyEiIiQhKysrKytJISIqIiIhRi1GKiQiK01MTExMRiwkISttbW1tO0YsRi1GLUYwJCErLisrPzhGLCQiKydvbW0xKiEjNUYtRi1GNCQiK29tbScpPkYsNiI=