module Eigen

! Modul odpovidajici pascalske UNITe tehoz jmena

use std_type

implicit none
private

public :: Jacobi, Tred, Tql


contains


! ===================== jacobi =================================    

subroutine Jacobi (dimA, A, V, etap)
    integer,intent(in) :: dimA
    real(kind=DP),dimension(:,:),intent(inout) :: A
    real(kind=DP),dimension(:,:),intent(out) :: V
    integer,intent(inout) :: etap
    intrinsic :: abs,sqrt

    real(kind=DP),parameter :: proMR = 100.0_DP
    integer :: i,j,MaxEtap,m,k
    real(kind=DP) :: p,q,t,c,s
    
    V = 0.0_DP           !jednotkova matice do V
    do i = 1,dimA
        V(i,i) = 1.0_DP
    end do      
    
    MaxEtap = etap      ! uchovani vstupu pro testy
    etap = 0            ! nastaveni pocitadla etap

    etapa: do              !start etapy
        etap = etap+1
        m = 0               ! pocitadlo preskocenych rotaci v etape
        cykl_i: do i = 1, dimA-1
            cykl_j: do j = i+1,dimA
                p = (A(i,j)+A(j,i))/2.0_DP
                    ! aritm. stred nebot nasobeni rotacni matici zleva
                    !  a zprava je vlivem zaokrouhlovacich chyb nekomutativni
                q = A(i,i)-A(j,j)
                t = sqrt(4.0_DP*p*p+q*q)
                if (t == 0.0_DP) then
                    m = m+1     ! preskocit rotaci
                    cycle cykl_j     ! dalsi j
                end if
                if (q<0) then   ! rotaci nutno provest, aby se usporadaly
                                ! diag. prvky a tim i vlastni hodnoty
                    s = sqrt((t-q)/(2*t))
                    if (p<0) s = -s
                    c = p/(t*s)
                else
                    if ((abs(A(i,i)) == (abs(A(i,i))+proMR*abs(p))).and.  &
                        (abs(A(j,j)) == (abs(A(j,j))+proMR*abs(p)))) then
                            ! faktor proMR se vybere tak, aby se provedla
                            ! "mala, ale ne prilis mala" rotace
                        m = m+1     ! preskocit rotaci
                        cycle cykl_j      ! jit na dalsi j
                    else            ! vypocet uhlu rotace, c=cos,s=sin
                        c = sqrt((t+q)/(2*t))
                        s = p/(t*c)
                    end if
                end if
                if (1+abs(s)==1) then   ! mala rotace, preskocit
                    m = m+1
                    cycle cykl_j        ! dalsi j
                end if
                
                do k = 1,dimA       ! rotace, nasobeni matice zleva
                    q = A(i,k)
                    A(i,k) = c*q+s*A(j,k)
                    A(j,k) = -s*q+c*A(j,k)
                end do
                do k = 1,dimA       ! nasobeni matice zprava
                    q = A(k,i)
                    A(k,i) = c*q+s*A(k,j)
                    A(k,j) = -s*q+c*A(k,j)
                end do
                do k = 1,dimA       ! vypocitat vlastni vektory
                    q = V(k,i)      ! nasobeni V zprava
                    V(k,i) = c*q+s*V(k,j)
                    V(k,j) = -s*q+c*V(k,j)
                end do
            end do cykl_j      
        end do cykl_i    
        
        if ((etap == MaxEtap) .or. (m==dimA*(dimA-1)/2)) return
        
    end do etapa        
     
end subroutine Jacobi

