New extract_dep_list from fnd_owner

base/internals/psi_extrct_dl.F90

@@ -126,6 +126,7 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,&
   use psb_const_mod
   use psb_error_mod
   use psb_desc_mod
+  use psb_sort_mod
   implicit none
 #ifdef MPI_H
   include 'mpif.h'
@@ -145,7 +146,8 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,&
   integer(psb_ipk_) :: i,pointer_dep_list,proc,j,err_act
   integer(psb_ipk_) :: err
   integer(psb_ipk_) :: debug_level, debug_unit
-  integer(psb_mpk_) :: iictxt, icomm, me, np,  dl_mpi,  minfo
+  integer(psb_mpk_) :: iictxt, icomm, me, np, minfo
+  logical, parameter :: dist_symm_list=.false., print_dl=.false.
   character  name*20
   name='psi_extrct_dl'
 
@@ -156,13 +158,14 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,&
   info = psb_success_
 
   call psb_info(iictxt,me,np)
-  allocate(itmp(np+1),length_dl(0:np),stat=info)
+  allocate(itmp(2*np+1),length_dl(0:np),stat=info)
   if (info /= psb_success_) then 
     info=psb_err_alloc_dealloc_
     goto 9999
   end if
   do i=0,np 
     length_dl(i) = 0
+    itmp(i+1)    = -1
   enddo
   i=1
   if (debug_level >= psb_debug_inner_)&
@@ -249,16 +252,13 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,&
   endif
 
   length_dl(me)=pointer_dep_list-1
-  dl_lda = max(length_dl(me),1)
+  if (dist_symm_list) then 
+    call psb_realloc(length_dl(me),itmp,info)
+    call psi_symm_dep_list(itmp,ictxt,info)       
+    dl_lda = max(size(itmp),1)
     call psb_max(iictxt, dl_lda)
-  !
+    call psb_realloc(dl_lda,itmp,info)
-  ! This doubling of DL_LDA is not 100% safe,
+    !  dl_lda = min(np,2*dl_lda)
-  ! but should work most of the time.
-  ! Will need to be improved later, perhaps move
-  ! from a 2D allocation (ellpack style) to
-  ! a 1D allocation (csr like). 
-  !
-  dl_lda = min(2*dl_lda,np+1)
     allocate(dep_list(dl_lda,0:np),stat=info)
     if (info /= psb_success_) then 
       call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
@@ -275,6 +275,93 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,&
       info=psb_err_alloc_dealloc_
       goto 9999
     endif
+  else
+    block
+      integer(psb_ipk_), allocatable :: list1(:,:), ldl2(:), list2(:,:)
+      integer(psb_ipk_) :: i,j,ip,dlsym, ldu, mdl, l1, l2
+
+      dl_lda = max(length_dl(me),1)
+      call psb_max(iictxt, dl_lda)
+      allocate(dep_list(dl_lda,0:np),stat=info)
+      if (info /= psb_success_) then 
+        call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
+        goto 9999      
+      end if
+      call psb_sum(iictxt,length_dl(0:np))
+      icomm = psb_get_mpicomm(iictxt)
+      call mpi_allgather(itmp,dl_lda,psb_mpi_ipk_,&
+           & dep_list,dl_lda,psb_mpi_ipk_,icomm,minfo)
+      info = minfo  
+      if (info /= psb_success_) then 
+        info=psb_err_alloc_dealloc_
+        goto 9999
+      endif
+      allocate(ldl2(0:np),stat=info)
+      ldl2 = 0 
+      do j=0, np-1
+        do i=1,length_dl(j)
+          ip = dep_list(i,j)
+          ldl2(ip) = ldl2(ip) + 1
+        end do
+      end do
+      dlsym = maxval(ldl2)
+      allocate(list2(dlsym,0:np),stat=info)
+      call move_alloc(dep_list,list1)
+      if (info /= psb_success_) then 
+        call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
+        goto 9999      
+      end if
+      ldl2 = 0 
+      do j=0, np-1
+        do i=1,length_dl(j)
+          ip       = list1(i,j)
+          ldl2(ip) = ldl2(ip) + 1
+          list2(ldl2(ip),ip) = j
+        end do
+      end do
+      mdl = 0 
+      do j = 0, np-1
+        l1 = length_dl(j)
+        l2 = ldl2(j) 
+        itmp(1:l1)       = list1(1:l1,j)
+        itmp(l1+1:l1+l2) = list2(1:l2,j)
+        ldu = l1 + l2
+        !if (me == 0) write(0,*) 'Iter ',j,':',l1,l2,':',itmp(1:l1),':',itmp(l1+1:l1+l2)
+        call psb_msort_unique(itmp(1:l1+l2),ldu)
+        mdl = max(mdl, ldu)
+        !if (me == 0) write(0,*) 'Iter ',j,':',ldu,':',itmp(1:ldu)
+      end do
+      dl_lda = mdl
+      allocate(dep_list(dl_lda,0:np),stat=info)
+      dep_list = -1
+      if (info /= psb_success_) then 
+        call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
+        goto 9999      
+      end if
+      do j = 0, np-1
+        l1 = length_dl(j)
+        l2 = ldl2(j) 
+        itmp(1:l1)       = list1(1:l1,j)
+        itmp(l1+1:l1+l2) = list2(1:l2,j)
+        ldu = l1 + l2 
+        call psb_msort_unique(itmp(1:l1+l2),ldu)
+        length_dl(j) = ldu
+        dep_list(1:ldu,j) = itmp(1:ldu)
+      end do
+      
+    end block
+  endif
+  if (print_dl) then
+    if (me == 0) then
+      write(0,*) ' Dep_list '
+      do i=0,np-1
+        j = length_dl(i) 
+        write(0,*) 'Proc ',i,':',dep_list(1:j,i)
+      end do
+      flush(0)
+    end if
+    call psb_barrier(ictxt)
+  end if
 
   call psb_erractionrestore(err_act)
   return

