you can use PROC IML like this:
PROC IML;
A = {1 2, 2 1};
call eigen(val, rvec, A) vecl="lvec";
print val;
print rvec;
QUIT;
Or you can define a function in fcmp which do the work. It becomes quite complicated, but work without using IML.
option cmplib=work.func;
proc fcmp outlib=work.func.matrix;
subroutine copy(a[*,*],b[*,*]);
outargs b;
do i=1 to dim(a,1);
do j=1 to dim(a,2);
b[i,j]=a[i,j];
end;
end;
endsub;
subroutine rotate(m[*,*],i,j,k,l,s,tau);
outargs m,i,j,k,l,s,tau;
*this routine is used by the jacobi-macro;
g=m[i,j];
h=m[k,l];
m[i,j]=g-s*(h+g*tau);
m[k,l]=h+s*(g-h*tau);
endsub;
subroutine jacobi(a[*,*],d[*],v[*,*],nrot);
array b{1} _temporary_;
array z{1} _temporary_;
array copy{1,1} _temporary_;
call dynamic_array(copy, dim(a,1),dim(a,1));
call dynamic_array(z, dim(a,1));
call dynamic_array(b, dim(a,1));
outargs d,v,nrot;
n=dim(a,1);
call identity(v);
call copy(a,copy);
do ip=1 to n;
b[ip]=copy[ip,ip];
d[ip]=copy[ip,ip];
z[ip]=0;
end;
nrot=0;
do i=1 to 50;
sm=0;
do ip=1 to (n-1);
do iq=(ip+1) to n;
sm=sm+abs(copy[ip,iq]);
end;
end;
if sm<0.000000001 then do;
i=51;
end;
else if i<4 then do;
tresh=0.2*sm/(n**2);
end;
else tresh=0;
do ip=1 to (n-1);
do iq=(ip+1) to n;
g=100*abs(copy[ip,iq]);
if ((i>4)*(abs(d[ip])+g=abs(d[ip]))*(abs(d[iq])+g=abs(d[iq]))) then copy[ip,iq]=0;
else if (abs(copy[ip,iq]) > tresh) then do;
h=d[iq]-d[ip];
if (abs(h)+g) =abs(h) then t=copy[ip,iq]/h;
else do;
theta=0.5*h/copy[ip,iq];
t=1/(abs(theta)+sqrt(1+theta**2));
if (theta<0) then t=-t;
end;
_c_=1/sqrt(1+t**2);
s=t*_c_;
tau=s/(1.0+_c_);
h=t*copy[ip,iq];
z[ip] =z[ip]- h;
z[iq] =z[iq]+ h;
d[ip] = d[ip]-h;
d[iq] = d[iq]+h;
copy[ip,iq]=0;
do j=1 to (ip-1);
call rotate(copy,j,ip,j,iq,s,tau);
end;
do j=(ip+1) to (iq-1);
call rotate(copy,ip,j,j,iq,s,tau);
end;
do j=(iq+1) to n;
call rotate(copy,ip,j,iq,j,s,tau);
end;
do j=1 to n;
call rotate(v,j,ip,j,iq,s,tau);
end;
nrot=nrot+1;
end;
end;
end;
do ip=1 to n;
b[ip] = b[ip]+ z[ip];
d[ip]=b[ip];
z[ip]=0;
end;
end;
endsub;
subroutine show(m[*,*]);
do i=1 to dim(m,1);
do j=1 to dim(m,2);
put m[i,j] @@;
end;
put;
end;
endsub;
quit;
data _NULL_;
array A{2,2} _temporary_ (1,2,2,1);
array vectors{2,2} _temporary_;
array values{2} _temporary_;
call show(A);
*eigenvectors;
call jacobi(A,values,vectors,nrot);
call show(vectors);
run;
good luck!
... View more
