Subroutine to transform the elemental entries into assembled CC vectors
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(MAT_SOLV), | intent(inout), | target | :: | mat | high level system type |
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