'**********************************************************************
'*        Eigenvalues and eigenvectors of a real symmetric            *
'*                square matrix by JACOBI'S METHOD                    *
'* ------------------------------------------------------------------ *
'* Ref.: "NUMERICAL RECIPES, Cambridge University Press, 1986,        *
'*        chap. 11, pages 346-348" [BIBLI 08].                        *
'*                                                                    *
'*                                     Basic version by J-P Moreau    *
'*                                          (www.jpmoreau.fr)         *
'* ------------------------------------------------------------------ *
'* SAMPLE RUN:                                                        *
'* Input file (elpro.dat):                                            *
'*                                                                    *
'*  4                                                                 *
'*   4 -2 -1  0                                                       *
'*      4  0 -1                                                       *
'*         4 -2                                                       *
'*            4                                                       *
'*                                                                    *
'* Output file (elpro.lst):                                           *
'*                                                                    *
'* ---------------------------------------------------                *
'*  Eigenvalues and eigenvectors of a real, symmetric                 *
'*     square matrix by Jacobi's iterative method                     *
'* ---------------------------------------------------                *
'*                                                                    *
'*  Dimension of matrix: 4                                            *
'*                                                                    *
'*  Symmetric matrix is:                                              *
'*                                                                    *
'*     4.000000    -2.000000    -1.000000     0.000000                *
'*    -2.000000     4.000000     0.000000    -1.000000                *
'*    -1.000000     0.000000     4.000000    -2.000000                *
'*     0.000000    -1.000000    -2.000000     4.000000                *
'*                                                                    *
'*  Eigenvalues:                                                      *
'*     7.00000000                                                     *
'*     1.00000000                                                     *
'*     3.00000000                                                     *
'*     5.00000000                                                     *
'*                                                                    *
'*  Eigenvectors (in lines):                                          *
'*     1.000000    -1.000000    -1.000000     1.000000                *
'*     1.000000     1.000000     1.000000     1.000000                *
'*     1.000000     1.000000    -1.000000    -1.000000                *
'*    -1.000000     1.000000    -1.000000     1.000000                *
'*                                                                    *
'*  Number of rotations: 25                                           *
'*                                                                    *
'**********************************************************************
'PROGRAM Test_Jacobi
defdbl a-h,o-z
defint i-n

' integer idim                           Dimension of square matrix
'                                        (called N in subroutine 1000)
' integer nrot                           Number of matrix rotations
' double  vmax                           Normalization factor

dim A(20,20)                            'given square symmetric matrix
dim D(20)                               'eigenvalues solution vector
dim V(20,20)                            'eigenvectors solution matrix

' open input/output files
  open "elpro.dat" for input as #1       'read mode
  open "elpro.lst" for output as #2      'write mode

' read size of matrix from input file
  input #1, idim

' read data from input file
  for i=1 to idim                       'read only upper half triangle
    for j=i to idim
      input #1, A(i,j)
      A(j,i)=A(i,j)
    next j
  next i
  close #1

' print title and input matrix
  cls
  print #2,"---------------------------------------------------"
  print #2," Eigenvalues and eigenvectors of a real, symmetric "
  print #2,"    square matrix by Jacobi's iterative method    "
  print #2,"---------------------------------------------------"
  print #2,""
  print #2," Dimension of matrix: ";  idim
  print #2,""
  print #2," Symmetric matrix is:"
  for i=1 to idim
    for j=1 to idim
      print #2, using "####.####"; A(i,j);
    next j
    print #2,""
  next i

' call jacobi subroutine
  gosub 1000

' print results to output file
  print #2,""
  print #2," Eigenvalues:"
  for i=1 to idim
    print #2, using "####.######"; D(i)
  next i
  print #2,""
  print #2," Eigenvectors (in lines):"
  for j=1 to idim
    vmax=V(1,j)
    for i=2 to idim
      if(abs(V(i,j))>abs(vmax)) then vmax=V(i,j)
    next i
    for i=1 to idim
      print #2, using " ####.######"; V(i,j)/vmax;
    next i
    print #2,""
  next j
  print #2,""
  print #2," Number of rotations: "; nrot
  print #2,"---------------------------------------------------"
  close #2
  print
  print " Program done."
  print " Results in file elpro.lst."
  print

END

'*************************************************************
'* This subroutine computes all eigenvalues and eigenvectors *
'* of a real symmetric square matrix A(N,N). On output, ele- *
'* ments of A above the diagonal are destroyed. D(N) returns *
'* the eigenvalues of matrix A. V(N,N) contains, on output,  *
'* the eigenvectors of A by columns. THe normalization to    *
'* unity is made by main program before printing results.    *
'* NROT returns the number of Jacobi matrix rotations which  *
'* were required.                                            *
'*************************************************************
1000 'Subroutine Jacobi1(A,N,D,V,NROT)
dim B(20), Z(20)
  N=idim
  for ip=1 to N    'initialize V to identity matrix
    for iq=1 to N
      V(ip,iq)=0#
    next iq
    V(ip,ip)=1#
  next ip
  for ip=1 to N
    B(ip)=A(ip,ip)
    D(ip)=B(ip)
    Z(ip)=0#
  next ip
  nrot=0
  for i=1 to 50
    sm=0#
    for ip=1 to N-1     'sum off-diagonal elements
      for iq=ip+1 to N
        sm=sm+ABS(A(ip,iq))
      next iq
    next ip
    if sm=0# then return  'normal return
    if i < 4 then
      tresh=0.2#*sm*sm
    else
      tresh=0#
    end if
    for ip=1 to N-1
      for iq=ip+1 to N
        g=100#*ABS(A(ip,iq))
' after 4 sweeps, skip the rotation if the off-diagonal element is small
        if i>4 and ABS(D(ip))+g=ABS(D(ip)) and ABS(D(iq))+g=ABS(D(iq)) then
	  A(ip,iq)=0#
        elseif ABS(A(ip,iq)) > tresh then
          h=D(iq)-D(ip)
          if ABS(h)+g=ABS(h) then
	    t=A(ip,iq)/h
          else
	    theta=0.5#*h/A(ip,iq)
            t=1#/(ABS(theta)+SQR(1#+theta*theta))
	    if theta < 0# then t=-t
          end if
	  c=1#/SQR(1#+t*t)
	  s=t*c
          tau=s/(1#+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.d0
	  for j=1 to ip-1
	    g=A(j,ip)
	    h=A(j,iq)
	    A(j,ip)=g-s*(h+g*tau)
	    A(j,iq)=h+s*(g-h*tau)
          next j
	  for j=ip+1 to iq-1
	    g=A(ip,j)
	    h=A(j,iq)
	    A(ip,j)=g-s*(h+g*tau)
	    A(j,iq)=h+s*(g-h*tau)
          next j		
	  for j=iq+1 to N
	    g=A(ip,j)
	    h=A(iq,j)
	    A(ip,j)=g-s*(h+g*tau)
	    A(iq,j)=h+s*(g-h*tau)
          next j		
	  for j=1 to N
	    g=V(j,ip)
	    h=V(j,iq)
	    V(j,ip)=g-s*(h+g*tau)
	    V(j,iq)=h+s*(g-h*tau)
          next j		
          nrot=nrot+1
        end if  'if ((i.gt.4)...
      next iq   'main iq loop
    next ip     'main ip loop
    for ip=1 to N
      B(ip)=B(ip)+Z(ip)
      D(ip)=B(ip)
      Z(ip)=0#
    next ip
  next i        'main i loop
  print " 50 iterations !"
return

' end of file tujacobi.bas