psblas-dev

base/modules/psb_inter_desc_type.f90
 base/modules/psb_tools_mod.f90
 base/tools/psb_zspins.f90

Defined functions for general linear mappings among index spaces:
  psb_linmap_init
  psb_linmap_ins
  psb_linmap_asb

base/modules/psb_inter_desc_type.f90

@@ -36,8 +36,8 @@
 !    to different spaces.
 !
 module psb_inter_descriptor_type
-  use psb_spmat_type
+  use psb_spmat_type, only : psb_dspmat_type, psb_zspmat_type
-  use psb_descriptor_type
+  use psb_descriptor_type, only : psb_desc_type
   
   
 
@@ -66,7 +66,7 @@ module psb_inter_descriptor_type
     integer, allocatable :: itd_data(:)
     type(psb_desc_type), pointer :: desc_1=>null(), desc_2=>null()
     integer, allocatable :: exch_fw_idx(:), exch_bk_idx(:)
-    type(psb_desc_type)  :: desc_ext_1, desc_ext_2
+    type(psb_desc_type)  :: desc_fw, desc_bk
     type(psb_d_map_type) :: dmap
     type(psb_z_map_type) :: zmap
   end type psb_inter_desc_type
@@ -97,6 +97,9 @@ module psb_inter_descriptor_type
          & psb_d_map_sizeof, psb_z_map_sizeof
   end interface
 
+  interface psb_linmap
+    module procedure  psb_d_apply_linmap, psb_z_apply_linmap
+  end interface
 
 contains
  
@@ -177,6 +180,8 @@ contains
 
 
   logical function psb_is_asb_inter_desc(desc)
+    use psb_descriptor_type
+    implicit none 
     type(psb_inter_desc_type), intent(in) :: desc
 
     psb_is_asb_inter_desc = .false.
@@ -189,6 +194,8 @@ contains
   end function psb_is_asb_inter_desc
 
   logical function psb_is_ok_inter_desc(desc)
+    use psb_descriptor_type
+    implicit none 
     type(psb_inter_desc_type), intent(in) :: desc
 
     psb_is_ok_inter_desc = .false.
@@ -209,6 +216,7 @@ contains
 
 
   function psb_d_map_sizeof(map)
+    use psb_spmat_type
     implicit none
     type(psb_d_map_type), intent(in) :: map
     Integer                      :: psb_d_map_sizeof
@@ -222,6 +230,7 @@ contains
   end function psb_d_map_sizeof
 
   function psb_z_map_sizeof(map)
+    use psb_spmat_type
     implicit none
     type(psb_z_map_type), intent(in) :: map
     Integer                      :: psb_z_map_sizeof
@@ -235,8 +244,9 @@ contains
   end function psb_z_map_sizeof
 
   function psb_itd_sizeof(desc)
-    
+    use psb_spmat_type
-    implicit none
+    use psb_descriptor_type
+    implicit none 
     type(psb_inter_desc_type), intent(in) :: desc
     Integer                      :: psb_itd_sizeof
     integer :: val
@@ -246,14 +256,16 @@ contains
     if (allocated(desc%itd_data))    val = val + 4*size(desc%itd_data)
     if (allocated(desc%exch_fw_idx)) val = val + 4*size(desc%exch_fw_idx)
     if (allocated(desc%exch_bk_idx)) val = val + 4*size(desc%exch_bk_idx)
-    val = val + psb_sizeof(desc%desc_ext_1)
+    val = val + psb_sizeof(desc%desc_fw)
-    val = val + psb_sizeof(desc%desc_ext_2)
+    val = val + psb_sizeof(desc%desc_bk)
     val = val + psb_sizeof(desc%dmap)
     val = val + psb_sizeof(desc%zmap)
     psb_itd_sizeof = val 
   end function psb_itd_sizeof
   
   function psb_d_inter_desc(map_kind,desc1,desc2,map_fw,map_bk,idx_fw,idx_bk)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_serial_mod
     use psi_mod
     implicit none 
@@ -299,6 +311,8 @@ contains
   end function psb_d_inter_desc
   
   function psb_d_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_serial_mod
     use psi_mod
     implicit none 
@@ -350,6 +364,8 @@ contains
   end function psb_d_inter_desc_noidx
 
   function psb_z_inter_desc(map_kind,desc1, desc2, map_fw, map_bk, idx_fw, idx_bk)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_serial_mod
     use psi_mod
     implicit none 
