diff --git a/base/modules/psb_inter_desc_type.f90 b/base/modules/psb_inter_desc_type.f90 index f7f85465..2fd5dcea 100644 --- a/base/modules/psb_inter_desc_type.f90 +++ b/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_descriptor_type + use psb_spmat_type, only : psb_dspmat_type, psb_zspmat_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) - - implicit none + use psb_spmat_type + 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_ext_2) + val = val + psb_sizeof(desc%desc_fw) + 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) - if (info == 0) call psb_csmm(alpha,desc%dmap%map_bk,x,beta,y,info) - + call psb_linmap(alpha,x,beta,y,desc%dmap%map_bk,& + & 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) - if (info == 0) call psb_csmm(alpha,desc%zmap%map_fw,x,beta,y,info) + call psb_linmap(alpha,x,beta,y,desc%zmap%map_fw,& + & 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) - if (info == 0) call psb_csmm(alpha,desc%zmap%map_bk,x,beta,y,info) + call psb_linmap(alpha,x,beta,y,desc%zmap%map_bk,& + & 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 diff --git a/base/modules/psb_tools_mod.f90 b/base/modules/psb_tools_mod.f90 index 48a3ef3e..06ec9a12 100644 --- a/base/modules/psb_tools_mod.f90 +++ b/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 diff --git a/base/tools/psb_zspins.f90 b/base/tools/psb_zspins.f90 index aa7bfc8b..dc700530 100644 --- a/base/tools/psb_zspins.f90 +++ b/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 +