! $Id: utils_mod.f90 190 2011-02-08 19:21:37Z pamestoy $
!     Tools:
!       _parameters,
!       _Matrix Market parser,
!       _wrappers for ordering, factorization and solve,
!       _cleaning.

!     Remark: every object (matrix, permutation, etc.) is provided
!     to the user with "Fortran subscripts" (i.e. from 1 to N).
!     LDL uses C subscripts, thus every input of an LDL routine is
!     modified before call to LDL (my_array(i) = my_array(i) - 1),
!     and every ouput of an LDL routine is modified after a call to
!     LDL (my_array(i) = my_array(i) + 1).

module utils
  implicit none

  type csc_mat
     integer                       :: n  ! n:  dimension 
     integer                       :: nz ! nz: number of non-zeros in the lower part

     integer, allocatable          :: jptr(:)  ! pointers on columns
     integer, allocatable          :: irn(:)   ! row indices of non-zeros
     real(kind(1.d0)), allocatable :: val(:)   ! values of non-zeros
  end type csc_mat


contains


  subroutine read_parameters(mat_file, want_amd, itmax, eps, prec, level, trace)
    ! Read parameters file

    character*72, intent(out)         :: mat_file
    integer, intent(out)              :: want_amd, itmax, prec, level, trace
    real(kind(1.d0))                  :: eps

    integer :: i, allocok

    write(6,*)'Reading parameters...'
    open(4,file='parameters', status='old', action='read')
    read(4,*)mat_file
    read(4,*)want_amd
    read(4,*)itmax
    read(4,*)eps
    read(4,*)prec
    read(4,*)level
    read(4,*)trace
    close(4)
    return

  end subroutine read_parameters

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


  subroutine csc_read_matrix(mat_file, A, Ierror, Trace)
    !     Purpose:
    !      Reads a Matrix Market file (lower part), and converts to
    !      full matrix (upper + lower part) in Compressed column format
    !      (CSC). (CSC is the format used for Harwell Boeing files)
    !      Note that on output entries in the columns are not sorted.
    !     Input parameters:
    !      mat_file : file name for input matrix in matrix market format
    !      Trace    : level of printing
    !     Output parameters:
    !      Full matrix (lower and upper parts) in Compressed column format
    !      N  : integer, order of matrix
    !      NZ : number of non zeros in matrix
    !      Ap : array of size N+1 of integers, pointer in Ai and Ax to beginning of column i
    !      Ai : row subscripts
    !      Ax : numerical values
    !      Ierror : Integer holding the error return
    !      = 0   No error
    !      = 1   Pb when reading the input file
    !      = 2   Pb during memeory allocation
    !      = 3   Out of range indices (<0 or >N) in input matrix

    character*72, intent(in) :: mat_file
    integer, intent(in)      :: trace
    type(csc_mat)            :: A
    integer, intent(out)     :: ierror

    character*72 :: rep
    integer      :: i, iz, j, k, nerr, jpos, nzmm

    integer                       :: allocok
    integer, allocatable          :: mm_ai(:), mm_aj(:)
    double precision, allocatable :: mm_ax(:)

    integer, parameter :: infile=10     ! logical unit to read matrix
    integer, parameter :: mp=6          ! logical unit to write messages

    if (mp.ge.0) write(mp,*)'Reading matrix... '
    ierror = 0

    ! Read file with matrix in matrix market MM) format

    ! Open and then read all comments in the header of the input MM file
    open(Infile,file=mat_file, status='old', action='read', err=100)
    read(infile,*)rep
    do
       if(rep(1:1) .ne. '%') exit
       read(infile,*)rep
    end do
    backspace(infile)

    ! Read sizes and perform allocations
    read(infile,*,err=200) A%n,A%n,nzmm
    allocate(mm_ai(nzmm), mm_aj(nzmm), mm_ax(nzmm), stat=allocok)
    if (allocok.ne.0) goto 300

    ! Read file in Matrix Market format
    do iz = 1, nzmm
       read(infile,*, err=200) mm_ai(iz), mm_aj(iz), mm_ax(iz)
    end do
    close(infile)


    ! Convert triangular matrix in MM format to
    ! Full (upper and lower parts) matrix in CSC format
    ! (does not detect dupplicates)

    ! Allocate matrix in CSC format
    allocate(A%jptr(A%n+1), stat=allocok)
    if (allocok.ne.0) goto 300

    ! Scan the input, count the elements in each column in Ap
    ! and ignore elements with indices out of range.
    nerr            = 0    ! number of out of range indices
    A%jptr(1:A%n+1) = 0    ! Initialize all pointers
    A%nz            = nzmm ! We need to count twice off-diag entries
    do k=1,nzmm
       i =mm_ai(k)
       j =mm_aj(k)
       if ( (i.le.0 .or. i.gt.A%n) .or. (j.le.0 .or. j.gt.A%n) .or. (j.gt.i) ) then
          ! Out of range row or column index
          nerr = nerr+1
          cycle               ! Skip entry
       endif
       ! Count entry in column J
       A%jptr(j) = A%jptr(j) +1
    end do

    allocate(A%irn(A%nz),A%val(A%nz), stat=allocok)
    if (allocok.ne.0) goto 300

    if ((nerr.gt.0).and.(mp.ge.0)) then
       write(mp,*) '  Error : Out of range entry in orginal MM matrix '
       write(mp,*) nerr, '  Out of range indices in matrice '
       goto 400
    end if

    ! Prepare the row pointers in such a way that pointer A%jptr(J)
    ! points to the first position of column J+1.
    A%jptr(1)=A%jptr(1)+1
    do k=2,A%n+1
       A%jptr(k)= A%jptr(k)+A%jptr(k-1)
    end do

    ! Set the column pointers and store the matrix by column in one pass
    do k=1,nzmm
       j          = mm_aj(k)
       i          = mm_ai(k)
       jpos       = A%jptr(j)
       ! Store new entry in col J position JPOS-1
       A%irn(jpos-1) = i
       A%val(jpos-1) = mm_ax(k)
       A%jptr(j)     = jpos-1   ! Update pointer of column J
    end do

    if ((trace.eq.1).and.(mp.gt.0)) then
       write(mp,'(a,i9)') '   Matrix of size', A%n
       write(mp,'(a,i9)') '   Number of entries in a :', A%nz
    end if
    if ((trace.eq.2).and.(mp.gt.0)) then
       write(mp,*) '  Input: Triangular matrix in mm format:'
       write(mp,'(a,i9)') '   Matrix of size                      :', A%n
       write(mp,'(a,i9)') '   Number of entries in the lower part :', nzmm
       write(mp,*) '  row  col    value '
       do k=1,nzmm
          write(mp,'(i5, i5, f12.5)') mm_ai(k), mm_aj(k), mm_ax(k)
       enddo
       write(mp,*)'  Output: Full matrix in csc format:'
       write(mp,'(a,i9)') '   Matrix of size                      :', A%n
       write(mp,'(a,i9)') '   Number of entries in A (lower part) :', A%nz
