BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
DoumbiaS
Quartz | Level 8

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

1 ACCEPTED SOLUTION

Accepted Solutions
JacobSimonsen
Barite | Level 11

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 solution in original post

3 REPLIES 3
DoumbiaS
Quartz | Level 8

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

JacobSimonsen
Barite | Level 11

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! 

 

DoumbiaS
Quartz | Level 8

Thank you ! that's great !

 

Best Regards 

sas-innovate-2026-white.png



April 27 – 30 | Gaylord Texan | Grapevine, Texas

Registration is open

Walk in ready to learn. Walk out ready to deliver. This is the data and AI conference you can't afford to miss.
Register now and lock in 2025 pricing—just $495!

Register now

What is ANOVA?

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.

Discussion stats
  • 3 replies
  • 4056 views
  • 2 likes
  • 2 in conversation