program main
  implicit none
!=========================================================

    character(10) :: filename
    integer       :: n, i
    integer, allocatable :: Ap(:), Ai(:)
    real(8), allocatable :: Ax(:), right(:), result(:)

    integer iwk
    integer lfil

    real(8), allocatable :: alu(:), w(:), wk(:)
    integer, allocatable :: jlu(:), ju(:), levs(:), jw(:)

    integer :: ipar(16), ierr
    real(8) :: fpar(16)

    real(8) :: startpoint, finishpoint

!=========================================================

    read(*, *) filename
    open(1, file=filename)

    read(1, *) n

    allocate(Ap(n + 1))
    read(1,*) (Ap(i), i = 1, n + 1)

    allocate(Ai(Ap(n + 1) - 1))
    read(1,*) (Ai(i), i = 1, Ap(n + 1) - 1)

    allocate(Ax(Ap(n + 1) - 1))
    read(1,*) (Ax(i), i = 1, Ap(n + 1) - 1)

    allocate(right(n))
    do i = 1, n
        right(i) = sin(real(i, 8))
    end do

    allocate(result(n))
    do i = 1, n
        result(i) = 0
    end do

!=======================================================

    iwk = 10 * (Ap(n + 1) - 1)
    lfil = 1
    ierr = 0

    allocate(alu(iwk), jlu(iwk), ju(n), levs(iwk), w(n), jw(3 * n))

!========================================================

    call cpu_time(startpoint)
    call iluk(n, Ax, Ai, Ap, lfil, alu, jlu, ju, levs, iwk, w, jw, ierr)
    call cpu_time(finishpoint)

    if (ierr .eq. 0) then
        print*, 'Preconditioner finished successfully in ', finishpoint - startpoint, ' seconds'
    else
        print*, 'Preconditioner failed with status: ', ierr
    end if

!========================================================

    ipar(1) = 0
    ipar(2) = 1
    ipar(3) = 1
    ipar(4) = 7 * n
    ipar(5) = 10
    ipar(6) = 1000
    fpar(1) = 1.0E-6
    fpar(2) = 1.0E-8
    fpar(11) = 0.0

    allocate(wk(7 * n))

    call cpu_time(startpoint)
    10  call bcg(n, right, result, ipar, fpar, wk)

    if (ipar(1).eq.1) then
        call amux(n, wk(ipar(8)), wk(ipar(9)), Ax, Ai, Ap)
        goto 10
    else if (ipar(1).eq.2) then
        call atmux(n, wk(ipar(8)), wk(ipar(9)), Ax, Ai, Ap)
        goto 10
    else if (ipar(1).eq.3 .or. ipar(1).eq.5) then
        call lusol(n, wk(ipar(8)),wk(ipar(9)), alu, jlu, ju)
        goto 10
    else if (ipar(1).eq.4 .or. ipar(1).eq.6) then
        call lutsol(n,wk(ipar(8)),wk(ipar(9)), alu, jlu, ju)
        goto 10
    else if (ipar(1).le.0) then
        if (ipar(1).eq.0) then
            print *, 'Iterative solver finished successfully '
        else
            print *, 'Iterative solver failed with status: ', ipar(1)
        endif
    endif
    call cpu_time(finishpoint)

    print*, 'in ', finishpoint - startpoint, ' seconds'
    print*, 'Iteration number: ', ipar(7)

end

real(8) function distdot(n,x,ix,y,iy)
    implicit none
    real(8), external :: ddot
    integer n, ix, iy
    real(8) x(1+(n-1)*ix), y(1+(n-1)*iy)
    distdot = ddot(n, x, ix, y, iy)
end function