9999   format (a,10i6/(8x, 10i6))
9998   format (a,5f12.5/(8x, 5f12.5))
       write(mp,9999) '   A%jptr  =', A%jptr
       write(mp,9999) '   A%irn   =', A%irn
       write(mp,9998) '   A%val   =', A%val
    end if

    goto 500  ! skip error processing

    ! error processing
100 continue
    ierror = 4
    goto 500
200 continue
    ierror = 1
    goto 500
300 continue
    ierror = 2
    goto 500
400 continue
    ierror = 3
500 continue

    ! free allocated memory
    if (allocated(mm_ai)) deallocate(mm_ai)
    if (allocated(mm_aj)) deallocate(mm_aj)
    if (allocated(mm_ax)) deallocate(mm_ax)

    return

  end subroutine csc_read_matrix

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


  subroutine csc_ordering(want_amd, A, trace)
    ! Ordering: AMD if requested, identity otherwise
    ! The permutation is applied explicitely to the matrix
    ! At the end, each column is sorted.

    integer, intent(in)          :: want_amd, trace
    type(csc_mat)                :: A


    integer, allocatable   :: P(:), Pinv(:), tmpptr(:), tmpirn(:), l(:)
    real(kind(1.d0)), allocatable :: tmpval(:)
    integer                :: amd_order, check_amd, amd_control, amd_info, allocok
    integer                :: i, j, jpos
    parameter(amd_control=5, amd_info=20)
    double precision       :: control(amd_control), info(amd_info)
    real                   :: t0, tf


    write(6,*)'Ordering...'
    call cpu_time(t0)
    ! Fortran to C subscripts to call AMD and LDL
    A%irn = A%irn-1
    A%jptr = A%jptr-1

    ! Ordering: AMD or nothing (identity)
    allocate(P(A%n), stat=allocok)
    if (allocok.ne.0) then
     write(*,*)'Problem with allocation of P in ordering', allocok
     stop
    endif
    if(want_amd.eq.1) then
      check_amd = amd_order (%val(A%N),A%jptr, A%irn, P, control, info)
      if((check_amd.ne.0).and.(check_amd.ne.1)) then
        write(6,*)'  Ordering with AMD failed !'
        stop
      else
        write(6,*)'  Ordering with AMD: OK.'
      end if
    else
      do i=1,A%N
        P(i)=i-1
      end do
      write(6,*)'  Ordering: identity.'
    end if

    ! C to Fortran subscripts for the user
    A%jptr = A%jptr+1
    A%irn = A%irn+1
    P = P+1

    call cpu_time(tf)
    write(6,'(A,F7.3,A)')'   Time for ordering =    ',tf-t0,' s'
    if(trace.eq.2) then
      9999 format (A,10I6/(8X, 10I6))
      9998 format (a,5f12.5/(8x, 5f12.5))
      write(6,9999)'   P   = ',P
    end if

    ! P is applied explicitely to the matrix. We need its inverse.

    write(*,*)'Permuting...'
    allocate(Pinv(A%n), stat=allocok)
    if (allocok.ne.0) then
      write(*,*)'Problem with allocation of Pinv in ordering', allocok
      stop
    endif
    do i=1,A%n ! Pinv
      Pinv(P(i)) = i
    end do
    do i=1,A%nz ! row indices
      A%irn(i)= Pinv(A%irn(i))
    end do

    ! the matrix is not triangular anymore, so we need to move elements
    ! which are in the upper part to the lower part, recompute the column pointers
    ! and rebuild the CSC as in read_matrix.
    allocate(tmpptr(A%n+1),tmpirn(A%nz),tmpval(A%nz), stat=allocok)
    if (allocok.ne.0) then
      write(*,*)'Problem with allocation of tmp arrays in ordering', allocok
      stop
    endif
    tmpptr = 0
    !do i=1,A%n ! first, tmptr will contain the number of elements per column: same as in A (sym perm)
    !  tmpptr(i) = A%jptr(i+1)-A%jptr(i)
    !end do
    do i=1,A%n ! count elements in each column
      do j=A%jptr(i),A%jptr(i+1)-1
        if(A%irn(j).lt.Pinv(i)) then ! element in the upper part: moved to the lower part (col irn(j))
          tmpptr(A%irn(j)) = tmpptr(A%irn(j))+1
        else ! element is in the lower part (col P(i))
          tmpptr(Pinv(i)) = tmpptr(Pinv(i))+1
        end if
      end do
    end do
    tmpptr(1)=tmpptr(1)+1 ! prepare pointers so that they point onto the first element of the next column
    do i=2,A%n+1
       tmpptr(i)= tmpptr(i)+tmpptr(i-1)
    end do
    tmpptr(A%n+1) = A%nz+1
    do i=1,A%n ! scan elements and put them at the right place
      do j=A%jptr(i),A%jptr(i+1)-1
        if(A%irn(j).lt.Pinv(i)) then ! element in the upper part: moved to the lower part (col irn(j))
          jpos             = tmpptr(A%irn(j))
          tmpirn(jpos-1)   = Pinv(i)
          tmpval(jpos-1)   = A%val(j)
          tmpptr(A%irn(j)) = jpos- 1
        else ! element is in the lower part (col P(i))
          jpos            = tmpptr(Pinv(i))
          tmpirn(jpos-1)  = A%irn(j)
          tmpval(jpos-1)  = A%val(j)
          tmpptr(Pinv(i)) = jpos- 1
        end if
      end do
    end do
    A%jptr = tmpptr
    A%irn  = tmpirn
    A%val  = tmpval
    deallocate(tmpptr, tmpirn, tmpval, P, Pinv)
    !write(*,9999) '   A%jptr  =', A%jptr
    !write(*,9999) '   A%irn   =', A%irn
    !write(*,9998) '   A%val   =', A%val

    ! Sort each column
    write(*,*)'Sorting...'
    allocate(l(A%n+2), stat=allocok)
    if (allocok.ne.0) then
      write(*,*)'Problem with allocation of l in ordering', allocok
      stop
    endif
    do i=1,A%n
      l = 0
      call mergesort(A%jptr(i+1)-A%jptr(i),  A%irn(A%jptr(i):A%jptr(i+1)-1), l, 1)
      call mergeswap(A%jptr(i+1)-A%jptr(i),l,A%irn(A%jptr(i):A%jptr(i+1)-1), A%val(A%jptr(i):A%jptr(i+1)-1))
    end do
    if(trace.eq.2) then
      write(*,9999) '   A%jptr  =', A%jptr
      write(*,9999) '   A%irn   =', A%irn
      write(*,9998) '   A%val   =', A%val
    end if
    deallocate(l)

    return

  end subroutine csc_ordering

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

  subroutine mergeswap(n, l, a1, a2)

    integer          :: i, lp, iswap, n
    integer          :: l(0:n+1), a1(n)
    real(kind(1.d0)) :: a2(n), dswap

    lp = l(0)
    i  = 1
    do
       if ((lp .eq. 0).or.(i>n)) exit
       do
          if (lp >= i) exit
          lp = l(lp)
       end do
       iswap    = a1(lp)
       a1(lp)   = a1(i)
       a1(i)    = iswap
       dswap    = a2(lp)
       a2(lp)   = a2(i)
       a2(i)    = dswap
       iswap    = l(lp)
       l(lp) = l(i)
       l(i)  = lp
       lp = iswap
       i  = i + 1
    enddo

    return

  end subroutine mergeswap

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

  subroutine mergesort(n, k, l, order)

    !      Plain implementation of the merge-sort algorithm
    !      as described in:

    !      D. E. Knuth "The Art of Computer Programming,"
    !      vol.3: Sorting and Searching, Addison-Wesley, 1973
    implicit none
    integer           :: n
    integer, intent(inout)           :: k(n), l(0:n+1)
    integer, optional :: order

    integer    :: p, q, s, t, i, iord

    iord = 1
    if(present(order)) iord = order

    if((iord .ne. 1) .and. (iord .ne. -1)) then
       write(*,'("Wrong input in mergesort")')
       return
    end if