@@ -395,6 +411,8 @@ contains
   end function psb_z_inter_desc
   
   function psb_z_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_serial_mod
     use psi_mod
     implicit none 
@@ -454,6 +472,8 @@ contains
   ! due to exch_fw_idx
   !
   subroutine psb_d_forward_map(alpha,x,beta,y,desc,info,work)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_comm_mod
     use psb_serial_mod
     use psi_mod
@@ -501,8 +521,9 @@ contains
       end if
 
     case(psb_map_gen_linear_)
-      call psb_halo(x,desc%desc_ext_1,info,work=work)
+
-      if (info == 0) call psb_csmm(alpha,desc%dmap%map_fw,x,beta,y,info)
+      call psb_linmap(alpha,x,beta,y,desc%dmap%map_fw,&
+           & desc%desc_fw,desc%desc_1,desc%desc_2)
         
       if (info /= 0) then 
         write(0,*) trim(name),' Error from inner routines',info
@@ -525,6 +546,8 @@ contains
   ! due to exch_bk_idx
   !
   subroutine psb_d_backward_map(alpha,x,beta,y,desc,info,work)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_comm_mod
     use psb_serial_mod
     use psi_mod
@@ -572,9 +595,8 @@ contains
       
 
     case(psb_map_gen_linear_)
-      call psb_halo(x,desc%desc_ext_2,info,work=work)
+      call psb_linmap(alpha,x,beta,y,desc%dmap%map_bk,&
-      if (info == 0) call psb_csmm(alpha,desc%dmap%map_bk,x,beta,y,info)
+           & desc%desc_bk,desc%desc_2,desc%desc_1)
-        
       if (info /= 0) then 
         write(0,*) trim(name),' Error from inner routines',info
         info = -1
@@ -595,6 +617,8 @@ contains
   ! due to exch_fw_idx
   !
   subroutine psb_z_forward_map(alpha,x,beta,y,desc,info,work)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_comm_mod
     use psb_serial_mod
     use psi_mod
@@ -641,8 +665,8 @@ contains
       end if
       
     case(psb_map_gen_linear_)
-      call psb_halo(x,desc%desc_ext_1,info,work=work)
+      call psb_linmap(alpha,x,beta,y,desc%zmap%map_fw,&
-      if (info == 0) call psb_csmm(alpha,desc%zmap%map_fw,x,beta,y,info)
+           & desc%desc_fw,desc%desc_1,desc%desc_2)
         
       if (info /= 0) then 
         write(0,*) trim(name),' Error from inner routines',info
@@ -664,6 +688,8 @@ contains
   ! due to exch_bk_idx
   !
   subroutine psb_z_backward_map(alpha,x,beta,y,desc,info,work)
+    use psb_spmat_type
+    use psb_descriptor_type
     use psb_comm_mod
     use psb_serial_mod
     use psi_mod
@@ -711,8 +737,8 @@ contains
       
 
     case(psb_map_gen_linear_)
-      call psb_halo(x,desc%desc_ext_2,info,work=work)
+      call psb_linmap(alpha,x,beta,y,desc%zmap%map_bk,&
-      if (info == 0) call psb_csmm(alpha,desc%zmap%map_bk,x,beta,y,info)
+           & desc%desc_bk,desc%desc_1,desc%desc_2)
         
       if (info /= 0) then 
         write(0,*) trim(name),' Error from inner routines',info
@@ -728,4 +754,61 @@ contains
   end subroutine psb_z_backward_map
 
 
