anchor.diff:        do i=clrstarts(clr),clrstarts(clr+1)-1
           clrnode=clrorder(i)
           do j=A%M_bound(clrnode),A%M_bound(clrnode+1)-1
  -          clrneigh=M%eptnmap(A%indi(j))
  +          node=A%indi(j)
  +          clrneigh=M%eptnmap(node)
             if (clrneigh/=clr) then
               helper(clrneigh)=1
  -            if (onfront(A%indi(j))==0) then
  +            interfsz=interfsz+1
  +            if (onfront(clrnode)==0) then
                 nfront=nfront+1
  -              front(nfront)=A%indi(j)
  -              onfront(A%indi(j))=1
  +              front(nfront)=clrnode
  +              onfront(clrnode)=-1
               endif
             endif
           enddo
  
positive.code:     endif
    end subroutine SpMtx_aggregate
  
    subroutine SpMtx_DistributeAssembled(A,M)
      use Graph_class
      use Mesh_class
      implicit none
  
      type(SpMtx),intent(inout) :: A
      type(Mesh)                :: M
      integer :: i,ierr
      integer, dimension(:), pointer :: xadj
      integer, dimension(:), pointer :: adjncy
      integer                        :: nedges
      type(Graph) :: G
      integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
      integer,dimension(:),pointer       :: partorder,partstarts
      integer,dimension(4)           :: buf
      
      if (ismaster()) then ! Here master simply splits the matrix into pieces
  