base/internals/psi_symm_dep_list.F90

@@ -37,13 +37,13 @@
 ! 
 ! Arguments: 
 ! 
-subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
+subroutine psi_symm_dep_list_inrv(rvsz,adj,ictxt,info)
   use psb_serial_mod
   use psb_const_mod
   use psb_error_mod
   use psb_penv_mod
   use psb_realloc_mod
-  use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list
+  use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list_inrv
 #ifdef MPI_MOD
   use mpi
 #endif
@@ -54,30 +54,23 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
 #endif
   integer(psb_mpk_), intent(inout)   :: rvsz(0:)
   integer(psb_ipk_), allocatable, intent(inout) :: adj(:)
-  class(psb_indx_map), intent(in) :: idxmap
+  integer(psb_ipk_), intent(in)      :: ictxt
   integer(psb_ipk_), intent(out)     :: info
-  integer(psb_ipk_), intent(in), optional :: flag
   
   !
   integer(psb_ipk_), allocatable :: ladj(:) 
-  integer(psb_mpk_), allocatable :: sdsz(:)
+  integer(psb_mpk_) :: icomm, minfo
-  integer(psb_mpk_) :: icomm, minfo, iictxt
   integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,&
        & last_ih, last_j, nidx, nrecv, nadj, flag_
   integer(psb_lpk_) :: mglob, ih
-  integer(psb_ipk_) :: ictxt,np,me
+  integer(psb_ipk_) :: np,me
   character(len=20)   :: name
 
   info = psb_success_
   name = 'psi_symm_dep_list'
   call psb_erractionsave(err_act)
 
-  ictxt   = idxmap%get_ctxt()
+  icomm = psb_get_mpicomm(ictxt)
-  icomm   = idxmap%get_mpic()
-  mglob   = idxmap%get_gr()
-  n_row   = idxmap%get_lr()
-  n_col   = idxmap%get_lc()
-  iictxt = ictxt 
 
   call psb_info(ictxt, me, np)
 
@@ -87,38 +80,92 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
     goto 9999
   endif
 
-  if (.not.(idxmap%is_valid())) then 
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap')
-    goto 9999      
-  end if
-
-
   nadj = size(adj)
+
+  !
+  ! Am getting this from the outside, possibly from the wrapper call
+  !
   if (size(rvsz)<np) then 
     call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid rvsz')
     goto 9999      
   end if
+  nrecv = 0 
+  do ip=0, np-1
+    if (rvsz(ip) > 0) nrecv = nrecv + 1
+  end do
 
-  if (present(flag)) then
+  !
-    flag_ = flag
+  !  Now fix adj to be symmetric
-  else
+  !
-    flag_ = psi_symm_flag_norv_
+  call psb_realloc(nadj+nrecv,ladj,info)
+  ladj(1:nadj) = adj(1:nadj)
+  do ip=0, np-1
+    if (rvsz(ip)>0) then
+      nadj  = nadj + 1 
+      ladj(nadj) = ip
     end if
+  end do
+  call psb_msort_unique(ladj,nadj)
+  call psb_realloc(nadj,adj,info)
+  if (info == 0) adj(1:nadj) = ladj(1:nadj) 
   
-  select case(flag_)
+  call psb_erractionrestore(err_act)
-  case(psi_symm_flag_norv_, psi_symm_flag_inrv_)
+  return
-    ! Ok, keep going
-  case default
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid flag')
-    goto 9999      
-  end select
 
