csr_mod.F90

Go to the documentation of this file.

module csr_mod

  use real_type
  use algebra_set, only: shellsort_dec
  use basic_tools
  !$ use OMP_LIB
 
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  !! TYPE
  public :: csr

  !! GENERIC
  public :: clear
  public :: write
  public :: print
  public :: valid          !! TESTED
  public :: transpose      !! TESTED
  public :: matvecprod     !! TESTED
  public :: setrow         !! TESTED
  public :: getrow         !! TESTED
  public :: addentry       !! TESTED
  public :: equal
  public :: copy 
  public :: add
  public :: scale
  public :: reorder  
  public :: sortrows
  public :: invllt

  !! NON GENERIC
  public :: addentry_2    !! TESTED
  public :: getdiag
  public :: add_submatrix
  
  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%

  type csr
     integer :: nl     = 0
     integer :: nnz    = 0
     logical :: sorted = .false.
     integer , dimension(:), allocatable :: row 
     integer , dimension(:), allocatable :: col 
     integer :: width  = 0
     real(RP), dimension(:), allocatable :: a   
     integer , dimension(:), allocatable :: diag

   contains
     
     final :: csr_clear

  end type csr

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

  !! constructor
  interface csr
     module procedure csr_create
  end interface csr

  interface print
     module procedure csr_print
  end interface print

  interface valid
     module procedure csr_valid        !! TESTED
  end interface valid

  interface getrow
     module procedure csr_getrow       !! TESTED
  end interface getrow

  interface setrow
     module procedure csr_setrow       !! TESTED
  end interface setrow

  interface sortrows
     module procedure csr_sortrows
  end interface sortrows

  interface addentry
     module procedure csr_addentry      !! TESTED
  end interface addentry

  interface equal
     module procedure csr_equal
  end interface equal

  interface write
     module procedure csr_write
  end interface write

  interface matvecprod
     module procedure csr_matvecprod    !! TESTED
  end interface matvecprod

  interface add
     module procedure csr_add_1, csr_add_2
  end interface add

  interface add_submatrix
     module procedure add_submatrix_1, add_submatrix_2
  end interface add_submatrix

  interface scale
     module procedure csr_scale
  end interface scale

  interface reorder
     module procedure csr_reorder
  end interface reorder

  interface copy
     module procedure csr_copy
  end interface copy

  interface transpose
     module procedure csr_transpose  !! TESTED
  end interface transpose

  interface invllt
     module procedure csr_invllt
  end interface invllt

contains

  subroutine csr_clear(m)
    type(csr), intent(inout) :: m

    call freemem(m%row)
    call freemem(m%col)
    call freemem(m%a  )
    call freemem(m%diag)
    m%nl     = 0
    m%nnz    = 0
    m%sorted = .false.
    m%width  = 0

  end subroutine csr_clear

  function csr_create(nnz) result(m)
    type(csr)                            :: m
    integer, dimension(:), intent(in)    :: nnz

    integer :: nz, ii, jj, wdt, length

    call clear(m)
    m%nl  = size(nnz,1)
    nz    = sum(nnz)
    m%nnz = nz
    wdt   = 0
    call allocmem(m%row,m%nl+1)
    call allocmem(m%col, nz   )
    call allocmem(m%a  , nz   )

    m%row = 0
    m%col = 0
    m%a   = 0._rp

    ! Construction of m%row
    jj = 1
    do ii = 1, m%nl
       m%row(ii) = jj
       length = nnz(ii)
       if (length>wdt) wdt = length
       jj = jj + length
    end do
    m%row(m%nl+1) = jj

    m%width = wdt
    if (.NOT.valid(m)) call quit( &
         & 'csr_mod: csr_create')

  end function csr_create


  function csr_valid(m) result(b)
    type(csr), intent(in) :: m
    logical :: b

    integer :: ii, wdt, length

    b = .false.
    if (.NOT. (allocated(m%row)) ) return
    if (.NOT. (allocated(m%row)) ) return
    if (.NOT. (allocated(m%a))   ) return

    b =  ( m%nl>0 ) .AND. ( m%nnz>0 )
    b =  b .AND. ( size(m%row,1) == m%nl+1 )
    b =  b .AND. ( size(m%col,1) == m%nnz  )
    b =  b .AND. ( size(m%a  ,1) == m%nnz  )
    b =  b .AND. ( m%row(m%nl+1) == size(m%col)+1 )

    if (allocated(m%diag)) then
       b = b.AND.(size(m%diag,1)  == m%nl)
       b = b.AND.( maxval(m%diag) <= m%nnz )
    end if

    wdt = 0
    do ii = 1, m%nl
       length = m%row(ii+1) - m%row(ii)
       if (length>wdt) wdt = length
    end do
    b = ( wdt == m%width )

  end function csr_valid



  subroutine csr_getrow(sz, v, p, n, g, ii)
    real(RP)   , dimension(n), intent(out):: v
    integer    , dimension(n), intent(out):: p
    integer                  , intent(out):: sz
    type(csr)                , intent(in) :: g
    integer                  , intent(in) :: ii, n

    integer :: j1, j2

