anchor.diff:            if (numprocs==1) then 
               A_interf_%nnz=0
             endif
  -write(stream,*)'################### calling:multi_subsolve'
  -call flush(stream)
             call multi_subsolve(               &
                    nids=A%nsubsolves,          &
                    ids=A%subsolve_ids,         &
  
positive.code:                  val_interf=A_interf_%val)
          else
            if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
  write(stream,*)'##########B######## calling:multi_subsolve'
  call flush(stream)
              call multi_subsolve(               &
                     nids=A%nsubsolves,          &
                     ids=A%subsolve_ids,         &
  
negative.code:     if (associated(subsol)) deallocate(subsol)
      if (associated(subrhs)) deallocate(subrhs)
    end subroutine multi_subsolve
  
    !subroutine multiplicative_multi_subsolve(A,M,nids,ids,sol,rhs,res,subd,nfreds,nnz,indi,indj,val, &
    subroutine multiplicative_multi_subsolve(A,M,ids,sol,rhs,res,subd,nfreds,nnz,indi,indj,val, &
      nnz_interf,indi_interf,indj_interf,val_interf, &
      nagr1,starts1,nodes1, & ! fine level aggregate list
      nagr2,starts2,nodes2, & ! coarse level aggregate list
      starts3,nodes3)         ! restrict op. starts & nodelist
      use SpMtx_mods
      use Mesh_class
  !chk use SpMtx_util
      implicit none
      type(SpMtx),intent(inout) :: A
      type(Mesh),intent(in) :: M ! Mesh
      !integer,intent(inout) :: nids
      integer :: nids
      integer,dimension(:),pointer :: ids
      real(kind=rk),dimension(:),pointer :: sol,rhs,res
      type(indlist),dimension(:),pointer :: subd
      integer,intent(in),optional :: nfreds
      integer,intent(in),optional :: nnz
      integer,dimension(:),intent(in),optional :: indi,indj
      real(kind=rk),dimension(:),optional :: val
      integer,intent(in),optional :: nnz_interf
      integer,dimension(:),intent(in),optional :: indi_interf,indj_interf
      real(kind=rk),dimension(:),optional :: val_interf
      integer,intent(in),optional :: nagr1,nagr2
      integer,dimension(:),pointer,optional :: &
                                          starts1,starts2,starts3
      integer,dimension(:),pointer,optional :: nodes1,nodes2,nodes3
      ! ----- local: ------
      integer :: i,ii,j,n,nod1,nod2,sd,agr1,agr2,snnz,nnz_per_fred
      integer,save :: nnz_est=0
      integer :: nselind ! number of selected indeces
      integer,dimension(:),allocatable :: floc ! freedom location in selind array
      integer,dimension(:),pointer :: selind ! selected indeces
      integer,dimension(:),allocatable :: sindi,sindj ! selected indeces
      real(kind=rk),dimension(:),allocatable :: sval ! selected values
      real(kind=rk),dimension(:),pointer :: subrhs,subsol
      real(kind=rk) :: t1
      logical :: dofactorise
      logical,save :: upward=.true.,full
  
  
      !chk real(kind=rk),dimension(:),pointer,save :: vmp,vsval
      !chk integer,dimension(:),allocatable,save :: vsindi,vsindj ! selected indeces
      
      call SpMtx_arrange(A,arrange_type=D_SpMtx_ARRNG_ROWS,sort=.true.)   !
      if (present(nagr1).or.present(nagr2).or.&
          present(val).or.present(nnz_interf)) then
        dofactorise=.true.
      else
        dofactorise=.false.
      endif
      if (upward) then
        if (dofactorise) then
          !-----------------------------------------
          res=rhs
          sol=0.0_rk
          !-----------------------------------------
          if (present(nagr2)) then ! Need still to collect together fine aggregates
            nids=nagr2
            A%nsubsolves=nids
            if (.not.present(nagr1)) then
              call DOUG_abort('[multi_subsolve] agr1 must be present!',-1)
            endif
          elseif (present(nagr1)) then ! form subdomains based on fine aggregates
            nids=nagr1
            A%nsubsolves=nids
          elseif (present(nnz_interf)) then ! form subdomains based on fine aggregates
            nids=1
            A%nsubsolves=nids
          endif
          if (.not.associated(subrhs)) then
            allocate(subrhs(nfreds)) ! todo: could be somewhat economised...
          endif
          if (.not.associated(subsol)) then
            allocate(subsol(nfreds))
          endif
          setuptime=0.0_rk
          factorisation_time=0.0_rk
          t1 = MPI_WTIME()
          nselind=0
          if (.not.present(nfreds)) then
            call DOUG_abort('[multi_subsolve] nfreds must be present!',-1)
          endif
          allocate(floc(nfreds))
          allocate(selind(nfreds))
          if (present(nnz_interf)) then
            allocate(sindi(nnz+nnz_interf)) ! todo: this is in most cases far too much
            allocate(sindj(nnz+nnz_interf))
            allocate(sval(nnz+nnz_interf))
          else
            allocate(sindi(nnz)) ! todo: this is in most cases far too much
            allocate(sindj(nnz))
            allocate(sval(nnz))
          endif
          if (present(nagr2)) then !{
            if (present(starts3)) then !{
              do agr2=1,nagr2 ! loop over Restriction TODO here!
                floc=0
                nselind=0
                do agr1=starts2(agr2),starts2(agr2+1)-1 ! look through fine agrs.
                  nod2=nodes2(agr1) ! actual aggregate numbers
                  do i=starts3(nod2),starts3(nod2+1)-1 ! loop over nodes in aggr.
                    nod1=nodes3(i) ! node number
                    if (floc(nod1)==0) then
                      nselind=nselind+1
                      floc(nod1)=nselind
                      selind(nselind)=nod1
                    endif
                  enddo
                enddo
                ! now we have gathered information for subd agr2...
                snnz=0
                do i=1,nnz
                  ! todo:
                  !   need to be avoided going it all through again and again?
                  ! idea: to use linked-list arrays as in old DOUG
                  if (floc(indi(i))>0) then ! wheather in the subdomain?
                    if (floc(indj(i))>0) then
                      snnz=snnz+1
                      sindi(snnz)=floc(indi(i))
                      sindj(snnz)=floc(indj(i))
                      sval(snnz)=val(i)
                    elseif (sctls%overlap>1) then ! connection to the overlap
                      snnz=snnz+1
                      sindi(snnz)=floc(indi(i))
                      if (floc(indj(i))==0) then ! at the overlap node the 1st time
                        nselind=nselind+1
                        selind(nselind)=indj(i)
                        sindj(snnz)=nselind
                        floc(indj(i))=-nselind
                      else ! node already added to the overlap
                        sindj(snnz)=-floc(indj(i))
                      endif
                      sval(snnz)=val(i)
                    !elseif (sctls%overlap==1) then
                    endif
                  endif
                enddo
                subrhs(1:nselind)=res(selind(1:nselind))
                if (nnz_est==0) then
                  nnz_per_fred=snnz/nselind+1
                  nnz_est=nfreds*nnz_per_fred
                endif
                if (agr2<10.or.agr2>nagr2-10) then
                  write(stream,*)'Subd.',agr2,' size,nnz:',nselind,snnz,&
                    ' #fine_aggrs:',starts2(agr2+1)-starts2(agr2)
                endif
                setuptime=setuptime+(MPI_WTIME()-t1) ! switchoff clock
        
                call factorise_and_solve(ids(agr2),subsol,subrhs,nselind,snnz,sindi,sindj,sval)
        
                t1 = MPI_WTIME() ! switch on clock
                sol(selind(1:nselind))=sol(selind(1:nselind))+subsol(1:nselind)
                !-----------------------------------------
                ! now update also the residual: !TODO -- compute this only on the
                !                                        affected freedoms.
                if (agr2==1) then
                  full=.true.
                else
                  full=.false.
                endif
                call find_sub_residual(res,A,sol,rhs,nselind,selind,full)
                !-----------------------------------------
                allocate(subd(agr2)%inds(nselind))
                subd(agr2)%ninds=nselind
                subd(agr2)%inds(1:nselind)=selind(1:nselind)
              enddo
            else !}{
              do agr2=1,nagr2 ! loop over coarse aggregates
                floc=0
                nselind=0
                do agr1=starts2(agr2),starts2(agr2+1)-1 ! look through fine agrs.
                  nod2=nodes2(agr1) ! actual aggregate numbers
                  do i=starts1(nod2),starts1(nod2+1)-1 ! loop over nodes in aggr.
                    nod1=nodes1(i) ! node number
                    if (floc(nod1)==0) then
                      nselind=nselind+1
                      floc(nod1)=nselind
                      selind(nselind)=nod1
                    endif
                  enddo
                enddo
                ! now we have gathered information for subd agr2...
                snnz=0
                do i=1,nnz
                  ! todo:
                  !   need to be avoided going it all through again and again?
                  ! idea: to use linked-list arrays as in old DOUG
                  if (floc(indi(i))>0) then ! wheather in the subdomain?
                    if (floc(indj(i))>0) then
                      snnz=snnz+1
                      sindi(snnz)=floc(indi(i))
                      sindj(snnz)=floc(indj(i))
                      sval(snnz)=val(i)
                    elseif (sctls%overlap>1) then ! connection to the overlap
                      snnz=snnz+1
                      sindi(snnz)=floc(indi(i))
                      if (floc(indj(i))==0) then ! at the overlap node the 1st time
                        nselind=nselind+1
                        selind(nselind)=indj(i)
                        sindj(snnz)=nselind
                        floc(indj(i))=-nselind
                      else ! node already added to the overlap
                        sindj(snnz)=-floc(indj(i))
                      endif
                      sval(snnz)=val(i)
                    !elseif (sctls%overlap==1) then
                    endif
                  endif
                enddo
                if (sctls%overlap>1) then ! connections from nodes on overlap
                  do i=1,nnz
                    if (floc(indi(i))<0) then ! from overlap
                      if(floc(indj(i))>0) then ! to inner
                        snnz=snnz+1
                        sindi(snnz)=-floc(indi(i))
                        sindj(snnz)=floc(indj(i))
                        sval(snnz)=val(i)
                      elseif(floc(indj(i))<0) then ! to overlap
                        snnz=snnz+1
                        sindi(snnz)=-floc(indi(i))
                        sindj(snnz)=-floc(indj(i))
                        sval(snnz)=val(i)
                      endif
                    endif
                  enddo
                endif
                subrhs(1:nselind)=res(selind(1:nselind))
                if (nnz_est==0) then
                  nnz_per_fred=snnz/nselind+1
                  nnz_est=nfreds*nnz_per_fred
                endif
                if (agr2<10.or.agr2>nagr2-10) then
                  write(stream,*)'Subd.',agr2,' size,nnz:',nselind,snnz,&
                    ' #fine_aggrs:',starts2(agr2+1)-starts2(agr2)
                endif
                setuptime=setuptime+(MPI_WTIME()-t1) ! switchoff clock
        
                call factorise_and_solve(ids(agr2),subsol,subrhs,nselind,snnz,sindi,sindj,sval)
        
                t1 = MPI_WTIME() ! switch on clock
                sol(selind(1:nselind))=sol(selind(1:nselind))+subsol(1:nselind)
                !-----------------------------------------
                ! now update also the residual: !TODO -- compute this only on the
                !                                        affected freedoms.
                if (agr2==1) then
                  full=.true.
                else
                  full=.false.
                endif
                call find_sub_residual(res,A,sol,rhs,nselind,selind,full)
                !-----------------------------------------
                ! keep indlist:
                allocate(subd(agr2)%inds(nselind))
                subd(agr2)%ninds=nselind
                subd(agr2)%inds(1:nselind)=selind(1:nselind)
              enddo
            endif !}
          elseif (present(nagr1)) then !}{
            do agr1=1,nagr1 ! look through fine agrs.
              floc=0
              nselind=0
              do i=starts1(agr1),starts1(agr1+1)-1 ! loop over nodes in aggr.
                nod1=nodes1(i) ! node number
                if (floc(nod1)==0) then
                  nselind=nselind+1
                  floc(nod1)=nselind
                  selind(nselind)=nod1
                endif
              enddo
              ! now we have gathered information for subd agr1...
              snnz=0
              do i=1,nnz
                ! todo:
                !   need to be avoided going it all through again and again?
                ! idea: to use linked-list arrays as in old DOUG
                if (floc(indi(i))>0) then ! wheather in the subdomain?
                  if (floc(indj(i))>0) then
                    snnz=snnz+1
                    sindi(snnz)=floc(indi(i))
                    sindj(snnz)=floc(indj(i))
                    sval(snnz)=val(i)
                  elseif (sctls%overlap>1) then ! connection to the overlap
                    snnz=snnz+1
                    sindi(snnz)=floc(indi(i))
                    if (floc(indj(i))==0) then ! at the overlap node the 1st time
                      nselind=nselind+1
                      selind(nselind)=indj(i)
                      sindj(snnz)=nselind
                      floc(indj(i))=-nselind
                    else ! node already added to the overlap
                      sindj(snnz)=-floc(indj(i))
                    endif
                    sval(snnz)=val(i)
                  !elseif (sctls%overlap==1) then
                  endif
                endif
              enddo
              if (sctls%overlap>1) then ! connections from nodes on overlap
                do i=1,nnz
                  if (floc(indi(i))<0) then ! from overlap
                    if(floc(indj(i))>0) then ! to inner
                      snnz=snnz+1
                      sindi(snnz)=-floc(indi(i))
                      sindj(snnz)=floc(indj(i))
                      sval(snnz)=val(i)
                    elseif(floc(indj(i))<0) then ! to overlap
                      snnz=snnz+1
                      sindi(snnz)=-floc(indi(i))
                      sindj(snnz)=-floc(indj(i))
                      sval(snnz)=val(i)
                    endif
                  endif
                enddo
              endif
              subrhs(1:nselind)=res(selind(1:nselind))
              if (nnz_est==0) then
                nnz_per_fred=snnz/nselind+1
                nnz_est=nfreds*nnz_per_fred
              endif
              if (agr1<20.or.agr1>nagr1-20) then
                write(stream,*)'Subdomain ',agr1,' size:',nselind
              endif
              setuptime=setuptime+(MPI_WTIME()-t1) ! switchoff clock
        
              call factorise_and_solve(ids(agr1),subsol,subrhs,nselind,snnz,sindi,sindj,sval)
        
              t1 = MPI_WTIME() ! switch on clock
              sol(selind(1:nselind))=sol(selind(1:nselind))+subsol(1:nselind)
              !-----------------------------------------
              ! now update also the residual: !TODO -- compute this only on the
              !                                        affected freedoms.
              if (agr1==1) then
                full=.true.
              else
                full=.false.
              endif
              call find_sub_residual(res,A,sol,rhs,nselind,selind,full)
              !-----------------------------------------
              ! keep indlist:
              allocate(subd(agr1)%inds(nselind)) !todo: siin mingi jama?
              subd(agr1)%ninds=nselind
              subd(agr1)%inds(1:nselind)=selind(1:nselind)
            enddo
          else !}{ ! this case is not considered with multiplicative yet now...
            call DOUG_Abort('1-subdomain-per-proc not considered yet',5464)
          endif !}
          deallocate(sval)
          deallocate(sindj)
          deallocate(sindi)
          deallocate(selind)
          deallocate(floc)
          setuptime=setuptime+(MPI_WTIME()-t1) ! switchoff clock
        else ! perform backsolves:
          res=rhs
          nids=A%nsubsolves
          sol=0.0_rk
          if (.not.associated(subrhs)) then
            allocate(subrhs(maxval(subd(1:nids)%ninds)))
          endif
          if (.not.associated(subsol)) then
            allocate(subsol(maxval(subd(1:nids)%ninds)))
          endif
          do sd=1,nids
            n=subd(sd)%ninds
            subrhs(1:n)=res(subd(sd)%inds(1:n))
            call factorise_and_solve(ids(sd),subsol,subrhs,n)
        
            sol(subd(sd)%inds(1:n))=sol(subd(sd)%inds(1:n))+subsol(1:n)
            !-----------------------------------------
            ! now update also the residual: !TODO -- compute this only on the
            !                                        affected freedoms.
            if (sd==1) then
              full=.true.
            else
              full=.false.
            endif
            call find_sub_residual(res,A,sol,rhs,n,subd(sd)%inds,full)
            !-----------------------------------------
          enddo
        endif
        upward=.false. ! next time downward
      else
        nids=A%nsubsolves
        if (.not.associated(subrhs)) then
          allocate(subrhs(maxval(subd(1:nids)%ninds)))
        endif
        if (.not.associated(subsol)) then
          allocate(subsol(maxval(subd(1:nids)%ninds)))
        endif
        nids=A%nsubsolves
        ! now the other-way-round for symmetry:
        do sd=nids,1,-1
          n=subd(sd)%ninds
          subrhs(1:n)=res(subd(sd)%inds(1:n))
        
          call factorise_and_solve(ids(sd),subsol,subrhs,n)
        
          sol(subd(sd)%inds(1:n))=sol(subd(sd)%inds(1:n))+subsol(1:n)
          if (sd>1) then
            !-----------------------------------------
            ! now update also the residual: !TODO -- compute this only on the
            !                                        affected freedoms.
            call find_sub_residual(res,A,sol,rhs,n,subd(sd)%inds,full)
            !-----------------------------------------
          endif
        enddo
        upward=.true. ! next time upward
      endif
      if (associated(subsol)) deallocate(subsol)
      if (associated(subrhs)) deallocate(subrhs)
    end subroutine multiplicative_multi_subsolve
  
    ! the very first, far from optimised version:
    subroutine find_sub_residual(res,A,sol,rhs,nselind,selind,full)
      use SpMtx_mods
      implicit none
      type(SpMtx) :: A
      real(kind=rk),dimension(:),pointer :: res,sol,rhs
      integer :: nselind
      integer,dimension(:),pointer :: selind
      logical :: full
      integer,dimension(:),allocatable :: isin
      integer :: i,ii,jj
      if (full) then
        res=rhs
        do i=1,A%nnz
          ii=A%indi(i)
          jj=A%indj(i)
          res(ii)=res(ii)-A%val(i)*sol(jj)
        enddo
      else
        allocate(isin(size(rhs)))
        isin=0
        isin(selind(1:nselind))=1
        do i=1,A%nnz
          ii=A%indi(i)
          jj=A%indj(i)
          if (isin(ii)==1.or.isin(jj)==1) then
            if (isin(ii)==0) then
              isin(ii)=2
            endif
            if (isin(jj)==0) then
              isin(jj)=2
            endif
          endif
        enddo
        do i=1,A%nnz
          ii=A%indi(i)
          jj=A%indj(i)
          if (isin(ii)>0) then
            res(ii)=rhs(ii)
            isin(ii)=-1
          endif
          !if (isin(ii)/=0.or.isin(jj)/=0) then
          if (isin(ii)/=0) then
            res(ii)=res(ii)-A%val(i)*sol(jj)
          endif
        enddo
        deallocate(isin)
      endif
  
             !isin=.false.
             !isin(selind(1:nselind))=.true.
             !do i=1,nselind
             !  ii=selind(i)
             !  !!rhs(ii)=0.0_rk
             !  !!!rhs(ii)=0.0_rk
             !  !!!!rhs(ii)=-rhs_in(ii)
             !  rhs(ii)=-b(ii)
             !  do j=A%M_bound(ii),A%M_bound(ii+1)-1
             !    if (isin(A%indj(j))) then
             !      rhs(ii)=rhs(ii)+A%val(j)*sol(A%indj(j))
             !      !!!rhs(ii)=-rhs(ii)+A%val(j)*sol(A%indj(j))
             !    endif
             !  enddo
             !enddo
    end subroutine find_sub_residual
  
  
    subroutine sparse_singlesolve(id,sol,rhs,nfreds,nnz,indi,indj,val,tot_nfreds,nnz_est)
      ! For adding a new factorised matrix id must be 0.
  