jacobview_v.F File Reference

#include "implicit_f.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	jacobiew_v (dim1, dim2, a, ew, ev, nrot)

Function/Subroutine Documentation

jacobiew_v()

subroutine jacobiew_v	(	integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		intent(inout)	a,
		intent(inout)	ew,
		intent(inout)	ev,
		integer, dimension(dim1)	nrot )

Definition at line 29 of file jacobview_v.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
        INTEGER, INTENT(in) :: DIM1
        INTEGER, INTENT(in) :: DIM2
        INTEGER, DIMENSION(DIM1) :: NROT
        my_real, DIMENSION(DIM1,DIM2), intent(inout) :: ew 
        my_real, DIMENSION(DIM1,DIM2,DIM2), intent(inout) :: a,ev
 
        my_real, DIMENSION(:,:), ALLOCATABLE :: b,z
        INTEGER IZ,IS,ITER,J,IJK,I
        INTEGER :: NINDX
        INTEGER, DIMENSION(DIM1) :: INDX
        my_real, DIMENSION(DIM1) :: sumrs,eps
        my_real :: g,h,t,c,s,tau,theta      
C----------------------------------------------------------------
 
        ! ---------------------
        ALLOCATE( b(dim1,dim2) )
        ALLOCATE( z(dim1,dim2) )
        DO iz=1,dim2
            DO is=1,dim2
                IF(iz>is) THEN 
                    DO i=1,dim1
                        a(i,is,iz) = a(i,iz,is) 
                    ENDDO
                ENDIF
                ev(1:dim1,iz,is)=zero
            ENDDO
            b(1:dim1,iz)=a(1:dim1,iz,iz)
            ew(1:dim1,iz)=b(1:dim1,iz)
            z(1:dim1,iz)=0.
            ev(1:dim1,iz,iz)=one
        ENDDO
        nrot(1:dim1)=0
        ! ---------------------
 
 
        ! ---------------------
        DO iter = 1,50
            sumrs(1:dim1) = zero
            ! ---------------------  
            ! diagonal sum          
            DO iz=1,dim2-1
                DO is=iz+1,dim2
                    sumrs(1:dim1)=sumrs(1:dim1)+abs(a(1:dim1,iz,is))
                ENDDO
            ENDDO
            ! ---------------------
            ! check if the computation is mandatory : sumrs(i) = 0 --> skip the element
            nindx = 0
            DO i=1,dim1
                IF (sumrs(i)/=zero) THEN
                    nindx = nindx + 1
                    indx(nindx) = i
                ENDIF
            ENDDO
            ! ---------------------
            IF (iter > 4)   THEN
                eps(1:dim1) = zero
            ELSE
                eps(1:dim1) = one_fifth*sumrs(1:dim1)/dim2**2
            ENDIF
            ! ---------------------
            DO iz=1,dim2-1
                DO is=iz+1,dim2
#include "vectorize.inc"                       
                    DO ijk=1,nindx
                        i = indx(ijk)
                        g = 100. * abs(a(i,iz,is))
                        IF ( iter>4 .AND. abs(ew(i,iz))+g==abs(ew(i,iz))
     &                          .AND. abs(ew(i,is))+g==abs(ew(i,is))) THEN
                            a(i,iz,is)=zero
                        ELSEIF (abs(a(i,iz,is)) > eps(i)) THEN
                            h = ew(i,is)-ew(i,iz)
                            IF (abs(h)+g==abs(h)) THEN
                                t = a(i,iz,is)/h
                            ELSE
                                theta = half*h/a(i,iz,is)
                                t=one/(abs(theta)+sqrt(one+theta**2))
                                IF (theta < zero) t=-t
                            ENDIF
                            c=one/sqrt(one+t**2)
                            s=t*c
                            tau=s/(one+c)
                            h=t*a(i,iz,is)
                            z(i,iz)=z(i,iz)-h
                            z(i,is)=z(i,is)+h
                            ew(i,iz)=ew(i,iz)-h
                            ew(i,is)=ew(i,is)+h
                            a(i,iz,is)=zero
                            DO j=1,iz-1
                                g=a(i,j,iz)
                                h=a(i,j,is)
                                a(i,j,iz)=g-s*(h+g*tau)
                                a(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            DO j=iz+1,is-1
                                g=a(i,iz,j)
                                h=a(i,j,is)
                                a(i,iz,j)=g-s*(h+g*tau)
                                a(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            DO j=is+1,dim2
                                g=a(i,iz,j)
                                h=a(i,is,j)
                                a(i,iz,j)=g-s*(h+g*tau)
                                a(i,is,j)=h+s*(g-h*tau)
                            ENDDO
                            DO j=1,dim2
                                g=ev(i,j,iz)
                                h=ev(i,j,is)
                                ev(i,j,iz)=g-s*(h+g*tau)
                                ev(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            nrot(i)=nrot(i)+1                        
                        ENDIF
                    ENDDO
                ENDDO
            ENDDO
            ! ---------------------
            ! update b/ew/z
            DO iz=1,dim2
#include "vectorize.inc"                    
                DO ijk=1,nindx
                    i = indx(ijk)
                    b(i,iz)=b(i,iz)+z(i,iz)
                    ew(i,iz)=b(i,iz)
                    z(i,iz)=zero
                ENDDO
            ENDDO
            ! ---------------------
        ENDDO
 
        DEALLOCATE( b )
        DEALLOCATE( z )
 
        RETURN