! ************************************************************************
! Fájl: QRAlg.f08
! Leírás: QR felbontás és QR módszer mátrix sajátértékeinek meghatározására
! ************************************************************************

program QRAlg
    implicit none
    ! Rövid főprogram amely teszt adatokon futtatja a QR módszer algoritmusát
    double precision, dimension(4,4) :: A,eig
    double precision :: eps = 0.00001 ! epsilon a megengedett hiba

    ! Az A 4x4-es mátrix, elemei oszlopfolytonosan megadva
    A = reshape( [2.0, 2.0, 0.0, 0.0, 1.0, 3.0, 0.0, 0.0, -1.0, 1.0, 4.0, -1.0, 1.0, -2.0, 2.0, 1.0], shape(A) )
    ! eig egy 4x4-es felsőháromszögmátrix,
    ! a diagonálisában a számított sajátértékekkel
    eig = qr_method(A, eps)

    ! Input és az eredmény kiíratása
    call print_matrix(A)
    print *, '--------------------------'
    call print_matrix(eig)

    ! A lényegi munkát végző függvények
    contains

        ! ********************************************************
        ! Egy adott mátrix QR felbontását készíti el
        ! Input: A = négyzetes valós elemű mátrix
        ! Output: Q,R = négyzetes valós mátrixok,
        ! A = Q*R, Q ortogonális, R felsőháromszög mátrix
        ! A megvalósítás Gram-Schmidt ortogonalizációval történik
        ! Segédfüggvények: normalize(), projection()
        ! ********************************************************
        subroutine qr_decomposition(A,Q,R)
            implicit none
            double precision, intent(in), dimension(:,:) :: A
            double precision, intent(out), dimension(:,:) :: Q,R

            ! segédváltozó, az A n x n-es esetén egy n hosszú vektor
            double precision, dimension(size(A,1)) :: sum_proj
            integer :: n,k,j
            n = size(A,1) ! mátrix mérete

            ! Elég csupán a Q mátrixot elkészíteni, az R utána szorzással megkapható
            ! Az algoritmus inicializálása, az input mátrix első oszlopát normalizáljuk

            ! Megjegyzés: A(1:n,1) syntaxis: array slice, A első oszlopát adja
            Q(1:n,1) = normalize(A(1:n,1))
            ! A ciklusban Q további oszlopait ortogonalizáljuk minden előzőre
            do k=2,n
                ! Gram-Schmidt féle projekciók elkészítése és összegzése
                sum_proj = 0.0
                do j=1,k-1
                    ! A k. lépésben az eredeti mátrix k. oszlopát vetítjük a már meglévőkre
                    sum_proj = sum_proj + projection(vector=A(1:n,k),project_to=Q(1:n,j))
                end do
                ! A vetítések összegét levonva az eredeti vektorból,
                ! egy az összes a rendszerben őt megelőzőre merőleges vektort kapunk,
                ! majd az elkészült új oszlopot normalizáljuk
                Q(1:n, k) = normalize(A(1:n, k) - sum_proj)
            end do

            ! R = Q'*A, a két intrinsic a transzponálást és mátrixszorzást végzi
            R = matmul(transpose(Q), A)

        end subroutine qr_decomposition

        ! ********************************************************
        ! Merőleges vetület számítása
        ! Input: vector = n hosszú valós vektor, amit vetítünk
        ! project_to, ugyanakkora méretű vektor, amire vetítünk
        ! Result: vetület vektor
        ! ********************************************************
        function projection(vector, project_to)
            implicit none
            double precision, intent(in), dimension(:) :: vector, project_to
            double precision, dimension(size(vector)) :: projection
            ! proj a on b = <a,b>/<a,a> * a
            ! a vetület vektor tehát két skalár szorzat hányadosa * a vektor amire vetítünk
            projection = dot_product(project_to, vector)/dot_product(project_to, project_to) * project_to
        end function projection

        ! ********************************************************
        ! Vektor normalizálása (2-es normában)
        ! Input: vector, n hosszú valós vektor
        ! Result: a normalizált vektor
        ! ********************************************************
        function normalize(vector)
            implicit none
            double precision, intent(in), dimension(:) :: vector
            double precision, dimension(size(vector)) :: normalize
            ! A vektort elosztjuk a hosszával, vagyis sqrt(<v,v>) - vel
            normalize = vector / sqrt(dot_product(vector, vector))
        end function normalize

        ! ********************************************************
        ! Mátrix sajátértékeinek meghatározása QR módszerrel
        ! Input: matrix, négyzetes valós mátrix
        ! eps: numerikus hibakorlát
        ! Result: R felsőháromszög mátrix,
        ! ami a sajátértékeket a diagonálisában tartalmazza
        ! X(k+1) = R(k) * Q(k) k=1,2,... az iteráció
        ! segédfüggvény: qr_decomposition(),check_convergence()
        ! ********************************************************
        function qr_method(matrix, eps)
            implicit none
            double precision, intent(in), dimension(:,:) :: matrix
            double precision, dimension(size(matrix,1),size(matrix,2)) :: qr_method
            double precision, intent(in) :: eps
            double precision, dimension(size(matrix,1),size(matrix,2)) :: Q,R,eigenvalues

            ! inicializáljuk az algoritmust
            eigenvalues = matrix
            ! A fő iterációs ciklus
            do while ( check_convergence(eigenvalues, eps) .eqv. .false. )
                ! Előállítjuk a QR felbontást,
                call qr_decomposition(eigenvalues, Q,R)
                ! majd fordított sorrendben szorozzuk őket össze
                eigenvalues = matmul(R,Q)
            end do
            qr_method = eigenvalues
        end function qr_method

        ! ********************************************************
        ! A QR módszer segéd függvénye a konvergencia ellenőrzésére
        ! Input: matrix: négyzetes valós mátrix
        ! eps: a numerikus hibakorlát
        ! ********************************************************
        function check_convergence(matrix, eps)
            implicit none
            logical :: check_convergence
            double precision, intent(in), dimension(:,:) :: matrix
            double precision, intent(in) :: eps
            double precision :: max_err
            integer :: n,i,j
            n = size(matrix,1)
            ! A mátrixnak felsőháromszög mátrixhoz kell tartania,
            ! vagyis az alsóháromszög beli elemeknek 0-hoz kell tartania

            ! Megkeressük az alsóháromszög rész abszolútértékben maximális elemét
            max_err = 0.0
            do i=2,n
                do j=1,i-1
                    if ( abs(matrix(i,j)) > max_err ) then
                        max_err = abs(matrix(i,j))
                    end if
                end do
            end do

            ! Ha a maximális hiba kisebb mint epsilon akkor ok
            if ( max_err < eps ) then
                check_convergence = .true.
            else
                check_convergence = .false.
            end if
        end function check_convergence


        ! Egyszerű kiírató segédfüggvény, a mátrixot soronként írja ki
        subroutine print_matrix(matrix)
            implicit none
            double precision, dimension(:,:) :: matrix
            integer :: i,j,n,m
            n = size(matrix,1)
            m = size(matrix,2)
            do i=1,n
                print *, ( matrix(i,j),  j=1,m )
            end do

        end subroutine print_matrix

end program QRAlg