precond_mod.F90

Go to the documentation of this file.

module precond_mod

  use choral_constants
  use choral_variables
  use real_type
  use algebra_set, only: shellsort_dec
  use basic_tools
  use csr_mod
 
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: precond
  public :: clear
  public :: jacobi, icc0, reverse_cuthill_mackey

  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type precond

     integer :: type

     real(RP), dimension(:), allocatable :: invd

     type(csr) :: llt

     integer , dimension(:), allocatable :: perm

     integer , dimension(:), allocatable :: perminv

     character(20) :: name

   contains
     
     final :: precond_clear

  end type precond
  
  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface precond
     module procedure precond_create
  end interface precond

  interface clear
     module procedure precond_clear
  end interface clear

contains

  subroutine precond_clear(pc)
    type(precond), intent(inout)   :: pc

    call clear(pc%LLT) 
    call freemem(pc%invD)
    call freemem(pc%perm)
    call freemem(pc%permInv)

  end subroutine precond_clear


  function precond_create(K, type) result(pc)
    type(precond)            :: pc
    type(csr), intent(inout) :: K
    integer  , intent(in)    :: type

    call clear(pc)
    if(.NOT.valid(k)) call quit("precond_mod: precond_create:&
         & matrix 'K' not valid")

    pc%type = type

    !! Name definition
    !!
    select case(type)

    case(pc_0)
       pc%name = "void"

    case(pc_jacobi)
       pc%name = "Jacobi"

    case(pc_icc0)
       pc%name = "iCC0"

    case default
       call quit( 'precond_mod: precond_create:&
            & preconditioning type incorrect')

    end select

    if (choral_verb>0) write(*,*) &
         &'precond_mod     : create         = ',&
         & trim(pc%name)

    !! assembling
    !!
    select case(type)

    case(pc_jacobi)
       call jacobi(pc%invD, k)

    case(pc_icc0)
       call reverse_cuthill_mackey(pc%perm, pc%permInv, k)
       call reorder(k, pc%perm, pc%permInv)
       call icc0(pc%LLT, k)

    end select

  end function precond_create

  subroutine  jacobi(invD, s)
    real(RP) , dimension(:), allocatable   :: invD
    type(csr)              , intent(inout) :: s

    integer :: ii, jj

    if (choral_verb>1) write(*,*) &
         & 'precond_mod     : jacobi'
    call getdiag(s)

    !! check diagonal entries
#if DBG>0
    if (minval(s%diag) <= 0 ) call quit( &
         & 'precond_mod: jacobi: negative diagonal entries')
#endif
    
    call allocmem(invd, s%nl)
    do ii=1, s%nl      
       jj = s%diag(ii)
       invd(ii) = 1._rp / s%a(jj)
    end do

  end subroutine jacobi



  subroutine reverse_cuthill_mackey(perm, permInv, s)   
    type(csr)              , intent(in)  :: s
    integer  , dimension(:), allocatable :: perm, permInv

    integer, dimension(s%nl) :: deg, P, LS, LS0, degLS, new_i
    logical, dimension(s%nl) :: marked
    integer  :: ii, jj, jDeb, jFin, kk
    integer  :: nL, nL0, node

    if (choral_verb>1) write(*,*) &
         & 'precond_mod     : reverse_cuthill_macKey'

    call allocmem(perm, s%nl)
    call allocmem(perminv, s%nl)

    ! Initialisation
    marked = .false.

    ! Degre des sommets
    do ii = 1,s%nl
       deg(ii) = s%row(ii+1)-s%row(ii)
    end do

    ! Initialisation de l'algo
    kk   = 1
    node = 1
    jj = minval(deg)
    do ii=1, s%nl
       if (deg(ii)==jj) then
          node = ii                  ! Premier noeud, degre minimal
          exit
       end if
    end do
    marked(node) = .true.            ! Marque

    nl0    = 1                       ! taille de LS0
    ls0(1) = node                    ! Premiere Level-set, 1 elt

    p(kk) = node                     ! Mise à jour de la permutation


    do while ( kk<s%nl )
       ! Construction de la level-set suivant LS0

       nl = 0 
       do ii = 1, nl0 

          jdeb = s%row( ls0(ii)   )
          jfin = s%row( ls0(ii)+1 ) - 1

          do jj = jdeb,jfin  

             node = s%col(jj)

             if (.not.marked(node)) then
 
                nl           = nl+1
                ls(nl)       = node
                marked(node) = .true.     ! Marquage du noeud

             end if

          end do
       end do

       ! actualisation NL0, LS0
       nl0       = nl
       ls0(1:nl) = ls(1:nl)

       do ii = 1,nl
          degls(ii) = deg(ls(ii))         ! degLS := degre(LS)
       end do

       ! On parcourt LS par ordre croissant de degre
       call shellsort_dec(degls(1:nl), new_i(1:nl), nl)    ! Tri ordre decroissant
       do ii = nl,1,-1
          kk    = kk + 1                   ! Noeud suivant
          p(kk) = ls(new_i(ii))            ! MAJ permutation
       end do

    end do

    do ii = 1,s%nl
       perm(ii) = p(s%nl+1-ii)
    end do

    do ii = 1, s%nl
       perminv(perm(ii)) = ii
    end do

  end subroutine reverse_cuthill_mackey



  subroutine  icc0(LLT, s)
    type(csr), intent(out)    :: LLT
    type(csr), intent(inout)  :: s

    integer :: i,k,j,m,mm,band, kDeb,kFin
    real(RP) :: int1,int2,sum, t1, t2

    if (choral_verb>1) write(*,'(A35,$)') &
         &' precond_mod     : icc0           ='
    call cpu_time(t1)

    if (.not.s%sorted) call sortrows(s)
    call getdiag(s)
    call copy(llt, s)

    ! Calcul de la largeur de bande
    band = 0
    do i = 1,s%nl
       band = max( band, s%col(s%row(i+1)-1)-i )
    end do
    if (choral_verb>1) write(*,*) &
         & int(band,4), ' :  band width'


    ! Calcul de la decomposition de Choleski LL^T incomplete. C'est a
    ! dire que l'on ne rempli L (et L^T) que la ou S avait un
    ! coefficient non nul. Dans l'algo ci-dessous on raisonne par
    ! ligne, c'est donc la partie L^T (triangulaire superieure) qui
    ! est remplie en premier.
    do  i = 1,s%nl
       do m = s%diag(i),s%row(i+1)-1 ! Termes sup-diagonaux
          j = s%col(m) !               j>=i
          sum = s%a(m)
          do  k = max(i-band-1,1),i-1 ! Lignes precedentes, k = 1 a i-1
             int1 = 0._rp
             int2 = 0._rp
             do  mm = s%diag(k),s%row(k+1)-1 ! On ne cherche que parmi les sup-diagonaux
                if  (s%col(mm)==i) int1 = llt%a(mm) ! C'est LT_{ki}, k<i
                if  (s%col(mm)==j) int2 = llt%a(mm) ! C'est LT_{kj}, k<i<=j
             end do
             sum = sum - int1*int2
          end do
          if (j==i) then
             llt%a(m) = sqrt(sum)
          else
             llt%a(m) = sum / llt%a(s%diag(i))
          end if
       end do
    end do

    ! On recopie les coefficients dans la partie triangulaire
    ! inferieure pour construire L.
    do i = 1,llt%nl-1          ! Parcours des lignes 1 a N-1
       kdeb = llt%diag(i)+1
       kfin = llt%row(i+1)-1
       do k = kdeb,kfin       ! Parcours des colonnes j de la ligne i avec i<j
          j = llt%col(k)
          ! Il faut recopier le coeff L_{ij} dans L_{ji}
          do m = llt%row(j), llt%diag(j)-1 ! Parcours de la ligne j colonnes m tq m<j
             if (llt%col(m)==i) then        ! On a trouve l'indice de la colonne _i_
                llt%a(m) = llt%a(k)         ! On fait LLT_{ji} := LLT_{ij}
                exit                        ! On sort de la boucle sur _m_
             end if
          end do
       end do
    end do

    call cpu_time(t2)

  end subroutine icc0


end module precond_mod