mumps_mod.F90

Go to the documentation of this file.

module mumps_mod

  use choral_variables
  use choral_constants
  use real_type
  use basic_tools
  use csr_mod

  implicit none
  private

#ifdef WMUMPS
  include 'mpif.h'
  include 'dmumps_struc.h'
#endif

  public :: mumps_init
  public :: mumps, clear, solve

  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%  
  type mumps
     logical :: created = .false.

#ifdef WMUMPS
     type(dmumps_struc) :: dmumps
#endif

   contains

     final :: mumps_clear

  end type mumps

  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure mumps_clear
  end interface clear

  interface mumps
     module procedure mumps_create
  end interface mumps

  interface solve
     module procedure mumps_solve
  end interface solve

contains

  
  subroutine mumps_init()
#ifdef WMUMPS
    logical :: ierr

    if (choral_verb>1) write(*,*)"mumps_mod       : mumps_init"

    if (rp /= dp) call quit(&
         & "mumps_mod: mumps_init: double precision only")

    call mpi_init(ierr)
#endif   
  end subroutine mumps_init


  subroutine mumps_clear(lu)
    type(mumps), intent(inout) :: lu

#ifdef WMUMPS

    if (lu%created) then

       if (associated(lu%dmumps%irn)) deallocate(lu%dmumps%irn)
       lu%dmumps%irn => null()
       if (associated(lu%dmumps%jcn)) deallocate(lu%dmumps%jcn)
       lu%dmumps%jcn => null()
       if (associated(lu%dmumps%rhs)) deallocate(lu%dmumps%rhs)
       lu%dmumps%rhs => null()
       if (associated(lu%dmumps%a  )) deallocate(lu%dmumps%a)
       lu%dmumps%a => null()
       
       ! Destroy the instance (deallocate internal data structures)
       lu%dmumps%job = -2
       call dmumps(lu%dmumps)
    end if
    
#endif

    lu%created = .false.

  end subroutine mumps_clear




  function mumps_create(g, type, verb) result(lu)
    type(mumps)                     :: lu
    type(csr), intent(in)           :: g
    integer  , intent(in)           :: type
    integer  , intent(in), optional :: verb

    integer  :: ii, jj, cpt
    real(DP) :: cpu

    cpu = clock()

    if (choral_verb>0) write(*,*)"mumps_mod       : create"
    if(.NOT.valid(g)) call quit( &
         & "mumps_mod: mumps_create:  CSR matrix 'g' not valid" )
    call clear(lu)

