jacobi.sa

 
---------------------------> Sather 1.1 source file <--------------------------
-- jacobi.sa: jacobi method from NR
-- Author: Alex Cozzi <cozzi@neurop2.ruhr-uni-bochum.de>
 
-- COPYRIGHT NOTICE: This code is provided WITHOUT ANY WARRANTY
-- and is subject to the terms of the SATHER LIBRARY GENERAL PUBLIC
-- LICENSE contained in the file: Sather/Doc/License of the
-- Sather distribution. The license is also available from ICSI,
-- 1947 Center St., Suite 600, Berkeley CA 94704, USA.


class JACOBI
class JACOBI
is
   jacobi(aa :MAT, out d :VEC, out v:MAT) :INT
   -- Compute all eigenvalues and eigenvectors of a real symmetric matrix 
   -- a[N][N]. On output elements of a above the disgonal are destroyed.   
   -- d[N] returns the eigenvalues of a. v[N][N] is a matrix whose columns 
   -- contain, on output, the normalized eigenvectors of a. nrot returns   
   -- the number of Jacobi rotations that were required.		       
      pre aa.size1 = aa.size2 -- and aa.is_simmetrical
   is
      thresh, theta, tau, t, sm, s, h, g, c :FLT;
      n ::= aa.size1;
      a ::= aa.copy;
      d := #(n);
      v := #(n, n);
      loop v.set!(0.0); end;

      loop ip ::= v.ind1!;
	 v[ip, ip] := 1.0;
      end;
       
      --Initialize b and d to the diagonals of a
      b ::= #VEC(n);
      z ::= #VEC(n);

      loop ip ::= v.ind1!;
	 b[ip] := a[ip, ip];
	 d[ip] := a[ip, ip];
	 z[ip] := 0.0;
      end;

      nrot ::= 0;
      
      loop i ::= 0.upto!(99);
	 sm := 0.0;
	 loop ip ::= 0.upto!(n-2);
	    loop iq ::= (ip+1).upto!(n-1);
	       sm := sm + a[ip, iq].abs;
	    end;
	 end;

	 -- Normal return, which relies on quadratic convergence to
	 -- machine underflow 
	 if sm = 0.0 then  return nrot; end;

	 if i<3 
	 then thresh := 0.2 * sm/(n*n).flt; -- on the first three sweeps
	 else thresh := 0.0; --  the rest of the sweeps 
	 end;

	 loop ip ::= 0.upto!(n-2);
	    loop iq ::= (ip+1).upto!(n-1);
	       g := 100.0 * a[ip, iq].abs;
	       -- After 4 sweeps skip the rotation if the
	       -- off diagonal element is small
	       if i>3 and d[ip].abs+g = d[ip].abs
			and d[iq].abs+g = d[iq].abs
	       then a[ip, iq] := 0.0;
	       elsif a[ip, iq].abs > thresh 
	       then
		  h := d[iq]-d[ip];
		  if h.abs+g = h.abs 
		  then t := a[ip, iq]/h;
		  else 
		     theta := 0.5*h/a[ip, iq];
		     t := 1.0/(theta.abs + (1.0 + theta.square).sqrt);
		     if theta < 0.0 then t := -t; end;
		  end;
		  c := 1.0/(1.0+t*t).sqrt;
		  s := t*c;
		  tau := s/(1.0+c);
		  h := t*a[ip, iq];
		  z[ip] := z[ip] - h;
		  z[iq] := z[iq] + h;
		  d[ip] := d[ip] - h;
		  d[iq] := d[iq] + h;
		  a[ip, iq] := 0.0;
		  loop j::= 0.upto!(ip-2);
		     rotate(a, j, ip, j, iq, inout h, inout g, s, tau);
		  end;
		  loop j::= (ip+1).upto!(iq-2);
		     rotate(a, ip, j, j, iq, inout h, inout g, s, tau);
		  end;
		  loop j::= (iq+1).upto!(n-1);
		     rotate(a, ip, j, iq, j, inout h, inout g, s, tau);
		  end;
		  loop j::= 0.upto!(n-1);
		     rotate(v, j, ip, j, iq, inout h, inout g, s, tau);
		  end;
	       nrot := nrot + 1;
	       end;
	    end;
	 end;
	 loop ip::=0.upto!(n-1);
	    b[ip] := b[ip] + z[ip];
	    d[ip] := b[ip];
	    z[ip] := 0.0;
	 end;
      end;
      raise "\nFatal error: too many iterations in jacobi\n";
   end;

   private rotate(a :MAT, i, j, k, l :INT, inout h, inout g, s, tau :FLT)
   -- Do rotations required by Jacobi Transformation		    
   is
      g := a[i, j];
      h := a[k, l];
      a[i, j] := g - s*(h+g*tau);
      a[k, l] := h + s*(g-h*tau);
   end;
end;


class TEST_JACOBI
class TEST_JACOBI
is
   include TEST;
   main is
      class_name("JACOBI");
      m ::= #MAT(||1.0, 0.0| , |0.0, 1.0||);
      v :VEC;
      o :MAT;
      n ::= JACOBI::jacobi(m, out v, out o);
      test ("eigenvalues 1", v.str, "|1,1|");
      test ("eigenvectors 1", o.str, "||1,0|,|0,1||");
      m2::= #MAT(| |1.0, 2.0, 0.0|, |2.0, 1.0, 0.0|, |0.0, 0.0, -1.0| |);
      n := JACOBI::jacobi(m2, out v, out o);
      test ("eigenvalues 2", v.str, "|-1,3,-1|");
      test ("eigenvectors 2", o.str, 
	    "||0.707107,0.707107,0|,|-0.707107,0.707107,0|,|0,0,1||");
      m3 ::= #MAT(|| 4.0, -1.0, -1.0, -1.0 |,| -1.0, 4.0, -1.0, -1.0 |,
    		   | -1.0, -1.0, 4.0, -1.0 |,| -1.0, -1.0, -1.0, 4.0 ||);
      n := JACOBI::jacobi(m3, out v, out o);
      test ("eigenvalues 3", v.str, "|5,1,5,5|");
      finish;
   end;
end;