TABLE OF CONTENTS

ABINIT/jacobi [ Functions ]

[ Top ] [ Functions ]

NAME

  jacobi

FUNCTION

     Computes all eigenvalues and eigenvectors of a real symmetric matrix a,
     which is of size n by n, stored in a physical np by np array. On output,
     elements of a above the diagonal are destroyed. d returns the
     eigenvalues of a in its first n elements. v is a matrix with the same
     logical and physical dimensions as a, whose columns contain, on output,
     the normalized eigenvectors of a. nrot returns the number of Jacobi
     rotations that were required.

COPYRIGHT

  Copyright (C) 2014-2017 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      conducti_nc,critic

CHILDREN

SOURCE

 36 #if defined HAVE_CONFIG_H
 37 #include "config.h"
 38 #endif
 39 
 40 #include "abi_common.h"
 41 
 42 subroutine jacobi(a,n,np,d,v,nrot)
 43 
 44  use defs_basis, only : std_out
 45 
 46 !This section has been created automatically by the script Abilint (TD).
 47 !Do not modify the following lines by hand.
 48 #undef ABI_FUNC
 49 #define ABI_FUNC 'jacobi'
 50 !End of the abilint section
 51 
 52  implicit none
 53 !Arguments
 54  integer :: n,np,nrot
 55  real*8 :: a(np,np),d(np),v(np,np)
 56 !Local variables
 57  integer, parameter :: NMAX=500
 58  integer i,ip,iq,j
 59  real*8 c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
 60  do ip=1,n 
 61    do iq=1,n
 62      v(ip,iq)=0.
 63    enddo
 64    v(ip,ip)=1.
 65  enddo
 66  do ip=1,n
 67    b(ip)=a(ip,ip) 
 68    d(ip)=b(ip)
 69    z(ip)=0. 
 70  enddo
 71  nrot=0
 72  do i=1,50
 73    sm=0.
 74    do ip=1,n-1 
 75      do iq=ip+1,n
 76        sm=sm+abs(a(ip,iq))
 77      enddo
 78    enddo
 79    if(sm.eq.0.)return 
 80    if(i.lt.4)then
 81      tresh=0.2*sm/n**2 
 82    else
 83      tresh=0. 
 84    endif
 85    do ip=1,n-1
 86      do iq=ip+1,n
 87        g=100.*abs(a(ip,iq))
 88        if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) &
 89 &          .and.(abs(d(iq))+g.eq.abs(d(iq))))then
 90             a(ip,iq)=0.
 91        else if(abs(a(ip,iq)).gt.tresh)then
 92          h=d(iq)-d(ip)
 93          if(abs(h)+g.eq.abs(h))then
 94            t=a(ip,iq)/h 
 95          else
 96            theta=0.5*h/a(ip,iq) 
 97            t=1./(abs(theta)+sqrt(1.+theta**2))
 98            if(theta.lt.0.)t=-t
 99          endif
100          c=1./sqrt(1+t**2)
101          s=t*c
102          tau=s/(1.+c)
103          h=t*a(ip,iq)
104          z(ip)=z(ip)-h
105          z(iq)=z(iq)+h
106          d(ip)=d(ip)-h
107          d(iq)=d(iq)+h
108          a(ip,iq)=0.
109          do j=1,ip-1
110            g=a(j,ip)
111            h=a(j,iq)
112            a(j,ip)=g-s*(h+g*tau)
113            a(j,iq)=h+s*(g-h*tau)
114          enddo
115          do j=ip+1,iq-1
116            g=a(ip,j)
117            h=a(j,iq)
118            a(ip,j)=g-s*(h+g*tau)
119            a(j,iq)=h+s*(g-h*tau)
120          enddo
121          do j=iq+1,n 
122            g=a(ip,j)
123            h=a(iq,j)
124            a(ip,j)=g-s*(h+g*tau)
125            a(iq,j)=h+s*(g-h*tau)
126          enddo
127          do j=1,n
128            g=v(j,ip)
129            h=v(j,iq)
130            v(j,ip)=g-s*(h+g*tau)
131            v(j,iq)=h+s*(g-h*tau)
132          enddo
133          nrot=nrot+1
134        endif
135      enddo
136    enddo
137    do ip=1,n
138      b(ip)=b(ip)+z(ip)
139      d(ip)=b(ip) 
140      z(ip)=0. 
141    enddo
142  enddo
143  write(std_out,*) 'too many iterations in jacobi'
144 
145 end subroutine jacobi