! $Id: icfacto_mod.f90 198 2011-03-24 14:27:20Z fhrouet $
! Compute an incomplete ("fill level") Cholesky factorization.

module icfacto
  use utils
  use graph
  use graph_facto
  implicit none

contains

! just for us: remove it before giving to students
  subroutine csc_example_facto(A, L, p, trace)
    ! p maximum fill level allowed
    implicit none
    type(csc_mat) :: A, L

    integer :: i, j, p, trace

    write(6,*)'Incomplete factorization...'
    if(trace.eq.2) then
      write(6,*)'  Max fill level allowed:',p
    end if

    allocate(L%irn(A%n), L%jptr(A%n+1), L%val(A%n))

    L%n  = A%n
    L%nz = A%n

    do j=1, A%n
       L%irn(j)=j
       L%jptr(j)=j
       do p= A%jptr(j), A%jptr(j+1)-1
          i = A%irn(p)
          if(i .eq. j) L%val(j) = sqrt(A%val(p))
       end do
    end do
    L%jptr(A%n+1)=A%n+1

    return
  end subroutine csc_example_facto
! /just for us

  !------------------------------------------------------------------------------

  subroutine csc_incomplete_facto(A, L, p, trace)
  ! ---------
  ! parameters
  ! ---------
  ! Input (not modified in routine):
    ! A     input matrix 
    ! p     maximum level of fill-in allowed
    ! trace level of printing
           ! 0 no trace, 1 to display some relevant scalar
           ! values (N, LNZ, ..), 2 to display vectors/matrices (Ap, Ai, Ax, ..).
  ! Output
    ! L Matrix of incomplete factors
    implicit none
    type(csc_mat), intent(in)  :: A
    type(csc_mat), intent(out) :: L
    type(cell), dimension(:), allocatable :: g
    integer, intent (in) :: p, trace
  !----------------
  ! Local variables
  !----------------

    integer :: i, j, md

    write(6,*)'Incomplete factorization...'
    if(trace.eq.2) then
      write(6,*)'  Max fill level allowed:',p
    end if

    ! create graph associated to matrix A
    call csc_to_graph(A, g, trace)

    ! incomplete factorization : modification of the graph
    ! to include fill-in up to level p
    call Chol_IC(g, A%n, p, md, trace)

    if(trace.eq.1) then
      write(6,*)md,' diagonal elements had to be modified.'
    end if

    ! save graph of the incomplete factors into matrix L
    call graph_to_csc(g, A%n, L, trace)

    ! graph is not needed anymore: destroy it.
    call destroy_graph(g, A%n)

    return
  end subroutine csc_incomplete_facto

  !------------------------------------------------------------------------------

end module icfacto
