jacobi.m

program for subroutine of Jacobi is given in this ex.

%=======================================================
%
%Subroutine Jacobi.m
%
%=======================================================
function [d,e]=jacobi(d,n,eps,itmax);
iter=0;
for i=1:n
   for j=1:n
      e(i,j)=0;
      e(i,i)=1.0;	
   end
end
while iter<itmax
zz=0.0;
nm1=n-1;
for i=1:nm1
   ip1=i+1;
   for j=ip1:n
      if abs(d(i,j))>zz
         zz=abs(d(i,j));
         ir=i;
         ic=j;
      end
   end
end
if iter==0
   yy=zz*eps;
end
if zz<yy
   return
end
dif=d(ir,ir)-d(ic,ic);
tanz=(-dif+sqrt(dif^2+4*zz^2))/(2.0*d(ir,ic));
cosz=1.0/sqrt(1.0+tanz^2);
sinz=cosz*tanz;
for i=1:n
   zzz=e(i,ir);
   e(i,ir)=cosz*zzz+sinz*e(i,ic);
   e(i,ic)=cosz*e(i,ic)-sinz*zzz;
end
i=1;
while i<ir
   yyy=d(i,ir);
   d(i,ir)=cosz*yyy+sinz*d(i,ic);
   d(i,ic)=cosz*d(i,ic)-sinz*yyy;
   i=i+1;
end
i=ir+1;
while i<ic
   yyy=d(ir,i);
   d(ir,i)=cosz*yyy+sinz*d(i,ic);
   d(i,ic)=cosz*d(i,ic)-sinz*yyy;
   i=i+1;
end
i=ic+1;
while i<=n
   zzz=d(ir,i);
   d(ir,i)=cosz*zzz+sinz*d(ic,i);
   d(ic,i)=cosz*d(ic,i)-sinz*zzz;
   i=i+1;
end
yyy=d(ir,ir);
d(ir,ir)=yyy*cosz^2+d(ir,ic)*2.0*cosz*sinz+d(ic,ic)*sinz^2;
d(ic,ic)=d(ic,ic)*cosz^2+yyy*sinz^2-d(ir,ic)*2*cosz*sinz;
d(ir,ic)=0.0;
iter=iter+1;
end