100 continue
    l(0) = 1
    t = n + 1
    do  p = 1,n - 1
       if (iord*k(p) .le. iord*k(p+1)) then
          l(p) = p + 1
       else
          l(t) = - (p+1)
          t = p
       end if
    end do
    l(t) = 0
    l(n) = 0

    if (l(n+1) .eq. 0) then
       return
    else
       l(n+1) = iabs(l(n+1))
    end if

200 continue
    s = 0
    t = n+1
    p = l(s)
    q = l(t)

    if(q .eq. 0) return

300 continue
    if(iord*k(p) .gt. iord*k(q)) goto 600

400 continue
    l(s) = sign(p,l(s))
    s = p
    p = l(p)
    if (p .gt. 0) goto 300

500 continue
    l(s) = q
    s = t
    do
       t = q
       q = l(q)
       if (q .le. 0) exit
    end do
    goto 800

600 continue
    l(s) = sign(q, l(s))
    s = q
    q = l(q)
    if (q .gt. 0) goto 300

700 continue
    l(s) = p
    s = t
    do
       t = p
       p = l(p)
       if (p .le. 0) exit
    end do

800 continue
    p = -p
    q = -q
    if(q.eq.0) then
       l(s) = sign(p, l(s))
       l(t) = 0
       goto 200
    end if

    goto 300

    return

  end subroutine mergesort

end module utils

