EigenVectors

[DoumbiaS]

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

[DoumbiaS]

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 · Accepted solution]

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]

Thank you ! that's great !

Best Regards