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 2025: Save the Date

 SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!

Save the date!

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
  • 3196 views
  • 2 likes
  • 2 in conversation