lab22b

module methods
    implicit none
    contains

    function f(x) result(y) !Функция исходная(с пограшностями) !ОК
        real(8), intent(in) ::  x
        real(8) :: y
        !y = x**2d0 + 1d0 - (4d0*datan(1d0)/2d0 - x - (x**3d0)/6d0)
        y = x**2d0 + 1d0 - dacos(x)
        !y = x**2d0
    end function

    function e(x,n) result(y) !Базисная функция !ОК
        integer, intent(in) :: n
        real(8), intent(in) ::  x
        real(8) :: y
        y = x**real(n-1,8)
        !y = dcos(real(n-1,8)*dacos(x))
    end function

    function pf(res,x,n) result(y) !Вычисление приближения в точке !ОК
        integer, intent(in) :: n
        real(8), intent(in) ::  x
        integer :: i
        real(8), dimension(:) :: res(n)
        real(8) :: y
        y=0d0
        do i = 1,n
            y = y + res(i)*e(x,i)
        end do
    end function
    subroutine cholesky1(A,b,n,res) !Разложение Холецкого
        real(8),dimension(:,:) :: A
        real(8),dimension(:) :: b, res
        real(8), allocatable :: y1(:)
        real(8) :: sum
        integer, intent(in) :: n
        integer :: i,j,P
        integer :: info
        allocate(y1(n))
        !L = 0d0
        !L = A
        !do i = 1, n !i столбец
            !do j = i, n !j строка
                !Вычисление суммы
                !sum = 0d0
                !do p = 1, i-1
                    !sum = sum + A(i,p)*A(j,p)
                !end do
                !Разложение
                !if (j==i)  then
                    !A(i,i) = dsqrt(A(i,i) - sum)
                !else
                    !A(j,i) = (A(j,i) - sum)/A(i,i)
                !end if
            !end do
        !end do
        call DPOTRF("L", n, A, n,info)
        write(*,*) "info", info
        do i = 1,n-1
            A(i,i+1:) = 0d0
        end do
        !Решаем Ly = b
        y1 = 0d0
        do i = 1,n
            do j = 1, i-1
                b(i) = b(i) - A(i,j)*y1(j)
            end do
            y1(i) = b(i)/A(i,i)
        end do
        !Решаем (L**T)x = y
        A = transpose(A)
        res = 0d0
        do i = n,1,-1
            do j = n, i+1,-1
                y1(i) = y1(i) - A(i,j)*res(j)
            end do
            res(i) = y1(i)/A(i,i)
        end do
        !write(*,*) "RESULT", norm2(MATMUL(A,res) - b)
        deallocate(y1)
    end subroutine

    subroutine gengrid(grid,n,a,b) !Генерация равномерной сетки !ОК
        integer, intent(in) :: n
        real(8),intent(in) :: a, b
        real(8), dimension(:,:) :: grid(n,4)
        real(8) :: AA
        integer :: i
        !Генерация сеток
        grid = 0d0
        do i = 1,n
            grid(i,1) = (b-a)*real(i,8)/real(n-1,8) + a - (b-a)/real(n-1,8)
            grid(i,2) = f(grid(i,1))
        end do
        grid(:,3) = grid(:,2)

        !Выбросы(промахи)
        !AA = dabs(maxval(grid(:,2)) - minval(grid(:,2))) !Амплитуда
        !grid(1,3) = grid(1,3)*2d0*AA
        !grid(2,3) = grid(2,3)*3d0*AA
        !grid(4,3) = grid(4,3)*4d0*AA
    end subroutine

    subroutine MNK(grid,w,n,m,res) !метод наименьших квадратов с конкректным весами.
        integer, intent(in) :: n, m !n - число узлов m - степень строящегося полинома
        integer :: i,j,k
        real(8), dimension(:,:) :: grid(n,4)
        real(8), dimension(:) :: w(n), res(n) !ВЕСА и решение
        real(8), allocatable :: A(:,:), b(:)
        allocate(A(m,m),b(m))
        !Составление матрицы СЛАУ с весами
        b = 0d0
        A = 0d0
        do k = 1, m
            do i = 1,n
                b(k) = b(k) + w(i)*grid(i,3)*e(grid(i,1),k)
            end do
            do j = 1,m
                do i = 1,n
                    A(k,j) = A(k,j) + w(i)*e(grid(i,1),j)*e(grid(i,1),k)
                end do
            end do
        end do
        !write(*,"(3f14.7)") (A(i,:),i=1,m)
        call cholesky1(A,b,m,res) !Решение СЛАУ
        deallocate(A, b)
    end subroutine

    subroutine coefget(grid,res,m,n)
        integer, intent(in) :: n, m !n - число узлов m - степень строящегося полинома
        integer :: i,j,k
        real(8), dimension(:,:) :: grid(n,4)
        real(8), dimension(:) :: w(n), res(n) !ВЕСА и решение
        real(8), allocatable :: epsi(:)
        real(8) :: eps
        allocate(epsi(n))
        eps = 1e-14

        w = 1d0
        do while(maxval(epsi) > eps)
            res = 0d0
            call MNK(grid,w,n,m,res)
            do i = 1,n
                epsi(i) = dabs(grid(i,3) - pf(res,grid(i,3),n))
                w(i) = 1d0/(epsi(i)*2d0 + eps)
            end do
        end do
        deallocate(epsi)
    end subroutine
end module