from_elemental_to_assembled Subroutine

public subroutine from_elemental_to_assembled(mat)

Arguments

Type IntentOptional AttributesName
type(MAT_SOLV), intent(inout), target:: mat

high level system type


Called by

proc~~from_elemental_to_assembled~~CalledByGraph proc~from_elemental_to_assembled from_elemental_to_assembled proc~convert_matrice_format convert_matrice_format proc~convert_matrice_format->proc~from_elemental_to_assembled program~test_solvers test_solvers program~test_solvers->proc~convert_matrice_format

Contents


Source Code

   subroutine from_elemental_to_assembled(mat)
   implicit none
   type(MAT_SOLV), intent(inout), target :: mat   !! *high level system type*
      integer(kind=I4), pointer :: solver, nb_elt, n, ntot, nz
      integer(kind=I4), dimension(:), pointer :: eltptr
      integer(kind=I4), dimension(:), pointer :: eltvar
      real(kind=R8),    dimension(:), pointer :: a_elt
      integer(kind=I4), dimension(:), pointer :: irow, jcol
      integer(kind=I4), dimension(:), pointer :: jptr
      integer(kind=I4) :: inelt, imatorder, i, j, ii, jj, i1, i2, ir1, ir2, jr1, jr2, itmp, innz, state

      solver   => mat%slv_t
      nb_elt   => mat%ne
      n        => mat%nn
      ntot     => mat%nt
      nz       => mat%nz

      if (solver==MUMP) return

      ! conversion from elemental form to triplet form, perhaps with null a_elts
      !--------------------------------------------------------------------------
      state = 0
      if (.not.allocated(mat%irow)) allocate( mat%irow(ntot), stat = state )
      if (state/=0) stop 'Memory allocation problem, FROM_ELEMENTAL_TO_ASSEMBLED'
      if (.not.allocated(mat%jcol)) allocate( mat%jcol(ntot), stat = state )
      if (state/=0) stop 'Memory allocation problem, FROM_ELEMENTAL_TO_ASSEMBLED'

      eltptr   => mat%eltptr
      eltvar   => mat%eltvar
      a_elt    => mat%a_elt
      irow     => mat%irow
      jcol     => mat%jcol
      jptr     => mat%jptr

      irow(1:ntot) = -1
      jcol(1:ntot) = -1

      ii = 1
      do inelt = 1, nb_elt
         imatorder = eltptr(inelt +1) -eltptr(inelt)
         do i = 1, imatorder
            irow(ii:(ii +imatorder -1)) = eltvar(eltptr(inelt):eltptr(inelt+1) -1)
            jcol(ii:(ii +imatorder -1)) = eltvar(eltptr(inelt) +i -1)
            ii = ii +imatorder
         enddo
      enddo

      ! where a_elt brings no contribution, rows and columns are zeroed
      !-----------------------------------------------------------------
      where( abs(a_elt) < EPS_R8 )
         irow = 0
         jcol = 0
      endwhere

      ! the triplet irow, jcol and a_elt is sorted according jcol
      !-----------------------------------------------------------
      call sort_int_1int_1real(g=1, d=ntot, itabref=jcol(1:ntot), itab1=irow(1:ntot), rtab2=a_elt(1:ntot))

      ! column pointer determination for each new value
      !----------------------------------------------
      if (allocated(mat%jptr)) deallocate(mat%jptr) ; allocate( mat%jptr(n +1) ) ; jptr => mat%jptr
      ii = 1
      do
         if (jcol(ii) > 0) exit  ! zeroed terms are ignored
         ii = ii +1
      enddo
      jptr(1) = ii

      do i = 1, n -1
         itmp = jcol(ii)
         do
            ii = ii +1
            if (jcol(ii) /= itmp) exit
         enddo
         jptr(i +1) = ii
      enddo
      jptr(n +1) = ntot +1

      ! columns are already sorted; rows are now sorted for each row value
      !--------------------------------------------------------------------
      do i = 1, n
         i1 = jptr(i)
         i2 = jptr(i +1) -1
         call sort_int_1real(g=1, d=i2 -i1 +1, itabref=irow(i1:i2), rtab1=a_elt(i1:i2))
      enddo

      ! assembly starting from the jcol, irow and a_elt top
      ! for identical matrix locations, a_elts are added
      !-----------------------------------------------------
      innz = 1
      jj  = jptr(1)  ! first non-zero element
      jr1 = jcol(jj)
      ir1 = irow(jj)
      jcol(innz)  = jr1
      irow(innz)  = ir1
      a_elt(innz) = a_elt(jj)

      do j = jj +1, ntot
         jr2 = jcol(j)
         ir2 = irow(j)
         if ( (jr2 /= jr1).or.(ir2 /= ir1) ) then  ! if row or column index has changed
            innz = innz +1                         ! a non-zero term is added
            jcol(innz) = jr2
            irow(innz) = ir2
            a_elt(innz) = a_elt(j)
         else
            a_elt(innz) = a_elt(innz) + a_elt(j)   ! row and column indexes are the same, stiffness terms are added
         endif
         jr1 = jr2                                 ! stores (i-1) and (j-1) for further comparison
         ir1 = ir2
      enddo
      nz = innz
      jcol( nz +1:ntot) = -1
      irow( nz +1:ntot) = -1
      a_elt(nz +1:ntot) = huge(1._R8)

      ! col pointer update
      !----------------------------------------------
      jj      = 1
      jptr(1) = 1
      do j = 1, n -1
         itmp = jcol(jj)
         do
            jj = jj +1
            if (jcol(jj) /= itmp) exit
         enddo
         jptr(j +1) = jj
      enddo
      jptr(n +1) = nz +1

      nullify(eltptr, eltvar, a_elt, irow, jcol, jptr)

   return
   endsubroutine from_elemental_to_assembled