+
+  subroutine psb_d_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_comm_mod
+    use psb_serial_mod
+    use psi_mod
+    implicit none 
+    real(kind(1.d0)), intent(in)      :: alpha,beta
+    real(kind(1.d0)), intent(inout)   :: x(:),y(:)
+    type(psb_dspmat_type), intent(in) :: a_map
+    type(psb_desc_type), intent(in)   :: cd_xt,descin, descout 
+
+    integer :: nrt, nct, info
+    real(kind(1.d0)), allocatable :: tmp(:)
+
+    nrt = psb_cd_get_local_rows(cd_xt)
+    nct = psb_cd_get_local_cols(cd_xt)
+    allocate(tmp(nct),stat=info)
+    if (info == 0) tmp(1:nrt) = x(1:nrt)
+    if (info == 0) call psb_halo(tmp,cd_xt,info) 
+    if (info == 0) call psb_csmm(alpha,a_map,tmp,beta,y,info)
+    if (info /= 0) then 
+      write(0,*) 'Error in apply_map'
+    endif
+    
+  end subroutine psb_d_apply_linmap
+
+
+  subroutine psb_z_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_comm_mod
+    use psb_serial_mod
+    use psi_mod
+    implicit none 
+    complex(kind(1.d0)), intent(in)      :: alpha,beta
+    complex(kind(1.d0)), intent(inout)   :: x(:),y(:)
+    type(psb_zspmat_type), intent(in) :: a_map
+    type(psb_desc_type), intent(in)   :: cd_xt,descin, descout 
+
+    integer :: nrt, nct, info
+    complex(kind(1.d0)), allocatable :: tmp(:)
+
+    nrt = psb_cd_get_local_rows(cd_xt)
+    nct = psb_cd_get_local_cols(cd_xt)
+    allocate(tmp(nct),stat=info)
+    if (info == 0) tmp(1:nrt) = x(1:nrt)
+    if (info == 0) call psb_halo(tmp,cd_xt,info) 
+    if (info == 0) call psb_csmm(alpha,a_map,tmp,beta,y,info)
+    if (info /= 0) then 
+      write(0,*) 'Error in apply_map'
+    endif
+    
+  end subroutine psb_z_apply_linmap
+
+
 end module psb_inter_descriptor_type

base/modules/psb_tools_mod.f90

@@ -450,6 +450,16 @@ Module psb_tools_mod
       integer, intent(out)                 :: info
       logical, intent(in), optional        :: rebuild
     end subroutine psb_zspins
+    subroutine psb_zspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info)
+      use psb_descriptor_type
+      use psb_spmat_type
+      type(psb_desc_type), intent(in)      :: desc_ar
+      type(psb_desc_type), intent(inout)   :: desc_ac
+      type(psb_zspmat_type), intent(inout) :: a
+      integer, intent(in)                  :: nz,ia(:),ja(:)
+      complex(kind(1.d0)), intent(in)      :: val(:)
+      integer, intent(out)                 :: info
+    end subroutine psb_zspins_2desc
   end interface
 
 
@@ -564,6 +574,17 @@ Module psb_tools_mod
     end subroutine psb_get_ovrlap
   end interface
 
+  interface psb_linmap_init
+    module procedure psb_dlinmap_init, psb_zlinmap_init
+  end interface
+
+  interface psb_linmap_ins
+    module procedure psb_dlinmap_ins, psb_zlinmap_ins
+  end interface
+
+  interface psb_linmap_asb
+    module procedure psb_dlinmap_asb, psb_zlinmap_asb
+  end interface
 
 
 contains
@@ -738,5 +759,138 @@ contains
   end subroutine psb_cdasb
 
 
