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-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

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