! ========================== tred ==============================
subroutine Tred (dimA,A1,d,e)
    integer,intent(in) :: dimA
    real(kind=DP),dimension(:),intent(inout) :: A1
    real(kind=DP),dimension(:),intent(out) :: d,e
    intrinsic :: epsilon,tiny,sqrt
    
    real(kind=DP),parameter :: tol = 10.0_DP*tiny(1.0_DP)/epsilon(1.0_DP)
    
    integer :: i,iz,j,jk,k,m
    real(kind=DP) :: f,g,h,hh
    
    cykl_i: do i = dimA,1,-1
        m = i-1
        iz = (i*m)/2
        h = 0
        do k = 1,m
            f = A1(iz+k)
            d(k) = f
            h = h+f*f
        end do
        if (h <= tol) then
            e(i) = 0        ! Householderova transformace se neprovadi, je-li
            h = 0           ! h tak male, ze nelze zarucit ortogonalitu
        else
            g = sqrt(h)
            if (f >= 0) g = -g
            e(i) = g
            h = h-f*g
            d(m) = f-g
            A1(iz+m) = d(m)
            f = 0
            do j = 1,m
                g = 0
                jk = (j*(j-1))/2+1
                do k = 1,m
                    g = g+A1(jk)*d(k)
                    if (k<j) then
                        jk = jk+1
                    else
                        jk = jk+k
                    end if
                end do
                g = g/h
                e(j) = g
                f = f+g*d(j)
            end do
            hh = f/(h+h)    ! vypocet parametru K
            jk = 0      
                
            do j = 1,m      ! transformace matice A
                f = d(j)
                g = e(j)-hh*f
                e(j) = g
                do k = 1,j
                    jk = jk+1
                    A1(jk) = A1(jk)-f*e(k)-g*d(k)
                end do
            end do
            
        end if      ! pripad h>tol
        
        d(i) = A1(iz+i)
        A1(iz+i) = h
    end do cykl_i
    return

end subroutine Tred

!========================= tql ==========================
subroutine Tql(dimA,d,e,iteraci)
    integer,intent(in) :: dimA
    real(kind=DP),dimension(:),intent(inout) :: d
    real(kind=DP),dimension(:),intent(inout) :: e
    integer,intent(inout),optional :: iteraci
    intrinsic :: abs,epsilon,present,maxval
    
    integer,parameter :: MaxIteraci = 30
    integer :: MaxIter,iter
    integer :: i,k,m
    real(kind=DP) :: b,c,f,g,h,p,r,s
    integer,dimension(dimA) :: aiter
    
    iF (present(iteraci)) then
        MaxIter = iteraci
    else
        MaxIter = MaxIteraci
    end if
    
    e(1:dimA-1) = e(2:dimA)
    e(dimA) = 0
    b = 0
    f = 0
    
    cykl_k: DO k = 1,dimA
        iter = 0
        h = epsilon(1.0_DP)*(abs(d(k))+abs(e(k)))
        if (b<h) b = h      ! zanedbatelne maly poddiagonalni prvek
        m = k
        do
            if ((abs(e(m)) <= b) .and. (m >= dimA)) exit
            m = m+1
        end do
        if (m>k) then
            iterace: do
                if (iter == MaxIter) then
                    if (present(iteraci)) iteraci = MaxIter
                    return
                end if
                iter = iter+1
                g = d(k)
                p = (d(k+1)-g)/(2*e(k))
                r = sqrt(p*p+1)
                if (p <= 0) then
                    d(k) = e(k)/(p-r)
                else
                    d(k) = e(k)/(p+r)
                end if
                h = g-d(k)
                d(k+1:dimA) = d(k+1:dimA)-h
                f = f+h
                p = d(m)
                c = 1
                s = 0
                do i = m-1,k,-1
                    g = c*e(i)
                    h = c*p
                    if (abs(p) >= abs(e(i))) then
                        c = e(i)/p
                        r = sqrt(c*c+1)
                        e(i+1) = s*p*r
                        s = c/r
                        c = 1/r
                    else
                        c = p/e(i)
                        r = sqrt(c*c+1)
                        e(i+1) = s*e(i)*r
                        s = 1/r
                        c = c/r
                    end if
                    p = c*d(i)-s*g
                    d(i+1) = h + s*(c*g+s*d(i))
                end do
                e(k) = s*p
                d(k) = c*p
                
                if (abs(e(k)) <= b) exit iterace
            end do iterace
            
        end if  ! m>k
        
        p = d(k)+f
        
        do i = k,2,-1       ! usporadani vlastnich hodnot vzestupne
            if (p < d(i-1)) then
                d(i) = d(i-1)
            else
                exit
            end if
        end do
        if (i==2) i = 1
        d(i) = p
        aiter(k) = iter
    end do cykl_k
    
    if (present(iteraci)) iteraci = maxval(aiter)
    return
    
end subroutine Tql


end module Eigen

