Hello,
It exists a JK METHOD for finding eigenvalues and eigenvector of a real symmetric matrix.
Is this method implemented in one sas procs ?
Regards
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!
Hello,
It exists a JK METHOD for finding eigenvalues and eigenvector of a real symmetric matrix.
Is this method implemented in one sas procs ?
Regards
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!
Thank you ! that's great !
Best Regards
SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.