#if DBG>0
    if ( .NOT. valid(g) )     call quit( &
         & "csr_mod: csr_getRow_1: CSR matrix 'g' not valid")

    if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
         & "csr_mod: csr_getRow_1: wrong line index 'ii'")
#endif

    j1 = g%row(ii)
    j2 = g%row(ii+1)
    sz = j2 - j1

#if DBG>0
    if (sz>n )  call quit( &
         & "csr_mod: csr_getRow_1: unsufficient size&
         &  for the array 'p'")
#endif

    p(1:sz) = g%col(j1:j2-1)
    v(1:sz) = g%a(j1:j2-1)
    
  end subroutine csr_getrow

  subroutine csr_setrow(m, val, col, ii)
    type(csr)             , intent(inout):: m
    real(RP), dimension(:), intent(in)   :: val
    integer , dimension(:), intent(in)   :: col
    integer               , intent(in)   :: ii

    integer :: j1, j2

    j1 = m%row(ii)
    j2 = m%row(ii+1)
    m%col(j1:j2-1) = col
    m%a(  j1:j2-1) = val

  end subroutine csr_setrow


  function csr_equal(g1, g2) result(b)
    logical               :: b
    type(csr), intent(in) :: g1, g2

    integer , dimension(g1%width) :: p1, p2, q1, q2
    real(RP), dimension(g1%width) :: a1, a2
    integer :: ii, jj, sz, wdt

    b = valid(g1) .AND. valid(g2)
    b = b .AND. ( g1%nl == g2%nl ) .AND. ( g1%nnz == g2%nnz )
    if(.NOT.b) return

    b = b .AND. ( g1%width == g2%width )
    if(.NOT.b) return

    b = b .AND. (all( g1%row == g2%row ) )
    if(.NOT.b) return

    wdt = g1%width
    do ii=1, g1%nl
       call getrow(sz, a1, p1, wdt, g1, ii)
       call getrow(sz, a2, p2, wdt, g2, ii)

       if (sz==0) cycle

       call shellsort_dec(p1(1:sz), q1(1:sz), sz)
       call shellsort_dec(p2(1:sz), q2(1:sz), sz)

       do jj=1, sz
          b = b .AND. ( p1(q1(jj)) == p2(q2(jj)) ) 
          b = b .AND. equal(a1(q1(jj)), a2(q2(jj))) 
       end do

       if (.NOT. b) return

    end do 

  end function csr_equal


  subroutine csr_addentry(m, rw, cl, val)
    type(csr), intent(inout) :: m
    integer , intent(in)    :: rw, cl
    real(RP), intent(in)    :: val

    integer :: ii

    do ii = m%row(rw),m%row(rw+1)-1

       if (m%col(ii)==0) then
          ! the entry m(rw,cl) does not exist yet
          m%col(ii)= cl
          m%a(ii)  = val
          goto 10
       else if (m%col(ii)==cl) then
          ! the entry m(rw,cl) already exists
          m%a(ii) = m%a(ii)+val
          goto 10
       end if
    end do
    ! the line rw is already filled
    ! write(*,*) m%col(m%row(rw):m%row(rw+1)-1) 
    ! write(*,*) cl 
    call quit("csr_mod: csr_addEntry: 1")

10  continue
  end subroutine csr_addentry


  subroutine addentry_2(m, rw, cl, val)
    type(csr), intent(inout) :: m
    integer , intent(in)    :: rw, cl
    real(RP), intent(in)    :: val

    integer :: ii

    ii = m%row(rw)

    do while(.true.)
       if (m%col(ii)==cl) then
          m%a(ii) = m%a(ii)+val
          return
       end if
       ii = ii + 1
    end do

  end subroutine addentry_2



  subroutine csr_write (m, fileName, fileFormat)
    type(csr)       , intent(in) :: m
    character(len=*), intent(in) :: fileName, fileFormat

    write(*,*)"csr_mod   : csr_write   = "//trim(filename) 

    select case(trim(fileformat))
    case("matlab")
       call matlab_write (m, filename)

    case default
       call quit("csr_mod: csr_write: unknown fileFormat")
    end select

  end subroutine csr_write

  subroutine matlab_write (m, fileName)
    type(csr)       , intent(in) :: m
    character(len=*), intent(in) :: fileName

    integer :: ii, jj, col, nnz

#if DBG>0
    if (.NOT.valid(m)) call quit( &
         & 'csr_mod: matlab_write: 1')      
#endif

    write(*,*)"  Marlab ASCII format"

    nnz = size(m%a,1)
    open(unit = 999,file = filename)
    do ii = 1,m%nl
       do jj=m%row(ii), m%row(ii+1)-1
          col = m%col(jj)
          if (col==0) exit
          write(999,*) ii, col, m%a(jj) 
       end do
    end do
    close(999)

  end subroutine matlab_write



  subroutine csr_copy(s2,s1)
    type(csr), intent(inout) :: s2
    type(csr), intent(in)    :: s1

#if DBG>0
    if (.NOT.valid(s1)) call quit( &
         & 'csr_mod: csr_copy: 1')      
#endif

    call clear(s2)

    s2%nl        = s1%nl
    s2%nnz       = s1%nnz
    s2%sorted    = s1%sorted
    
    call allocmem(s2%row, s2%nl+1)
    call allocmem(s2%col, s2%nnz )
    call allocmem(s2%a  , s2%nnz )
    s2%row   = s1%row
    s2%col   = s1%col
    s2%a     = s1%a  
    s2%width = s1%width

    if (allocated(s1%diag)) then
       call allocmem(s2%diag, s2%nl)
       s2%diag = s1%diag
    end if

  end subroutine csr_copy



  subroutine csr_matvecprod(v,m,u)
    real(RP), dimension(:), intent(out) :: v
    type(csr)             , intent(in)  :: m
    real(RP), dimension(:), intent(in)  :: u

    real(RP) :: entry
    integer  :: ii, jj

#if DBG>0
    if (.NOT.valid(m)) call quit( &
         & 'csr_mod: csr_matVecProd: CSR matrix not valid')      
    if (size(v,1) /= m%nl) call quit(&
         & "csr_mod: csr_matVecProd: output vectr size uncorrect")
    ! if (size(u,1) /= m%nc ) call quit(&
    !      & "csr_mod: csr_matVecProd: intput vectr size uncorrect2")
#endif

    !$OMP PARALLEL PRIVATE(entry) 
    !$OMP DO
    do ii = 1,m%nl
       !  m%row(ii)     = line ii storage strarting index
       !  m%row(ii+1)-1 = line ii storage ending    index

       entry = 0._rp

       do jj = m%row(ii), m%row(ii+1)-1
          entry = entry + m%a(jj)*u(m%col(jj))
       end do          

       v(ii) = entry

    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine csr_matvecprod


  subroutine csr_add_2(m, m1, a, m2, b)
    type(csr), intent(inout) :: m
    type(csr), intent(in)    :: m1, m2
    real(RP) , intent(in)    :: a, b

#if DBG>0      

    !! checks if patterns are identical

    logical    :: b1, b2

    if (.NOT.valid(m1)) call quit( &
         & "csr_mod: csr_add_2: 'm1' not valid")      
    if (.NOT.valid(m2)) call quit( &
         & "csr_mod: csr_add_2: m2' not valid")      

    b1 = (size(m1%row,1) == size(m2%row,1))
    b2 = (size(m1%col,1) == size(m2%col,1))
    if (b1) b1 = (all(m1%row==m2%row))
    if (b2) b2 = (all(m1%col==m2%col))
    b1 = (b1.AND.b2)
    if (.NOT.b1) call quit("csr_mod: csr_add_2:&
         & 'm1' and 'm2' have different patternsx")
#endif

    call copy(m, m1)
    m%a = a*m1%a + b*m2%a

  end subroutine csr_add_2



  subroutine csr_add_1(m1, a, m2, b)
    type(csr) , intent(inout) :: m1
    type(csr) , intent(in)    :: m2
    real(RP), intent(in)    :: a, b

#if DBG>0       

    !! checks if patterns are identical

    logical    :: b1, b2

    if (.NOT.valid(m1)) call quit( &
         & 'csr_mod: csr_add_1: 1')      
    if (.NOT.valid(m2)) call quit( &
         & 'csr_mod: csr_add_1: 2')      

    b1 = (size(m1%row,1) == size(m2%row,1))
    b2 = (size(m1%col,1) == size(m2%col,1))
    if (b1) b1 = (all(m1%row==m2%row))
    if (b2) b2 = (all(m1%col==m2%col))
    b1 = (b1.AND.b2)
    if (.NOT.b1) call quit('csr_mod : csr_add_3: 1')
#endif

    m1%a = a*m1%a + b*m2%a

  end subroutine csr_add_1


  subroutine csr_scale(m1, a)
    type(csr), intent(inout) :: m1
    real(RP) , intent(in)    :: a

#if DBG>0
    if (.NOT.valid(m1)) call quit( &
         & 'csr_mod: csr_scale: 1')      
#endif

    m1%a = m1%a *  a

  end subroutine csr_scale


  !       %----------------------------------------%
  !       |                                        |
  !       |          PRECONDITIONNING              |
  !       |                                        |
  !       %----------------------------------------%


  subroutine csr_reorder(s, perm, permInv)
    type(csr)            , intent(inout) :: s
    integer, dimension(:), intent(in)    :: perm, permInv

    type(csr)    :: s2
    integer      :: ii, jj, jDeb, jFin, kk, cpt

#if DBG>0
    if (.NOT.valid(s)) call quit( &
         & 'csr_mod: csr_reorder: 1')
    if (size(perm, 1)/=s%nl) call quit( &
         & 'csr_mod: csr_reorder: 2')
    if (size(perminv, 1)/=s%nl) call quit( &
         & 'csr_mod: csr_reorder: 3')
#endif

    call copy(s2, s)

    ! s2 <= reordering of s
    cpt = 1
    do ii=1, s%nl

       kk = perm(ii)
       s2%row(ii) = cpt

       jdeb = s%row(kk)
       jfin = s%row(kk+1) - 1

       do jj = jdeb,jfin  

          kk = s%col(jj)

          s2%col(cpt) = perminv(kk) 
          s2%a(cpt)   = s%a(jj)
 
          cpt = cpt + 1
       end do
    end do

    call copy(s, s2)

  end subroutine csr_reorder



  subroutine csr_sortrows(s)
    type(csr), intent(inout) :: s

    integer :: ll, j1, j2, sz, ii
    integer , dimension(s%width) :: new_col,old_col
    real(RP), dimension(s%width) :: old_a

    if (.NOT.valid(s)) call quit( &
         & "csr_mod: csr_sortRows: csr matrix 's' not valid")      

    do ll = 1,s%nl
       j1 = s%row(ll)  
       j2 = s%row(ll+1)
       sz = j2 - j1
       if (sz<=1) cycle

       j2 = j2 - 1
       old_col(1:sz) = s%col(j1:j2) 
       old_a(1:sz)   = s%a(j1:j2) 
       call shellsort_dec(old_col(1:sz), new_col(1:sz), sz) 

       do ii = 0, sz-1
          s%col(j1+ii) = old_col(new_col(sz-ii))
          s%a  (j1+ii) = old_a(new_col(sz-ii))
       end do

    end do

    s%sorted = .true.

  end subroutine csr_sortrows


  subroutine csr_invllt(x, LLT, b)
    type(csr)             , intent(in)  :: LLT
    real(RP), dimension(:), intent(in)  :: b
    real(RP), dimension(:), intent(out) :: x

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

#if DBG>0
    if (.NOT.valid(llt)) call quit( &
         & 'csr_mod: csr_invllt: 1')      
    if (size(b,1)/=llt%nl) call quit( &
         & 'csr_mod: csr_invllt: 2')      
    if (size(x,1)/=llt%nl) call quit( &
         & 'csr_mod: csr_invllt: 3')      
#endif
    
    call allocmem(z, llt%nl)

    call invl(z, llt, b)
    call invlt(x, llt, z)

  end subroutine csr_invllt


  subroutine invlt(x, LT, b)
    type(csr)               , intent(in)  :: LT
    real(RP), dimension(:), intent(in)  :: b
    real(RP), dimension(:), intent(out) :: x

    integer :: ii,jj, jDeb,jFin

#if DBG>0
    if (.NOT.valid(lt)) call quit( &
         & 'csr_mod: invLT: 1')      
    if (size(b,1)/=lt%nl) call quit( &
         & 'csr_mod: invLT: 2')      
    if (size(x,1)/=lt%nl) call quit( &
         & 'csr_mod: invLT: 3')      
#endif

    ! solve LT * x = b
    x(lt%nl) = b(lt%nl) / lt%a(lt%nnz) 
    do ii = lt%nl-1,1,-1

       jdeb = lt%diag(ii)
       jfin = lt%row(ii+1)-1

       x(ii) = b(ii)
       do jj = jdeb+1, jfin 
          x(ii) = x(ii) - lt%a(jj) * x(lt%col(jj))
       end do

       x(ii) = x(ii) / lt%a(jdeb)
    end do

  end subroutine invlt

  subroutine invl(x,L,b)
    type(csr)               , intent(in)  :: L
    real(RP), dimension(:), intent(in)  :: b
    real(RP), dimension(:), intent(out) :: x

    integer :: ii,jj, jDeb,jFin

#if DBG>0
    if (.NOT.valid(l)) call quit( &
         & 'csr_mod: invL: 1')      
    if (size(b,1)/=l%nl) call quit( &
         & 'csr_mod: invL: 2')      
    if (size(x,1)/=l%nl) call quit( &
         & 'csr_mod: invL: 3')      
#endif

    ! solve L * x = b
    x(1) = b(1) / l%a(1)                       
    do ii = 2, l%nl                            
       x(ii) = b(ii)                             
       jdeb = l%row(ii)
       jfin = l%diag(ii)

       do jj = jdeb, jfin-1                       
          x(ii) = x(ii) - l%a(jj) * x(l%col(jj)) 
       end do

       x(ii) = x(ii) / l%a(jfin)

    end do

  end subroutine invl


  subroutine getdiag(A)
    type(csr), intent(inout) :: A

    integer :: i, k, kDeb, kFin

#if DBG>0
    if (.NOT.valid(a)) call quit( &
         & 'csr_mod: getDiag: 1')      
#endif

    !! check if the diagonal has already been set:
    !! If so exit
    if (allocated( a%diag) ) then
       if (size(a%diag,1) ==  a%nl ) return
    end if

    !! initialisation
    call allocmem(a%diag, a%nl)
    a%diag = -1  

    do i = 1,a%nl
       kdeb = a%row(i)
       kfin = a%row(i+1)-1
       do k = kdeb,kfin
          if (a%col(k)==i) then 
             a%diag(i) = k
             exit 
          end if
       end do
    end do

  end subroutine getdiag


  subroutine csr_transpose(g2, g1)
    type(csr), intent(inout)  :: g2
    type(csr), intent(in)     :: g1

    integer, dimension(:), allocatable :: nnz
    integer :: nc, jj, ii, j1, j2

#if DBG>0
    if (.NOT.valid(g1)) call quit( &
         & 'csr_mod: csr_transpose: 1')      
#endif

    call clear(g2)

    nc = maxval(g1%col)

    call allocmem(nnz, nc)
    nnz = 0

    do ii=1, g1%nnz
       jj = g1%col(ii)
       nnz(jj) = nnz(jj) + 1
    end do
    g2 = csr(nnz)
    call freemem(nnz)

    do ii=1, g1%nl

       j1 = g1%row(ii)
       j2 = g1%row(ii+1)-1

       do jj=j1, j2
          call addentry(g2, g1%col(jj), ii, g1%a(jj))
       end do

    end do

  end subroutine csr_transpose


  subroutine add_submatrix_1(m, m2)
    type(csr) , intent(inout) :: m
    type(csr) , intent(in)    :: m2

    integer :: j1, j2, ii, jj

#if DBG>0
    if (.NOT.valid(m))  call quit( &
         & 'csr_mod: add_subMatrix1: 1')      
    if (.NOT.valid(m2)) call quit( &
         & 'csr_mod: add_subMatrix1: 2')      
#endif

    do ii=1, m2%nl

       j1 = m2%row(ii)
       j2 = m2%row(ii+1)

       do jj=j1, j2-1
          call addentry(m, ii, m2%col(jj), m2%a(jj))
       end do

    end do

  end subroutine add_submatrix_1


  subroutine add_submatrix_2(m, m2, a)
    type(csr) , intent(inout) :: m
    type(csr) , intent(in)    :: m2
    real(RP)  , intent(in)    :: a

    integer :: j1, j2, ii, jj

#if DBG>0
    if (.NOT.valid(m))  call quit( &
         & 'csr_mod: add_subMatrix1: 1')      
    if (.NOT.valid(m2)) call quit( &
         & 'csr_mod: add_subMatrix1: 2')      
#endif

    do ii=1, m2%nl

       j1 = m2%row(ii)
       j2 = m2%row(ii+1)

       do jj=j1, j2-1
          call addentry(m, ii, m2%col(jj), a*m2%a(jj))
       end do

    end do

  end subroutine add_submatrix_2



  subroutine csr_print(g)
    type(csr), intent(in) :: g

    write(*,*) 'csr_mod    : print'

    write(*,*) '  Valid                 = ', valid(g)
    write(*,*) '  Nb lines              =' , g%nl
    write(*,*) '  Nb non-zero entries   =' , g%nnz
    write(*,*) '  Sorted lines          = ', g%sorted
    write(*,*) '  Maximal line width    =' , g%width

  end subroutine csr_print

end module csr_mod