#ifdef WMUMPS

    !!
    !! MUMPS : INITIALISATION
    !!
    if (choral_verb>1) write(*,*)"  mumps structure initialisation"
    lu%dmumps%comm = mpi_comm_world
    lu%created = .true.
    lu%dmumps%job = -1

    select case(type)
    case(mumps_lu)
       lu%dmumps%sym = 0
       if (choral_verb>1) write(*,*)"  Non symmetric CSR"

    case(mumps_ldlt)
       lu%dmumps%sym = 2
       if (choral_verb>1) write(*,*)"  Symmetric CSR"

    case(mumps_ldlt_sdp)
       lu%dmumps%sym = 1
       if (choral_verb>1) write(*,*)"  Symmetric positive definite CSR"

    case default
       call quit("mumps_mod: mumps_create: &
            & unknown parameter 'type' " )

    end select

    lu%dmumps%par = 1
    call dmumps(lu%dmumps)
    !!
    !! check error
    if (lu%dmumps%infog(1)<0) then
       if (choral_verb>0) then
          write(*,*) "  dmumps error, code            =", &
               & lu%dmumps%infog(1), lu%dmumps%infog(2)

       end if
       call quit("mumps_mod: mumps_create: init. error")
    end if
    !!    
    !! transfer the CSR matrix entries
    !!
    if (lu%dmumps%myid==0) then
       
       !! number of lines
       lu%dmumps%n   = g%nl

       !! symmetric case : count entries a_ij with j <= i only
       if (lu%dmumps%sym > 0) then
          cpt = 0
          do ii=1, g%nl          
             do jj=g%row(ii), g%row(ii+1) - 1
                if (g%col(jj) <= ii ) cpt = cpt + 1
             end do
          end do
          lu%dmumps%nnz = cpt

       else  !! non symmetric case, all entries are counted
          lu%dmumps%nnz = g%nnz
       end if

       allocate(lu%dmumps%rhs(g%nl ))
       allocate(lu%dmumps%irn(g%nnz))
       allocate(lu%dmumps%jcn(g%nnz))
       allocate(lu%dmumps%a(  g%nnz))
       
       !! fill in the entries
       cpt = 0
       if (lu%dmumps%sym > 0) then  !! symmetric case

          do ii=1, g%nl 
          
             do jj=g%row(ii), g%row(ii+1) - 1

                if (g%col(jj) <= ii ) then
                   cpt = cpt + 1
                   
                   lu%dmumps%irn(cpt) = ii 
                   lu%dmumps%jcn(cpt) = g%col(jj)
                   lu%dmumps%a(cpt)   = g%a(jj)
                   
                end if
                
             end do
          end do
          
       else  !! non symmetric case

          do ii=1, g%nl
             
             do jj=g%row(ii), g%row(ii+1) - 1
                cpt = cpt + 1
                
                lu%dmumps%irn(cpt) = ii 
                lu%dmumps%jcn(cpt) = g%col(jj)
                lu%dmumps%a(cpt)   = g%a(jj)
                
             end do
          end do
       
       end if
    end if

    !!
    !! MUMPS VERBOSITY
    !!
    lu%dmumps%icntl(4) = -1
    !!
    if (present(verb)) then
       lu%dmumps%icntl(4) = verb
    end if

    !!
    !! MUMPS : ANALYSIS PHASE
    !!
    if (choral_verb>1) write(*,*)"  mumps analysis phase"
    lu%dmumps%job = 1
    call dmumps(lu%dmumps)
    !!
    !! check error
    if (lu%dmumps%infog(1)<0) then
       if (choral_verb>0) then
          write(*,*) "  dmumps error, code            =", &
               & lu%dmumps%infog(1), lu%dmumps%infog(2)

       end if
       call quit("mumps_mod: mumps_create: analysis phase error")
    end if

    !!
    !! MUMPS : FACTORISATION PHASE
    !!
    if (choral_verb>1) write(*,*)"  mumps factorisation phase"
    lu%dmumps%job = 2
    call dmumps(lu%dmumps)
    !!
    !! check error
    if (lu%dmumps%infog(1)<0) then
       if (choral_verb>0) then
          write(*,*) "  dmumps error, code            =", &
               & lu%dmumps%infog(1), lu%dmumps%infog(2)

       end if
       call quit("mumps_mod: mumps_create: factorisation phase error")
    end if

#else
    call quit("mumps_mod: mumps_create: MUMPS not available")

#endif

    cpu = clock() - cpu
    if (choral_verb>1) write(*,*) &
         & "  CPU                           =", real(cpu, sp)

  end function mumps_create


  subroutine mumps_solve(x, LU, b)
    real(RP), dimension(:), intent(inout) :: x
    type(mumps)           , intent(inout) :: LU
    real(RP), dimension(:), intent(in)    :: b

#ifdef WMUMPS

    lu%dmumps%rhs = b
    lu%dmumps%job = 3
    call dmumps(lu%dmumps)
    !!
    !! check error
    if (lu%dmumps%infog(1)<0) then

       write(*,*) "mumps_mod       : mumps_solve   = error ", &
            & lu%dmumps%infog(1), lu%dmumps%infog(2)
       call quit("mumps_mod: mumps_solve; error detected")
    end if

    !!  retrieve solution
    x = lu%dmumps%rhs 

#else
    call quit("mumps_mod: mumps_solve: MUMPS not available")

#endif

  end subroutine mumps_solve

end module mumps_mod