'********************************************************************** '* 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