Subject: Re: diagonalize a matrix

From: [email-masked] (Teerakiat Kerdcharoen)
Date: Tue, 11 Jan 1994 18:07:00 +0100 (NFT)

: 
:      Dear Netters
: 
: Has somebody a program, which allows to diagonalize a 3*3
: matrix as easily as possible
: 
:    thank you
: 
: Kempf Jean-Yves

Dear Dr. Kempf,
	I have enclosed a routine for matrix diagonalisation using
JACOBI method ( The code is free, but copyrighted by CCP5 group,
Daresbury, UK ). I would like to add an important note about the
algorithm. The Jacobi method, as one knows, modifies only the lower
triangle of the matrix. Once you use the code, you have to map the
lower triangle which is already modified to the upper which is left
unchanged, for example

	In the main program:

        call jacobi(mat1,mat2,3,1.0E-8)
        mat2(1,2) = mat2(2,1)
        mat2(1,3) = mat2(3,1)
        mat2(2,3) = mat2(3,2)

	Apart from the code, you may obtain other matrix numerical
codes from  [email-masked] by anonymously ftping to netlib2.cs.utk.edu
or by e-mail. Mail to [email-masked] having the content as follows,

	send index from eispack

	When you get the index of the code, if you want, i.e., diagonal
you can send command typically,

	send diagonal from eispack.
	
	( Note that the package name "diagonal" may not exist. it's just
          my example ).

	I believe you will enjoy other diagonalisation scheme beside
"jacobi".
							with best wishes,
							Teerakiat
----------------------------------------------------------------------------
Teerakiat Kerdcharoen (NOY)

Institute of General and Inorganic and Theoretical Chemistry
Innrain 52a, A-6020 Innsbruck AUSTRIA
e-mail:  [email-masked], [email-masked], [email-masked] 
      :  [email-masked], [email-masked] ( Bangkok )
Tel   :  (43)(512) 5073232 (office) , (43)(512) 2215201 (home)
FAX   :  (43)(512) 5073225
Research :  Molecular Dynamics simulations
         :  Computer Aided Molecular/Material Designs
-----------------------------------------------------------------------------
	
************************************************************************
*   U.3 JACOBI  Diagonalisation of symmetric matrices by jacobi method *
************************************************************************
      subroutine jacobi(a,v,n,eps)
      dimension a(n,n),v(n,n)
      logical any
C ----------------------------  Note  ---------------------------------
C ***      a = starting matrix
C ***      b = diagonalized matrix
C ***      n = dimension of matrix
C ***    eps = criteria for convergence
C ***          probably you can try  1.0E-6, 1.0E-8
C ***          The less the value is, the more accurate the outcome and
C ***          the more time spent.
C ***                                               Noted by:  Teerakiat
C ***
C ***      This part of the program is cited to CCP5, Daresbury
C ***      Collobarative Group, UK.
C ----------------------------------------------------------------------
*
* set v equal to identity
*
      do 10 j=1,n
      do 10 i=1,n
         v(i,j)=0.0
         if(i.eq.j) v(i,j)=1.0
   10 continue
*
* scale matrix a
*
      scl=0.0
      do 1 i=1,n
         scl=scl+a(i,i)**2
    1 continue
      scl=sqrt(scl)/3.0
      do 2 j=1,n
      do 2 i=1,n
         a(i,j)=a(i,j)/scl
    2 continue
*
* calculate quantity tes
*
      tes=0.0
      do 3 i=2,n
      do 3 j=1,i-1
         tes=tes+2.0*a(i,j)**2
    3 continue
      tes=sqrt(tes)
*
* loop while tes .gt. eps
*
      icnt=0
  100 continue
      icnt=icnt+1
*
* loop while array elements are changing
*
      jcnt=0
  110 continue
      jcnt=jcnt+1
      any=.false.
*
* loop over array elements
*
      do 120 i=2,n
      do 120 j=1,i-1
      if(abs(a(i,j)).le.tes) go to 120
      any=.true.
*
      v1=a(j,j)
      v2=a(i,j)
      v3=a(i,i)
      u=0.5*(v1-v3)
      if(abs(u).lt.eps) then
         omg=-1.0
      else
         omg=-v2/sqrt(v2**2+u**2)
         if(u.lt.0.0) omg=-omg
      endif
      s=omg/sqrt(2.0*(1.0+sqrt(1.0-omg**2)))
      c=sqrt(1.0-s**2)
*
      do 130 k=1,n
      if(k.ge.i) then
         temp=a(k,j)*c-a(k,i)*s
         a(k,i)=a(k,j)*s+a(k,i)*c
         a(k,j)=temp
      else
         if(k.lt.j) then
            temp=a(j,k)*c-a(i,k)*s
            a(i,k)=a(j,k)*s+a(i,k)*c
            a(j,k)=temp
         else
            temp=a(k,j)*c-a(i,k)*s
            a(i,k)=a(k,j)*s+a(i,k)*c
            a(k,j)=temp
         endif
      endif
*
      temp=v(k,j)*c-v(k,i)*s
      v(k,i)=v(k,j)*s+v(k,i)*c
      v(k,j)=temp
  130 continue
*
      a(j,j)=v1*c**2+v3*s**2-2.0*v2*s*c
      a(i,i)=v1*s**2+v3*c**2+2.0*v2*s*c
      a(i,j)=(v1-v3)*s*c+v2*(c**2-s**2)
  120 continue
*
      if(any) then
         if(jcnt.lt.50) then
            go to 110
         else
            print 101
         endif
      endif
*
      tes=tes/float(n)
      if(tes.gt.eps) then
         if(icnt.lt.50) then
            go to 100
         else
            print 101
         endif
      endif
*
* rescale matrix a
*
      do 4 j=1,n
      do 4 i=1,n
    4 a(i,j)=scl*a(i,j)
      return
  101 format(4x,'WARNING   Convergence problem in jacobi')
      end