negative.code:     !  ghost values in Ax operation
      !if ol>0, then the subdomain expands in ol layers (for solves)
  
      !TODO: it would be good to build also communication structures here.
      !        1. communication for Ax operation
      !           -- only the first layer of ghost nodes get communicated
      !        2. communication after applying the subdomain solves in Prec.
      !           -- all the ghost values get communicated
      !    but: seems that now a few datastructures are in SpMtx_operation
      !           (fexchindx,pid2indx) and we cannot do it right away...
      !        I think these need to get moved to Mesh and 
      !          need better namings as well! TODO !
      !TODO: A_interf and A_ghost need to get assembled here
  
      !Now, let's go on with next layers...:
      !if ol==0 we still need 1st layer for Ax operation...
      lastlayer=max(ol,1)
      ghostsz=0
      if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
        do layer=1,lastlayer
          frontstart(layer)=nfront+1
          ! look through the neighs of last layer's nodes:
          do i=frontstart(layer-1),frontend(layer-1)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indj(j)
              if (layer>1.and.onfront(node)>=0) then
  write(stream,*)'ghost: adding first:',front(i),node
                ghostsz=ghostsz+1
              endif
              if (onfront(node)==0) then
                clrneigh=M%eptnmap(node)
                if (clrneigh/=clr) then
                  nfront=nfront+1
                  front(nfront)=node
                  onfront(node)=layer !ie. if 1 => ghostnode for Ax op.
                endif
              endif
            enddo                                  
          enddo
          frontend(layer)=nfront
        enddo
      elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
        do layer=1,lastlayer
          frontstart(layer)=nfront+1
          ! look through the neighs of last layer's nodes:
          do i=frontstart(layer-1),frontend(layer-1)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indi(j)
  !            if (onfront(node)>=0) then
  !              ghostsz=ghostsz+1
  !            endif
              if (onfront(node)==0) then
                clrneigh=M%eptnmap(node)
                if (clrneigh/=clr) then
                  nfront=nfront+1
                  front(nfront)=node
                  onfront(node)=layer !ie. if 1 => ghostnode for Ax op.
                endif
              endif
            enddo                                  
          enddo
          frontend(layer)=nfront
        enddo
      else
        call DOUG_abort('SpMtx_keep_subd_wol: Matrix arrangment not done!',19)
      endif
      ! Add to the ghostsize also the connections within the outermost layer:
      if (ol>0) then
        if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
          ! look through the neighs of last layer's nodes:
          do i=frontstart(lastlayer),frontend(lastlayer)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indj(j)
              if (onfront(node)>0) then
  write(stream,*)'ghost: adding also:',front(i),node
                ghostsz=ghostsz+1
              endif
            enddo                                  
          enddo
        elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
          ! look through the neighs of last layer's nodes:
          do i=frontstart(lastlayer),frontend(lastlayer)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indi(j)
              if (onfront(node)>0) then
                ghostsz=ghostsz+1
              endif
            enddo                                  
          enddo
        endif
      endif
      ! by now we have counted only incoming (or outgoing) connections on the
      !   interface. And as we do not know yet, whether row or column storage
      !   will be used with A_interf, we take into account also the opposite
      !   connections...
      if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
        ! look through the neighs of first layer's nodes:
        do i=frontstart(1),frontend(1)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indj(j)
            if (onfront(node)==-1) then
              interfsz=interfsz+1
            endif
          enddo                                  
        enddo
      elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
        ! look through the neighs of first layer's nodes:
        do i=frontstart(1),frontend(1)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indi(j)
            if (onfront(node)==-1) then
              interfsz=interfsz+1
            endif
          enddo                                  
        enddo
      endif
      ! build A_interf and A_ghost based on layers...
      A_interf=SpMtx_newInit(interfsz)
      k=0
      if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
        ! at first connections from outside
  !write(stream,*)'frontstart(0)=',frontstart(0)
  !write(stream,*)'frontend(0)=',frontend(0)
  !write(stream,*)'front(',frontstart(0),',',frontend(0),')=',&
  !                front(frontstart(0):frontend(0))
        do i=frontstart(0),frontend(0)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indj(j)
            if (onfront(node)==1) then
              k=k+1
  !write(stream,*)'add to interf:',k,':',front(i),node
              A_interf%indi(k)=front(i)
              A_interf%indj(k)=node
              A_interf%val(k)=A%val(j)
            endif
          enddo                                  
        enddo
        ! connections to outside
  !write(stream,*)'frontstart(1)=',frontstart(1)
  !write(stream,*)'frontend(1)=',frontend(1)
  !write(stream,*)'front(',frontstart(1),',',frontend(1),')=',&
  !                front(frontstart(1):frontend(1))
        do i=frontstart(1),frontend(1)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indj(j)
            if (onfront(node)==-1) then
              k=k+1
  !write(stream,*)'add to interf:',k,':',front(i),node
              A_interf%indi(k)=front(i)
              A_interf%indj(k)=node
              A_interf%val(k)=A%val(j)
            endif
          enddo                                  
        enddo
      elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
        ! at first connections to outside
        do i=frontstart(0),frontend(0)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indi(j)
            if (onfront(node)==1) then
              k=k+1
              A_interf%indi(k)=node
              A_interf%indj(k)=front(i)
              A_interf%val(k)=A%val(j)
            endif
          enddo                                  
        enddo
        ! connections from outside
        do i=frontstart(1),frontend(1)
          do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
            node=A%indi(j)
            if (onfront(node)==-1) then
              k=k+1
              A_interf%indi(k)=node
              A_interf%indj(k)=front(i)
              A_interf%val(k)=A%val(j)
            endif
          enddo                                  
        enddo
      endif
      if (k/=interfsz) then
        write (stream,*)'k=',k,' interfsz=',interfsz,'ol=',ol,' ghostsz=',ghostsz
        write(stream,*)'onfront=',onfront
        call DOUG_abort('SpMtx_build_interfghost -- interfsz wrong!',66)
      endif
      if (ol>0) then
        A_ghost=SpMtx_newInit(ghostsz)
        k=0
        if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
          do i=frontstart(1),frontend(lastlayer)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indj(j)
              if (onfront(node)>0) then
                k=k+1
                A_ghost%indi(k)=front(i)
                A_ghost%indj(k)=node
                A_ghost%val(k)=A%val(j)
              endif
            enddo                                  
          enddo
        elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
          do i=frontstart(1),frontend(lastlayer)
            do j=A%M_bound(front(i)),A%M_bound(front(i)+1)-1
              node=A%indi(j)
              if (onfront(node)>0) then
                k=k+1
                A_ghost%indi(k)=node
                A_ghost%indj(k)=front(i)
                A_ghost%val(k)=A%val(j)
              endif
            enddo                                  
          enddo
        endif
        if (k/=ghostsz) then
          write (stream,*)'k=',k,' ghostsz=',ghostsz
          write(stream,*)'onfront=',onfront
          call DOUG_abort('SpMtx_build_interfghost -- ghostsz wrong!',67)
        endif
      endif
      !now clean up the original matrix from unneeded entries:
      nnz=0
      maxclrnode=0
      if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
        do i=clrstarts(clr),clrstarts(clr+1)-1
          clrnode=clrorder(i)
          if (clrnode>maxclrnode) then
            maxclrnode=clrnode
          endif
          do j=A%M_bound(clrnode),A%M_bound(clrnode+1)-1
            if (onfront(A%indj(j))<=0) then
              nnz=nnz+1
            endif
          enddo
        enddo
      elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
        do i=clrstarts(clr),clrstarts(clr+1)-1
          clrnode=clrorder(i)
          if (clrnode>maxclrnode) then
            maxclrnode=clrnode
          endif
          do j=A%M_bound(clrnode),A%M_bound(clrnode+1)-1
            if (onfront(A%indi(j))<=0) then
              nnz=nnz+1
            endif
          enddo
        enddo
      endif
      A%nnz=nnz
      allocate(itmp(nnz))
      allocate(jtmp(nnz))
      allocate(rtmp(nnz))
      allocate(btmp(maxclrnode+1))
      btmp=0
      nnz=0
      if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
        A%nrows=maxclrnode
        do i=clrstarts(clr),clrstarts(clr+1)-1
          clrnode=clrorder(i)
          do j=A%M_bound(clrnode),A%M_bound(clrnode+1)-1
            if (onfront(A%indj(j))<=0) then
              nnz=nnz+1
              btmp(clrnode+1)=btmp(clrnode+1)+1
              itmp(nnz)=A%indi(j)
              jtmp(nnz)=A%indj(j)
              rtmp(nnz)=A%val(j)
            endif
          enddo
        enddo
      elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
        A%ncols=maxclrnode
        do i=clrstarts(clr),clrstarts(clr+1)-1
          clrnode=clrorder(i)
          do j=A%M_bound(clrnode),A%M_bound(clrnode+1)-1
            if (onfront(A%indi(j))<=0) then
              nnz=nnz+1
              btmp(clrnode+1)=btmp(clrnode+1)+1
              itmp(nnz)=A%indi(j)
              jtmp(nnz)=A%indj(j)
              rtmp(nnz)=A%val(j)
            endif
          enddo
        enddo
      endif
      ! indi
      deallocate(A%indi)
      allocate(A%indi(1:A%nnz))
      A%indi(1:A%nnz)=itmp(1:nnz)
      deallocate(itmp)
      ! indj
      deallocate(A%indj)
      allocate(A%indj(1:A%nnz))
      A%indj(1:A%nnz)=jtmp(1:nnz)
      deallocate(jtmp)
      !M_bound
      btmp(1)=1
      do i=2,maxclrnode+1
        btmp(i)=btmp(i-1)+btmp(i)
      enddo
      deallocate(A%M_bound)
      allocate(A%M_bound(maxclrnode+1))
      A%M_bound=btmp
      deallocate(btmp)
      ! val
      deallocate(A%val)
      allocate(A%val(1:A%nnz))
      A%val(1:A%nnz)=rtmp(1:nnz)
      deallocate(rtmp)
      
      deallocate(frontend)
      deallocate(frontstart)
      deallocate(onfront)
      deallocate(front)
      deallocate(helper)
    end subroutine SpMtx_build_interfghost
  
    subroutine SpMtx_buildAdjncy(A,nedges,xadj,adjncy)
      use globals, only : stream, D_MSGLVL
  