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) !! Subroutine to transform the elemental entries into assembled CC vectors 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_array2(tab_inout = jcol(1:ntot), & ! tab1 = irow(1:ntot), & ! tab2 = a_elt(1:ntot), n = 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_array2(tab_inout = irow(i1:i2), & ! tab1 = a_elt(i1:i2), n = i2 - i1 +1 ) ! 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