+  subroutine psb_dlinmap_init(a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_serial_mod
+    use psb_penv_mod
+    use psb_error_mod
+    implicit none 
+    type(psb_dspmat_type), intent(out) :: a_map
+    type(psb_desc_type), intent(out)   :: cd_xt
+    type(psb_desc_type), intent(in)    :: descin, descout 
+
+    integer :: nrow_in, nrow_out, ncol_in, info, ictxt
+
+    ictxt = psb_cd_get_context(descin)
+    call psb_cdcpy(descin,cd_xt,info)
+    if (info ==0) call psb_cd_reinit(cd_xt,info)
+    if (info /= 0) then 
+      write(0,*) 'Error on reinitialising the extension map'
+      call psb_error(ictxt)
+      call psb_abort(ictxt)
+      stop
+    end if
+
+    nrow_in  = psb_cd_get_local_rows(cd_xt)
+    ncol_in  = psb_cd_get_local_cols(cd_xt)
+    nrow_out = psb_cd_get_local_rows(descout)
+
+    call psb_sp_all(nrow_out,ncol_in,a_map,info)
+
+  end subroutine psb_dlinmap_init
+
+  subroutine psb_dlinmap_ins(nz,ir,ic,val,a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    implicit none 
+    integer, intent(in)                  :: nz
+    integer, intent(in)                  :: ir(:),ic(:)
+    real(kind(1.d0)), intent(in)         :: val(:)
+    type(psb_dspmat_type), intent(inout) :: a_map
+    type(psb_desc_type), intent(inout)   :: cd_xt
+    type(psb_desc_type), intent(in)      :: descin, descout 
+    integer :: info
+    call psb_spins(nz,ir,ic,val,a_map,descout,cd_xt,info)
+
+  end subroutine psb_dlinmap_ins
+
+  subroutine psb_dlinmap_asb(a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_serial_mod
+    implicit none 
+    type(psb_dspmat_type), intent(inout) :: a_map
+    type(psb_desc_type), intent(inout)   :: cd_xt
+    type(psb_desc_type), intent(in)      :: descin, descout 
+
+    integer :: nrow_in, nrow_out, ncol_in, info, ictxt
+
+    ictxt = psb_cd_get_context(descin)
+
+    call psb_cdasb(cd_xt,info)
+    a_map%k = psb_cd_get_local_cols(cd_xt)
+    call psb_spcnv(a_map,info,afmt='CSR')
+
+  end subroutine psb_dlinmap_asb
+
+  subroutine psb_zlinmap_init(a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_serial_mod
+    use psb_penv_mod
+    use psb_error_mod
+    implicit none 
+    type(psb_zspmat_type), intent(out) :: a_map
+    type(psb_desc_type), intent(out)   :: cd_xt
+    type(psb_desc_type), intent(in)    :: descin, descout 
+
+    integer :: nrow_in, nrow_out, ncol_in, info, ictxt
+
+    ictxt = psb_cd_get_context(descin)
+
+    call psb_cdcpy(descin,cd_xt,info)
+    if (info ==0) call psb_cd_reinit(cd_xt,info)
+    if (info /= 0) then 
+      write(0,*) 'Error on reinitialising the extension map'
+      call psb_error(ictxt)
+      call psb_abort(ictxt)
+      stop
+    end if
+
+    nrow_in  = psb_cd_get_local_rows(cd_xt)
+    ncol_in  = psb_cd_get_local_cols(cd_xt)
+    nrow_out = psb_cd_get_local_rows(descout)
+
+    call psb_sp_all(nrow_out,ncol_in,a_map,info)
+
+  end subroutine psb_zlinmap_init
+
+  subroutine psb_zlinmap_ins(nz,ir,ic,val,a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    implicit none 
+    integer, intent(in)                  :: nz
+    integer, intent(in)                  :: ir(:),ic(:)
+    complex(kind(1.d0)), intent(in)      :: val(:)
+    type(psb_zspmat_type), intent(inout) :: a_map
+    type(psb_desc_type), intent(inout)   :: cd_xt
+    type(psb_desc_type), intent(in)      :: descin, descout 
+    integer :: info
+
+    call psb_spins(nz,ir,ic,val,a_map,descout,cd_xt,info)
+
+  end subroutine psb_zlinmap_ins
+
+  subroutine psb_zlinmap_asb(a_map,cd_xt,descin,descout)
+    use psb_spmat_type
+    use psb_descriptor_type
+    use psb_serial_mod
+    implicit none 
+    type(psb_zspmat_type), intent(inout) :: a_map
+    type(psb_desc_type), intent(inout)   :: cd_xt
+    type(psb_desc_type), intent(in)      :: descin, descout 
+
+    integer :: nrow_in, nrow_out, ncol_in, info, ictxt
+
+    ictxt = psb_cd_get_context(descin)
+
+    call psb_cdasb(cd_xt,info)
+    a_map%k = psb_cd_get_local_cols(cd_xt)
+    call psb_spcnv(a_map,info,afmt='CSR')
+
+  end subroutine psb_zlinmap_asb
+
+
 
 end module psb_tools_mod

base/tools/psb_zspins.f90

@@ -240,3 +240,164 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
 
 end subroutine psb_zspins
 
+
+subroutine psb_zspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info)
+  use psb_tools_mod, psb_protect_name => psb_zspins_2desc
+  use psb_descriptor_type
+  use psb_spmat_type
+  use psb_serial_mod
+  use psb_const_mod
+  use psb_error_mod
+  use psb_penv_mod
+  implicit none
+
+  !....parameters...
+  type(psb_desc_type), intent(in)      :: desc_ar
+  type(psb_desc_type), intent(inout)   :: desc_ac
+  type(psb_zspmat_type), intent(inout) :: a
+  integer, intent(in)                  :: nz,ia(:),ja(:)
+  complex(kind(1.d0)), intent(in)         :: val(:)
+  integer, intent(out)                 :: info
+  !locals.....
+
+  integer :: nrow, err_act, ncol, spstate
+  integer                :: ictxt,np,me
+  logical, parameter     :: debug=.false.
+  integer, parameter     :: relocsz=200
+  integer, allocatable   :: ila(:),jla(:)
+  character(len=20)  :: name, ch_err
+
+  info = 0
+  name = 'psb_zspins'
+  call psb_erractionsave(err_act)
+
+
+  ictxt = psb_cd_get_context(desc_ar)
+
+  call psb_info(ictxt, me, np)
+
+  if (.not.psb_is_ok_desc(desc_ar)) then 
+    info = 3110
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  if (.not.psb_is_ok_desc(desc_ac)) then 
+    info = 3110
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  if (nz < 0) then 
+    info = 1111
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+  if (size(ia) < nz) then 
+    info = 1111
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  if (size(ja) < nz) then 
+    info = 1111
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+  if (size(val) < nz) then 
+    info = 1111
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+  if (nz==0) return
+
+  spstate = a%infoa(psb_state_)
+  if (psb_is_bld_desc(desc_ac)) then 
+
+    allocate(ila(nz),jla(nz),stat=info)
+    if (info /= 0) then
+      ch_err='allocate'
+      call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
+      goto 9999
+    end if
+        ila(1:nz) = ia(1:nz)
+
+    call psb_glob_to_loc(ia(1:nz),ila(1:nz),desc_ar,info,iact='I',owned=.true.)
+
+    call psb_cdins(nz,ja,desc_ac,info,jla=jla, mask=(ila(1:nz)>0))
+
+    if (info /= 0) then
+      ch_err='psb_cdins'
+      call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
+      goto 9999
+    end if
+
+    nrow = psb_cd_get_local_rows(desc_ar)
+    ncol = psb_cd_get_local_cols(desc_ac)
+
+    call psb_coins(nz,ila,jla,val,a,1,nrow,1,ncol,info)
+    if (info /= 0) then
+      info=4010
+      ch_err='psb_coins'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+  else if (psb_is_asb_desc(desc_ac)) then 
+
+    write(0,*) 'Why are you calling me on an assembled desc_ac?'
+!!$    if (psb_is_large_desc(desc_a)) then 
+!!$
+!!$      allocate(ila(nz),jla(nz),stat=info)
+!!$      if (info /= 0) then
+!!$        ch_err='allocate'
+!!$        call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
+!!$        goto 9999
+!!$      end if
+!!$
+!!$      ila(1:nz) = ia(1:nz)
+!!$      jla(1:nz) = ja(1:nz)
+!!$      call psb_glob_to_loc(ila(1:nz),desc_a,info,iact='I')
+!!$      call psb_glob_to_loc(jla(1:nz),desc_a,info,iact='I')
+!!$      nrow = psb_cd_get_local_rows(desc_a)
+!!$      ncol = psb_cd_get_local_cols(desc_a)
+!!$
+!!$      call psb_coins(nz,ila,jla,val,a,1,nrow,1,ncol,&
+!!$           & info,rebuild=rebuild_)
+!!$      if (info /= 0) then
+!!$        info=4010
+!!$        ch_err='psb_coins'
+!!$        call psb_errpush(info,name,a_err=ch_err)
+!!$        goto 9999
+!!$      end if
+!!$
+!!$    else
+!!$      nrow = psb_cd_get_local_rows(desc_a)
+!!$      ncol = psb_cd_get_local_cols(desc_a)
+!!$      call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,&
+!!$           & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_)
+!!$      if (info /= 0) then
+!!$        info=4010
+!!$        ch_err='psb_coins'
+!!$        call psb_errpush(info,name,a_err=ch_err)
+!!$        goto 9999
+!!$      end if
+!!$    end if
+  else
+    info = 1122
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act == psb_act_abort_) then
+    call psb_error(ictxt)
+    return
+  end if
+  return
+
+end subroutine psb_zspins_2desc
+