+9999 call psb_error_handler(ictxt,err_act)
 
-  if (flag_ == psi_symm_flag_norv_) then 
+  return
+
+end subroutine psi_symm_dep_list_inrv
+
+subroutine psi_symm_dep_list_norv(adj,ictxt,info)
+  use psb_serial_mod
+  use psb_const_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_realloc_mod
+  use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list_norv
+#ifdef MPI_MOD
+  use mpi
+#endif
+
+  implicit none 
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  integer(psb_ipk_), allocatable, intent(inout) :: adj(:)
+  integer(psb_ipk_), intent(in)      :: ictxt
+  integer(psb_ipk_), intent(out)     :: info
+  
+  !
+  integer(psb_mpk_), allocatable :: rvsz(:), sdsz(:) 
+
+  integer(psb_mpk_) :: icomm, minfo
+  integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,&
+       & last_ih, last_j, nidx, nrecv, nadj
+  integer(psb_ipk_) :: mglob, ih
+  integer(psb_ipk_) :: np,me
+  character(len=20) :: name
+
+  info = psb_success_
+  name = 'psi_symm_dep_list'
+  call psb_erractionsave(err_act)
+
+  icomm = psb_get_mpicomm(ictxt)
+
+  call psb_info(ictxt, me, np)
+
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  nadj = size(adj)
   
   ! write(0,*) me,name,' Going through ',nidx,nadj
   
-    Allocate(sdsz(0:np-1), stat=info)
+  Allocate(sdsz(0:np-1), rvsz(0:np-1), stat=info)
   !
   ! First, send sizes according to adjcncy list
   !
@@ -130,28 +177,12 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
   
   call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
        & rvsz,1,psb_mpi_mpk_,icomm,minfo)
+  if (minfo == 0)  call psi_symm_dep_list(rvsz,adj,ictxt,info)
+  if ((minfo /=0).or.(info /= 0)) then
+    call psb_errpush(psb_err_from_subroutine_,name,a_err='inner call symm_dep')
+    goto 9999      
   end if
 
-  nrecv = 0 
-  do ip=0, np-1
-    if (rvsz(ip) > 0) nrecv = nrecv + 1
-  end do
-
-  !
-  !  Now fix adj to be symmetric
-  !
-  call psb_realloc(nadj+nrecv,ladj,info)
-  ladj(1:nadj) = adj(1:nadj)
-  do ip=0, np-1
-    if (rvsz(ip)>0) then
-      nadj  = nadj + 1 
-      ladj(nadj+1:nadj+nrecv) = ip
-    end if
-  end do
-  call psb_msort_unique(ladj,nadj)
-  call psb_realloc(nadj,adj,info)
-  adj(1:nadj) = ladj(1:nadj) 
-  
   call psb_erractionrestore(err_act)
   return
 
@@ -159,4 +190,4 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
 
   return
 
-end subroutine psi_symm_dep_list
+end subroutine psi_symm_dep_list_norv

base/modules/desc/psb_indx_map_mod.f90

@@ -313,17 +313,23 @@ module psb_indx_map_mod
 
   integer, parameter :: psi_symm_flag_norv_ = 0
   integer, parameter :: psi_symm_flag_inrv_ = 1
-  interface 
+  interface psi_symm_dep_list
-    subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag)
+    subroutine psi_symm_dep_list_inrv(rvsz,adj,ictxt,info)
       import :: psb_indx_map, psb_ipk_, psb_lpk_, psb_mpk_
       implicit none 
-      integer(psb_mpk_), intent(inout) :: rvsz(:)
+      integer(psb_mpk_), intent(inout)   :: rvsz(0:)
-      integer(psb_ipk_), intent(in)    :: adj(:)
+      integer(psb_ipk_), allocatable, intent(inout) :: adj(:)
-      class(psb_indx_map), intent(in)  :: idxmap
+      integer(psb_ipk_), intent(in)      :: ictxt
       integer(psb_ipk_), intent(out)     :: info
-      integer(psb_ipk_), intent(in), optional :: flag
+    end subroutine psi_symm_dep_list_inrv
-    end subroutine psi_symm_dep_list
+    subroutine psi_symm_dep_list_norv(adj,ictxt,info)
-  end interface
+      import :: psb_indx_map, psb_ipk_, psb_lpk_, psb_mpk_
+      implicit none 
+      integer(psb_ipk_), allocatable, intent(inout) :: adj(:)
+      integer(psb_ipk_), intent(in)      :: ictxt
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_symm_dep_list_norv    
+  end interface psi_symm_dep_list
 
 contains